Trait-based model development to support breeding programs. A case study for salt tolerance and rice

Eco-physiological models are increasingly used to analyze G × E × M interactions to support breeding programs via the design of ideotypes for specific contexts. However, available crop models are only partly suitable for this purpose, since they often lack clear relationships between parameters and traits breeders are working on. Taking salt stress tolerance and rice as a case study, we propose a paradigm shift towards the building of ideotyping-specific models explicitly around traits involved in breeding programs. Salt tolerance is a complex trait relying on different physiological processes that can be alternatively selected to improve the overall crop tolerance. We developed a new model explicitly accounting for these traits and we evaluated its performance using data from growth chamber experiments (e.g., R2 ranged from 0.74 to 0.94 for the biomass of different plant organs). Using the model, we were able to show how an increase in the overall tolerance can derive from completely different physiological mechanisms according to soil/water salinity dynamics. The study demonstrated that a trait-based approach can increase the usefulness of mathematical models for supporting breeding programs.

dependencies play a key role in determining the optimal value of traits contributing to complex plant responses, like in the case of tolerance to abiotic stressors 11 .
With more than 830 million hectares of salt-affected soils, salinity is one of the major environmental stress limiting agricultural production worldwide 12 . Moreover, soil salinization is further increasing 13 because of both human activities (e.g., inappropriate irrigation practices) and natural causes (e.g., tsunamis), with the latter being exacerbated by climate change 14 . Despite rice is one of the most sensitive crop to salt stress 15,16 , its frequent cultivation in coastal areas and river deltas increases it exposure because of recurrent flooding and seawater intrusion. For these reasons, ongoing rice breeding activities target different tolerance traits: (i) reduction of shoot sodium (Na + ) uptake 17,18 ; (ii) sequestration of Na + into structural tissue 19,20 ; (iii) compartmentation of Na + into senescent leaves 21 , (iv) higher leaf tissue tolerance via sequestration of toxic ions into the vacuole and synthesis of osmoprotectants 22 , (v) and higher tolerance to salt-induced sterility 23 . These traits rely on different genetic basis 15 and can be combined in the same genotype or singularly introduced in commercial varieties even via non-GM technologies (i.e., using marker-assisted selection) 24 . Since the effectiveness of these traits varies according to the environmental context, breeding programs should target these "component traits" and not the overall performance at crop level 25,26 .
Although models reproducing crop response to salt stress are available 27,28 , they mainly focus on the effect of salinity on soil osmotic potential, without explicitly considering the toxic effect of Na + in plant tissues which, instead, is a key component of salt stress 15,17 . In these approaches, plant tolerance is accounted for via few empirical parameters directly linking plant response (in terms of yield or overall growth rate) to soil salinity, without taking into account the physiological traits at the basis of such response. Therefore, these models cannot be considered as suitable to design ideotypes actually relying on the real tolerance traits identified by breeders.
The objectives of this study were (i) building a new model for the toxic effect of Na + on rice by explicitly taking into account the tolerance traits breeders are working on, and (ii) presenting a case study on ideotype design targeting production districts in California and Greece.

Results and Discussion
A new model for salt stress built around actual plant traits. A new model for the salt stress on rice was developed by directly targeting tolerance traits breeders are working on (Fig. 1). In the following equations, terms with the prefix (T-) refer to traits (Table 1).
Shoot Na + uptake. The plant capability to reduce the shoot Na + uptake relies on a "root filter" 29 that prevents Na + from entering the roots and getting translocated, via the xylematic stream, to the shoot 26 . The root filter can be more or less pronounced, leading to the identification of "excluder" genotypes as opposite to varieties less effective in preventing Na + from entering the xylem stream 17 .
From a physiological standpoint, the root filter is made of two components: morphological barriers reducing the apoplastic entry of Na + (bypass flow), and channels in the plasma membrane of root epidermal/cortical cells that mediate the selective and non-selective transport of Na + 18, 30 . Although by-pass flow accounts for 22-35% of the total Na + uptake in rice 17,31 , it is a key responsible for the differences in the degree of salt tolerance in rice genotypes. This makes the reduction of bypass flow a promising trait for increasing salt tolerance in rice 17,32 . For this reason and despite its relative importance compared to the other pathway, we decided to focus on by-pass flow to model genotypic differences, given its prominent role in rice breeding activities 17 . The amount of water potentially flowing through the apoplastic pathway (Jv B , mm day −1 ) is alculated according to the following equation: where Tr act (mm day −1 ) is the actual transpiration and RJv B (%) is the percentage of water uptake through the apoplastic pathway, estimated as a function of air relative humidity (RH; %) according to Faiyue et al. 31 and Steudle et al. 33 : The development of Casparian bands and the deposition of suberin lamellae in the root exo-and endodermis reduce the water transport through the apoplast and thus the bypass flow-Na + uptake 34 . The relative bypass flow (RRBF, unitless, 0-1) is thus calculated as a function of the root suberin content 35 : where (T1)RRBF min (%) is the bypass flow when the suberin content is maximum; SC (mg g −1 ) is the root suberin content; (T1)SCmax (mg g −1 ) is the maximum suberin content; SC min (mg g −1 ) is the minimum root suberin concentration at which bypass flow starts to be reduced. The root suberin content is derived as a function of plant age (equation 4) and of the genotype sensitivity to salinity (equation 5) 34, 36 : where DVS (unitless; 0-2) is a SUCROS-type development stage code (0: emergence; 1: anthesis; 2: maturity); F sc (unitless; 0-1) is deposition of suberin in response to salinity; (T1)SubDepEff (unitless; 0-1) is the suberin deposition efficiency (genotype specific); [Na + ] ext and Max[Na + ] ext (mM) are, respectively, the actual and maximum (at which the suberin deposition is maximum) Na + concentrations in the external medium. Therefore, the amount of Na + actually delivered to the shoot via bypass-flow (NaUptake AP , mg ha −1 ) is: The fraction of bypass Na + on the total Na + uptake (RNaUptake AP ; %) is derived analogously to equation (2) 31 . Finally, the total amount of Na + daily entering the shoot (NaUptake, mg ha −1 ) is derived according to equation (7): = NaUptake NaUptake RNaUptake 100 (7) AP AP Sequestration of Na + into structural/senescent organs. To reduce the amount of Na + reaching the leaves, plants have developed mechanisms to accumulate toxic ions in the tissues of culm base and leaf sheath 19,20 . The former is estimated considering genotypic differences and feedback mechanisms triggered by the amount of Na + already stored in culms 29 . The actual culm sequestration rate (ActCulmSeqRate; mg ha −1 ) is derived as: 1 13 Figure 1. Flowchart of the new model for the impact of salt stress on rice growth. T1, T2, T3, T4, and T5 refer, respectively, to the traits involved with the capability of limiting shoot Na + uptake, sequestrating Na + into the culm base, reducing salt-induced spikelet sterility, storing toxic ions in senescent leaves, and decreasing the toxicity of Na + in photosynthetically active leaves. Grey items represent coupling points between the salinity model and the crop simulator. White items represent variables and parameters of the salinity model. Black items represent parameters corresponding to traits directly involved with tolerance to salinity. Symbols meaning: represents the main driving variable for salinity-related processes; identifies parameters; and refer to state and rate variables, respectively; indicates a process.
The amount of Na + daily reaching the panicle (NaPanicle, mg ha −1 ) is derived as a function of the Na + not sequestrated in culms and of the transport of photosynthates to panicles:

NaPanicle
ParP T Na ToPan NaUptake ActCulmSeqRate DVS elsewhere where ParP (unitless) is the relative amount of photosynthates daily partitioned to panicles and (T3)Na + ToPan (unitless; 0-1) is the factor for Na + translocation to panicles. The amount of Na + daily reaching leaves (Na leaves , mg ha −1 ) is the difference between total Na + uptake in shoots and the amounts of Na + sequestrated in the culms and partitioned to panicles. Plants tend to accumulate Na + in the oldest leaves to preserve photosynthetically active tissues from toxic ions 19,30 . This is represented using equation (11): where NaUptakeL i (mg ha −1 ) is the Na + delivered to the ith canopy layer (starting to count from the bottom); (T4)PartCap (unitless; 0-1) is the genotype capability of partitioning Na + to oldest leaves; N is the number of living canopy layers. To account for this heterogeneity in Na + accumulation in leaves of different ages, the model for salt stress should be coupled with a crop model providing a multilayer canopy structure.
Impact of Na + on photosynthesis, leaf senescence and spikelet sterility. The Na + effect on photosynthesis 37 and leaf senescence 38, 39 depends on the genotype ability to sequestrate toxic ions in the vacuole and to synthesize osmolytes to counterbalance the osmotic pressure. The stress factor for photosynthesis (RPn, unitless 0-1) -also used to increase senescence -is derived using equation (12): is the Na + concentration in leaves at the ith canopy layer; [Na + ] leaf min (mg g −1 ) is the Na + concentration in unstressed leaves; (T5)ThreshL (mg g −1 ) is the Na + concentration above which salt stress starts; (T5)CritL (mg g −1 ) is the Na + concentration at which photosynthesis becomes null; C is a shaping coefficient. Na + also increases maintenance respiration due to the high metabolic costs of the processes of ion exclusion, vacuolar compartmentation and synthesis of osmolytes 40 . The factor increasing maintenance respiration in leaves (MRespF, unitless, 0-1) is thus: leaves Thresh

Crit Thresh
where [Na + ] Crit (set to 3 mg g −1 ) is the Na + concentration at which respiration is double; [Na + ] Thresh (set to 0.5 mg g −1 ) is the Na + concentration at which maintenance respiration starts to be affected; [Na + ] leaves (mg g −1 ) is the average concentration in leaves (weighted for layers' biomass). The same function is used for the increased maintenance respiration in culms. Salt stress also affects spikelet sterility 23 according to panicle Na + concentration and plant susceptibility, the latter depending on the genotype and phenological stage (equation 14).
where SterilityF (unitless, 0-1) is the factor reducing spikelet fertility due to salt stress; (T3)SuscSt (unitless; 0-1) represents the genotype susceptibility; bellF (unitless; 0-1) is a factor modulating susceptibility according to the within-and between-plant heterogeneity phenological development 41 . The bellF is calculated considering two phenological stages of maximum susceptibility to abiotic stress-induced sterility: booting (microsporogenesis; DVS = 0.8) and flowering (DVS = 1). SterilityF is then used to reduce the amount of photosynthates daily partitioned to panicles. The reduction of growth due to the osmotic potential in the external medium 15  where SLAstress (unitless, 0-1) is the reduction factor for specific leaf area (SLA; m 2 kg −1 ); CVSstress (unitless, 0-1) is the reduction factor for culm growth; is the Na + concentration in the external medium. The salinity model was coupled to the WOFOST model as modified by Stella et al. 44 , to benefit from an explicit multi-layer canopy representation. This version of WOFOST, like the other versions available, does not include algorithms for salt stress. The resulting modelling solution was evaluated using the growth chamber datasets described below and then used for the ideotyping study.
Model evaluation. The agreement between observed and simulated values of aboveground biomass, yield, biomass of culms, leaves and panicles, leaf area index and plant sodium content is shown in Fig. 2. In general, the model showed good performances in reproducing the impact of salt stress on aboveground biomass accumulation and yield (Fig. 2a), with relative root mean square error (RRMSE; %; 0 to +∞, optimum 0) equal to 28.0% and 23.8%, respectively. Good values for these two variables were achieved also for R 2 (0.89 and 0.90) and modelling efficiency (EF; −∞ to +1, optimum +1) 45  Good performances were achieved also for the simulation of the biomass of culms, leaves, and panicles (Fig. 2b), with RRMSE never exceeding 36% and values of R 2 ranging between 0.74 (leaves) and 0.94 (panicles). Regardless of the organ, EF was always positive. Concerning the simulation of panicles weight, the model showed a slight tendency to underestimate the values at heading, likely because of the observed heterogeneity in tiller development. An opposite behavior (slight overestimation around heading) was instead observed for culm and leaf biomass. In this case, the reason is likely an underestimation of the impact of the osmotic component of salt stress on leaf area expansion and tiller development, which are only implicitly reproduced in the current version of the modelling solution, since the crop model used does not explicitly simulate tillering and leaf size. Future model improvements could thus refer to the implementation of dedicated approaches for the simulation of the Na + -induced reduction in cell turgor pressure and related decrease in tissue expansion.
The discussed overestimation of leaf biomass explains the similar behavior showed by the model for the simulation of leaf area index (Fig. 2c) (CRM = −0.33), although the model correctly reproduced the relative reduction in leaf area for increasing Na + concentrations (R 2 = 0.89; EF = 0.65). This is considered as particularly important, since leaf area index is a key state variable involved with the amount of water daily transpired and thus with the potential entry of Na + through the apoplastic pathway 17 (equations 1-6).
Scientific RepoRts | 7: 4352 | DOI:10.1038/s41598-017-04022-y Concerning Na + uptake (Fig. 2d), simulated Na + contents showed a good agreement with measured data, although the values of the performance metrics were slightly worse compared to those achieved for other outputs (RRMSE was around 50% and CRM was −0.20). However, the large portion of variance explained (R 2 = 0.72) and the largely positive value of EF (0.66) allow considering also this variable as satisfactorily simulated. Indeed, Na + content results from the simulation of both Na + dynamics (entry, translocation, sequestration, etc.) and the related effects on crop growth, and also by the uncertainty of the crop model itself. The interaction between the crop model and the salt stress model is deeply involved with the way plant growth drives Na + uptake through the simulation of actual transpiration (in turn driven by leaf area index) and Na + sink capacity (driven by organ biomass). For these reasons, the general overestimation of plant Na + content could be related to the overestimation in leaf area index and culm biomass ( Fig. 2c and b). However, the model was not able to reproduce the high plant Na + content measured at harvest for the cultivar Baldo at the highest Na level (150 mg plant −1 ). Another aspect able to partly explain the values for the agreement metrics obtained for the simulation of Na + content is the Na + accumulation for the control treatment (0 mM NaCl). In this case, the model simulates a null Na + plant content, whereas small values were measured in real plant tissues.

Identification of key traits in different scenarios.
Simulations performed during the SA (5632 for each site) led to 14% and 12% average yield losses in California and Greece, respectively, in agreement with the expected yield reduction under salinity levels similar to those explored 47  Results of the SA (Fig. 3) showed how different scenarios could allow defining contrasting breeding strategies to increase the overall cultivar salt tolerance and, in turn, further demonstrated the potential of trait-based modelling for designing district-specific ideotypes 6 . Indeed, simulations revealed that -despite salt tolerance is an issue in both scenarios -different traits would guarantee the highest increase in yields in California and Greece. This is in agreement with Roy et al. 26 , who observed how different traits could be exploited to increase salt tolerance under different salinity levels. Tissue tolerance (T5, indicated in blue in Fig. 3) was the most important trait in California (Fig. 3b), where high-salinity peaks occur for short periods in the first part of the season (Fig. 4). Indeed, although the reduction of shoot Na + uptake (T1, green in Fig. 3) also played a role (because of the relevance of a parameter involved with Na + exclusion at root level via suberin deposition), the sharp increase in salinity in a moment when root barriers are still not developed makes the response at leaf level more important to increase the overall plant tolerance. A similar peak occurring later in the season (with a higher root suberin content) would have led to different results.
In case of prolonged stressful conditions like those characterizing the Greek scenario (Fig. 4), instead, the most important tolerance traits were involved with the plant capability to prevent Na + from reaching leaf blades (Fig. 3d), i.e., reduction of Na + uptake (suberin deposition) and − to a lesser extent -higher sequestration in structural tissues (T2, indicated in red in Fig. 3). The trait involved with tissue tolerance was not considered as relevant for this scenario, because Na + accumulation in leaf tissues was so fast to rapidly overcome the capability of the plant to segregate toxic ions into the vacuole and to synthesize osmolytes. This is in agreement with Munns et al. 24 , who observed the same relationship between high salinity and the importance of excluding Na + at root level for wheat.
The comparison of the Sobol' first-and total-order effects highlighted strong interactions only for the Californian scenario ( Fig. 3a and b) for the two parameters modulating the toxic effect of Na + in leaf blades, i.e., (T5)ThreshL and (T5)CritL. This is partly due to the fact that these two parameters are involved in the same response function (equation 12), thus they interact with each other. However, the large differences in the values of the first-and total-order sensitivity metrics suggest that these two parameters interacted also with others, being tissue tolerance strictly depending on the whole chain of processes modulating the sodium uptake at root level and its translocation and accumulation in photosynthetic tissues.  Table 1). The Sobol' sensitivity metrics for a parameter indicate the portion of output variance explained by changes in the value of that parameter alone (first-order effect) and of that parameter in interaction with all the others (total-order effect). The analysis was carried out for two scenarios, one in California (a,b) and the other in Greece (c,d). The output for which sensitivity metrics were calculated was the final yield. Regardless of the scenario, parameters involved with tolerance to salt-induced spikelet sterility did not result important in affecting yields. The reason is that Na + concentration in the young panicles was not high enough during the two sensitive stages (i.e., booting and flowering). In California, indeed, the accumulation of Na + in the developing panicles before flowering was low because of the reduction of salinity due to fresh water inflow after the herbicide treatments, whereas in Greece Na + concentration in field water started to increase too late to generate high contents in the reproductive organs before flowering.
For both California and Greece, the parameter values for the ideotypes were those that led to achieve the highest yield (i.e., those leading to the highest overall tolerance) during the SA experiments. Parameter values for the two ideotypes are shown in Table 2, together with the parameter values for the two cultivars used in the growth chamber experiments to evaluate the salinity model.

Concluding remarks.
Taking salt stress tolerance and rice as a case study, we showed how the design of ideotypes can benefit from the availability of models explicitly developed starting from traits breeders are working on. Indeed, the model developed demonstrated its suitability for analyzing in depth key G × E × M interactions. Salt-induced rice yield losses depend on the extent at which the genotype tolerance traits are effective in dealing with the salt stress dynamics deriving from management practices (e.g., water management modulating salinity level of rice field water) and weather conditions (e.g., evapotranspiration affecting potential sodium entry via by-pass flow). As our results show, these G × E × M interactions would lead the same genotype being tolerant in one environment and susceptible in another, as the ideotype defined for California, which would fail in coping with the prolonged salt-stress period of the Greek scenario.
However, salt tolerance is a complex trait relaying on a variety of physiological mechanisms that can be alternatively selected to improve the overall crop tolerance. The availability of new tools for genetic improvement (e.g., Marker Assisted Selection) allows breeding for specific traits instead of targeting the overall crop tolerance 26 , thus increasing the efficiency of the breeding process, since -as explained above -the effectiveness of the changes in the values of different traits varies according to the agro-environmental context. For the first time, a modelling approach dedicated to ideotyping was developed by explicitly building algorithms around traits for which breeding activities are ongoing. We consider this strategy as the only one able to avoid inconsistencies between model parameters and plant traits, and thus between in silico ideotypes and the possibility of realizing them in vivo. Indeed, a high level of detail in the representation of physiological processes cannot be considered as a guarantee of direct relationships between model parameters and plant traits, given the same knowledge can be formalized in a variety of possible modelling structure 48 . Moreover, during model development, a pronounced process-based perspective was used to properly account for the key physiological processes and feedback mechanisms involved.
Results of model application to design district-specific ideotypes showed that, despite differences between scenarios were mainly limited to the seasonal dynamics of salt concentration in field water, putative traits to increase salt tolerance can rely on completely different physiological mechanisms. Results achieved thus encourage a paradigm shift towards the development of dedicated trait-based models and their use for supporting breeding programs at district level.
Limits of our study − and thus potential areas for model improvement − deal with the lack of approaches to simulate also the uptake and distribution of K + , since the toxicity of Na + seems to be related also to the Na + : K + ratio and not only to the concentration of Na + in itself 18 . Also, improved (i.e., more explicit) approaches would be needed for the simulation of the effect of the osmotic component of salt stress on tissue expansion and stomatal reaction, in order to better account for osmotic adjustment.

Methods
The growth chamber experiments. Two rice (Oryza sativa L. spp. japonica) cultivars differing in their level of salt tolerance, i.e., Baldo and Vialone Nano, were grown in dedicated hydroponics growth chamber Si (pH 6.5). Floating polystyrene foam sheets were used to hold seedlings and to allow renewing solutions without touching the plants to avoid any potential damage to the roots (which would make Na + entering directly from root lesions). Growing conditions were set to 14 h photoperiod; photosynthetically active radiation was supplied by fluorescent lamps (400 mmol m −2 s −1 ); day/night temperatures were 26 °C/18 °C; relative humidity ranged between 58 and 92%. Five NaCl treatments were applied from three weeks after sowing until maturity: 0 mM, 10 mM, 25 mM, 35 mM and 50 mM. In order to maintain NaCl concentrations nearly constant, solutions were renewed each three days. At late heading (BBCH code 51) and maturity (BBCH code 92) three plants for each combination cultivar × treatment were harvested and divided into stems, panicles and leaves. The latter were further separated into apical, medium and senescent leaves (referring respectively to the two youngest leaves, other green leaves and dead ones) to detect potential variation in Na + accumulation among leaves of different ages. Plant height, number of tillers and dry biomass of each organ were measured. Dry biomass samples were ground to a fine powder and digested by concentrated HNO 3 (10 mM) in a microwave digester (ETHOS D, milestone, Italy) at 100 °C 50 . The mineralized material was dissolved in 5 mL 0.1 M HNO 3 and Na + content was measured by inductively coupled plasma mass spectrometry (Bruker Aurora M90 ICP-MS, ICP Mass Spectrometer). Na + content and dry biomass of different plant organs were then used to calculate corresponding Na + concentrations. Immediately before the first sampling event (late heading), the impact of salt stress on net photosynthetic rate, stomatal conductance and transpiration rate was measured on the youngest fully expanded leaf using a CIRAS-3 Portable Photosynthesis System (PP Systems, Amsbury, MA, USA). Apical and medium leaves were then scanned to determine plant leaf area and SLA, the latter calculated as ratio between leaf area and leaf dry biomass. At harvest, spikelet sterility was determined.
The ideotyping study. Case studies and simulation scenarios. An ideotyping study was carried out to using the salt stress model to identify the traits breeders should focus on in two production areas differing for the salinity seasonal dynamics (Fig. 4). Colusa is one of the six counties at the north of Sacramento where the production of rice in Californiathe second largest rice-producing state in the USis concentrated. While most of the irrigation water has a low salt content, water holding periods for herbicide distribution and high temperatures promote evapo-concentration of salt in rice fields 51, 52 . This leads salinity in field water during the first part of the crop cycle to increase up to 3.5 dS m −1 , and to decrease rapidly once the flow of fresh water is restored. Although the mean seasonal salinity is not high, the rice susceptibility during early phenological phases can lead yield losses to exceed 10% 47 . The second scenario targets the southeastern region of Axios River plain, near Thessaloniki, one of the key regions for rice production in Greece 53 . Ninety percent of the soils in the area are saline, causing increases in the salts content of irrigation water during infiltration. Evapo-concentration of salts due to high temperatures also contributes to increase salinity, which progressively reaches values of 2.5 dS m −1 from mid-season to harvest. Dynamics of salinity in field water and information on management practices used for the two case studies were derived from Scardaci et al. 51 and Linquist et al. 52 (California), and from Lekakis et al. 54 (Greece).
Simulations were performed for the cultivar Thaibonnet (also known as L202), developed by the California Co-operate Rice Research Foundation and currently representing an important variety in Greece 53 .
Sensitivity analysis. The ideotyping study was performed using global sensitivity analysis (SA) techniques 4, 5 . In particular, the variance-based method of Sobol' 55 -considered as a references for SA 56 -was used, targeting yield as reference output. The analysis focused on first-and total-order effects, accounting, respectively, for the effects of variations in each parameter on simulated yield, and for the effects of variations in parameter including possible interactions with others. The sample size for the combinations of parameters was 5632, i.e., the lowest value of M | M > (γ · n), with M = 2 q+3 (2n + 2), q = {1, 2, 3, …, Q}, γ is the suggested number of model runs for each parameter, and n is the number of parameters in the sensitivity analysis. In this study, γ was set to 500 57 .
Parameterization of the crop model WOFOST-GT2 for Tropical Japonica rice varieties was derived from Stella et al. 44 . Concerning the parameters of the salt stress model, the values derived from the growth chamber experiments were used as means (Table 1), with the exception of a correction factor applied to the maximum suberin content to account for differences in root development between hydroponic and soil conditions 58 . Distributions for the SA were assumed as normal, and standard deviations were set to 5% of the mean values for the parameters 59 .
In order to avoid the risk of including in SA results the effect of a specific meteorological season, simulations for both sites were performed on 25-year series of weather data (1990-2015) retrieved from the European Centre for Medium-Range Weather Forecasts (ECMWF; ERA-Interim database; www.ecmwf.int).