In silico system analysis of physiological traits determining grain yield and protein concentration for wheat as influenced by climate and crop management

Highlight A global uncertainty and sensitivity analysis of the process-based model SiriusQuality2 was carried out to quantify the relationship between simple morpho-physiological traits and grain yield and protein concentration for wheat.


Introduction
It has been estimated that in order to ensure food security, crop grain yields (GYs) should be increased globally by 70-100% within the next 40 years (Bruinsma, 2009). This means that the past relative rate of world GY increase, estimated at 0.5-1.7% for the major crops (Cassman, 2001;Tilman et al., 2002;Ewert et al., 2005), needs to be increased by ~40% (Tester and Langridge, 2010). However, in contrast to past GY increases, future GY improvements will have to be achieved in the context of global warming and with reduction of the use of water and fertilizers, because of environmental issues (Tilman, 1999;Godfray et al., 2010). At the same time, grain protein concentration (GPC) will have to be maintained at its current level, since grain proteins are the major source of dietary proteins for humans, especially in developing countries, and since the economic value of most grain crops greatly depends on their GPC (Shewry, 2007). Maintaining GPC while increasing GY represents another challenge for plant breeders because of the commonly observed negative correlation between GY and GPC (e.g. Cooper et al., 2001;Oury et al., 2003;Aguirrezábal et al., 2014).
Genetic improvements of GY and GPC are also impeded by their low heritability and by large genotype×environment×management interactions (Cooper et al., 2001;Zheng et al., 2009). Considerable efforts have been made to identify environmentally stable and genetically heritable traits and processes related to GY that can be used to guide crop breeding programmes, but only a few successful cases have been reported where identified traits or processes have led to genetic improvements of GY (Sinclair et al., 2004;Yin and Struik, 2008;Passioura, 2010). The reasons lie in the many compensatory effects between traits (e.g. grain size versus grain number, light-saturated photosynthesis versus leaf surface area) and in the fact that complex traits such as GY and GPC are inherently determined at the canopy (community) level. Moreover improvements at the organ or plant level are often not translated at the canopy level (e.g. leaf level photosynthesis versus canopy radiation use efficiency; Long et al., 2006). Therefore, quantitative analyses of the 'trait hierarchy' leading to GY and GPC improvements are needed (Sinclair et al., 2004;Martre et al., 2014a).
The use of virtual crops represented by a simulation model has been proposed to help molecular biologists, physiologists, and breeders to identify relevant traits and processes and to analyse their genetic determinism (e.g. Tardieu, 2003;Hammer et al., 2006;Messina et al., 2009;Bertin et al., 2010). Crop simulation models can help develop hypotheses starting near the top of the 'trait hierarchy' leading to GY in the target environments (Sinclair et al., 2004). Before addressing the question of how to translate the information of model simulations into knowledge that can be used by physiologists or geneticists, it is appropriate to have a better understanding of the model properties and behaviour. One of the best ways to do that is to conduct a global uncertainty and sensitivity analyses of the model to investigate the behaviour of the model in response to variations in model inputs . By perturbing model parameters associated with simple physiological traits, uncertainty and sensitivity analyses allow investigation of crop responses and can help in identifying those traits that lead to a consistent high and stable GY or GPC in the target environments.
Most sensitivity analyses of crop simulation models have been performed in localized regions of the parameters' space by varying one input parameter at a time, keeping the others at their nominal values, (Asseng et al., 2002;Ruget et al., 2002;Semenov et al., 2009;Sinclair et al., 2010). This simple approach does not take into account interactions between input parameters (i.e. their context dependency) and nonlinear responses. Powerful numerical and statistical methods and tools for global model uncertainty and sensitivity analysis have been developed (Saltelli et al., 2000), but until recently they have been overlooked by crop modellers and users. In the few studies where global sensitivity methods have been used to analyse crop simulation models, only a selected number of parameters, usually <20% of the total number, thought to be the most influential ones, were studied (Makowski et al., 2006;Pathak et al., 2007;Confalonieri, 2010;Confalonieri et al., 2010a;Richter et al., 2010). Moreover, with the exception of the study of Confalonieri et al. (2010b), the uncertainty associated with climate (site) and weather (year) variability was not considered. Consideration of all the input parameters of a model and their interactions as well as variations of environments (e.g. soil characteristics and climate) and crop management [e.g. nitrogen (N) fertilizer application] is required for a robust model evaluation.
Here, the results of a global uncertainty and sensitivity analysis of SiriusQuality2, a process-based model of wheat growth, conducted at three contrasted European sites using long-term weather data and two N treatments, are reported. The results show that GY and GPC are influenced by several traits and processes, and that the ranking of the traits depended on both the environment and N supply. The high level of parameter interactions indicated that the expression of the effect of a trait at the crop level also depends on the value of the other traits.

Materials and methods
The wheat simulation model SiriusQuality2 The wheat simulation model SiriusQuality2 (http://www1.clermont.inra.fr/siriusquality/) used in this study is a revised version of SiriusQuality1 (Martre et al., 2006) and Sirius (Jamieson et al., 1998b;. SiriusQuality2 has been calibrated and evaluated for several modern wheat cultivars and tested in many environments and climates, including conditions of climate change, and crop management scenarios (Martre et al., 2006(Martre et al., , 2014bFerrise et al., 2010;Asseng et al., 2013Asseng et al., , 2014. SiriusQuality2 is a process-based model consisting of eight submodels that describe on a daily time step crop phenology (Phenology submodel), canopy development (Leaf Layer Expansion submodel), accumulation and partitioning of dry mass (DM; Light Interception and Use Efficiency and DM Allocation submodels), and N (N Allocation, and Root Growth and N Uptake submodels), including responses to shortage in the supply of soil water (Soil Drought submodel) and N (dealt with in the N Allocation submodel), and accumulation and partitioning of grain DM and N (Grain submodel). Two additional submodels describe crop evapotranspiration and soil N and water balances (Addiscott and Whitmore, 1991;Jamieson et al., 1998b). The parameters of these two submodels were not investigated in this study, since they are related to soil properties and physical constants.
The Phenology submodel calculates the duration of six development phases, comprising pre-emergence (sowing to emergence), leaf production (emergence to flag leaf ligule appearance), flag leaf ligule appearance to anthesis, anthesis to the beginning of grain fill, effective grain filling, and maturation (Jamieson et al., 1998a). The anthesis date is mainly determined by the rate of leaf production (1/ phyllochron) and the final leaf number. The final leaf number is calculated by daylength (photoperiod) and temperature (vernalization) response subroutines (Brooking et al., 1995;Jamieson et al., 1998a;He et al., 2012). Canopy development is simulated in the Leaf Layer Expansion submodel as a series of leaf layers associated with individual stems (Fournier et al., 2005;Lawless et al., 2005).
The Light Interception and Use Efficiency submodel calculates the amount of light intercepted by each leaf layer using the turbid medium approach (Monsi and Saeki, 2005) and uses it to produce biomass at an efficiency [radiation use efficiency (RUE)] calculated from temperature, air CO 2 concentration, soil water deficit, leaf N mass per unit surface area (specific leaf N, SLN), and the ratio of diffuse to direct radiation . The total aboveground biomass at any time is the sum of the daily rate of biomass accumulation of each leaf layer, which, in turn, is the product of RUE and intercepted photosynthetically active radiation by the leaf layers (Monteith, 1977). The canopy light extinction coefficient (K L ) is assumed to be independent of N and water shortages (Robertson and Giunta, 1994). Part of the biomass produced each day is allocated to the stem so that a target constant specific DM of the leaf laminae (SLW p ) and sheath (SSW p ) is maintained. After the end of the endosperm cell division period, all new biomass plus a constant proportion of the labile leaf and stem biomass at the end of the endosperm cell division phase are allocated to grains.
Crop N uptake is driven by canopy expansion. The vertical distribution of leaf N follows the light distribution (Bertheloot et al., 2008b), and the ratio of N to light extinction coefficients is determined by the crop N status (Moreau et al., 2012). As for the biomass, any N not allocated to the leaves is allocated to the stem if its N concentration is less than its maximum. After anthesis, the capacity of the root system to take up N from the soil decreases linearly with thermal time (Martre et al., 2006). After the end of the endosperm cell division stage, the rate of N transfer to grains is determined by the stem and leaf N concentrations, and follows first-order kinetics (Bertheloot et al., 2008a, b).
Seventy-five parameters are used to parameterize the processes described above (Table 1), which were subjected to the uncertainty and sensitivity analyses described below.
Simulation set-up Two French sites, Avignon (AV) and Clermont-Ferrand (CF), and one UK site, Rothamsted (RR), were selected for this study ( Table 2). The characteristics of the soil used at the three sites are summarized in Table 3. The soil water-holding characteristics at CF and RR were similar to that previously reported for these two sites (Jamieson et al., 1998b;Martre et al., 2006), while at AV parameters corresponding to a sandy soil were used to ensure water-limited conditions for most years. Simulations were conducted with historical weather data of 40 continuous years  at each site. Recommended sowing dates were used at each site ( Table 2). The sowing density was 250 seeds m -2 at the three sites. Soil N characteristics (organic and inorganic N content) were the same at all sites (Table 3).
Two N treatments were considered at each site, a high N treatment [HN, 20 g (N) m -2 ] corresponding to the commercial N rate in the investigated environments with a late N application to optimize GPC, and a low N treatment [LN, 4 g (N) m -2 ] corresponding to 20% of HN. For HN, N fertilizer was applied as four split dressings at the growth stages GS21 [4 g (N) m -2 ], GS31 [8 g (N) m -2 ], GS37 [6 g (N) m -2 ], and GS65 [2 g (N) m -2 ] described by the Zadoks scale (Zadoks et al., 1974). For LN, N fertilizer was applied once at GS30. An irrigation of 10 mm was applied at sowing to ensure good crop emergence conditions, and 5 mm irrigations were applied on the dates of N applications to ensure dissolution of N fertilizer in the soil solution. This limited amount of water did not have any other effect on crop development or growth.

Morris screening analysis
The 75 plant parameters of SiriusQuality2 were first analysed using the one-at-a-time (OAT) Morris method (Morris, 1991) to rank them in order of importance and to identify parameters involved in interactions with other parameters or parameters whose effect are non-linear.
The Morris method is based on individually randomized OAT experiments, in which the impact of varying each of the chosen parameters is evaluated in turn. In this method, each vector of parameters is first discretized into p levels corresponding to the quantiles of the parameter distribution and defining a k-dimensional and p-level region of experiment in the parameter space, where k is the number of parameters (Saltelli et al., 2004). The sampling strategy proposed by Campolongo et al. (2007) was used, where r elementary effects are sampled to construct r trajectories of (k+1) points in the region of the experiment, each providing k elementary effects, one per input parameter. The number of model executions is thus equal to r×(k+1). Here, p was set to 8, and 760 vectors of parameters corresponding to 10 trajectories were generated using the uncertainty and sensitivity analysis software SimLab version 2.2 (Saltelli et al., 2004). The parameters were assumed to be uniformly distributed in [0, 1] and then rescaled from the unit hypercube to their actual values using the following equation: where X i m and X i d are the modified and default values of the parameter x i , respectively; ∆ i min and ∆ i max are the minimum and maximum fractional variations of the parameter x i , respectively; and P is a random number following a uniform distribution in [0, 1]. The space of each parameter was set to be ±20% of its default value, which encompass the observed range of genetic variation for most of the parameters (Martre et al., 2007).
After running SiriusQuality2 for each of the parameter sets under each of the site-N treatment-year combinations, the model outputs were post-processed using SimLab. For each parameter, two sensitivity indices were obtained from the distribution of their elementary effects. The first is the mean of the distribution of the absolute values of the elementary effects (μ*), which estimates the overall importance of the parameter on the model outputs (Campolongo et al., 2007). The second is the standard deviation of the distribution of elementary effects (σ), which indicates either a parameter interacting with other parameters or a parameter whose effect is non-linear.
In order to plot heat maps of μ* and σ, they were rescaled in [0, 1] since different output variables may have different magnitudes of μ* and σ. For each site-N treatment-output variable combination, the median of μ* and σ over the 40 years of simulation was first calculated and then rescaled using the formula: where z i, j r , z i, j min , and z i, j max are the rescaled, minimum and maximum values of μ* or σ for the ith parameter and jth output variable, respectively. The sign of μ* was retrieved by conducting an independent OAT experiment. Each parameter was changed in turn by ±50% of its default value in 10% increments. Then for each parameter-site-N treatment combination, the slopes of the linear regression between the model outputs and the parameters were calculated and the sign of the slopes was given to the rescaled μ*. This provided useful information of the sign of the elementary effects, but it should be treated with caution since it might be significantly overestimated in the case of non-monotonic elementary effect distribution. For instance, a parameter with a large and negative μ* value should be considered as very significant, but the sign of the effect might be negative or positive depending on the position of the parameter input space where it is sampled, and its overall effect is therefore not necessarily negative. However, in a complex model such as Table 1. Name, symbol, definition, nominal value, and unit of the 75 parameters of the wheat simulation model SiriusQuality2 All of the parameters belong to eight submodels of Phenology, Leaf layer expansion, Light interception and use efficiency, Grain, DM allocation, N allocation, Root growth and N uptake, and Soil drought factors.    SiriusQuality2, non-monotonic elementary effect distributions are infrequent (Saltelli et al., 2008).

Extended Fourier amplitude sensitivity test
The influential parameters identified by the Morris method were further analysed using the quantitative extended Fourier amplitude sensitivity test (E-FAST; Saltelli et al., 1999). As for the Morris analysis, the space of each parameter was set to be ±20% of its default value assuming a uniform distribution. The k input parameters were coded so that their domain of variation was [0, 1] k . A total of 29 889 parameter vectors were generated using SimLab (i.e. 729 different positions in the input space were considered for each parameter), which is much higher than the minimum number of 2665 required by the E-FAST method with the number of parameters considered in this study. The generated random values of the parameters were then rescaled to their actual values using Equation 1, and SiriusQuality2 was run for all the site-N treatment-year combinations, resulting in 7.17 × 10 6 model runs from which the Fourrier spectrum was calculated in SimLab on the different model outputs to estimate the first-order (S i ) and total (S Ti ) effect of each of the 41 influential parameters for each site-N treatment combination.

Grain yield and protein concentration uncertainty and environment by N interactions
The three sites used for this study differed in terms of their climate (Table 2; Supplementary Fig. S1 available at JXB online) and soil (Table 3). AV has a Mediterranean climate with warm and dry summers and mild wet winters. CF has a Continental climate with hot summers and cold and relatively dry winters, while RR has an Oceanic climate with moderately cool summers and comparatively warm winters and evenly distributed precipitations. These sites are representative of the climate diversity of the wheat-growing regions in Europe. The site differences, as felt by the crops, are illustrated by the dynamics of the fraction of transpirable soil water (FTSW) during the growing season (Fig. 1). AV was characterized by high probabilities of mild water deficits early in the season and during the stem extension period, and of severe water deficits during grain filling. At this site, for 25% of the years, simulated FTSW decreased below 0.65, the threshold value for reduction of leaf expansion (LowerFSTWexpansion), after the beginning of stem extension, and for 55% of the years leaf senescence and stomatal conductance were significantly accelerated by water deficit (i.e. FTSW <0.5) during the grain-filling period.
Comparably, almost no water limitation occurred during the vegetative growth period at RR. At this site, mild water deficits occurred during the second half of grain filling for 55% of the years, and for most of the years it had no significant effect on GY and GPC (data not shown). The situation at CF was intermediate to that at the other two sites. This site was characterized by late water deficit. Although at a given site the patterns of FTSW among the simulated years were similar (data not shown), the levels of water deficits were very variable (as illustrated by the box plots in Fig. 1).
For example, at AV, the 10% and 90% percentiles of FTSW at anthesis were 0.31 and 0.87, respectively.
Simulated GY was significantly lower and more variable at AV than at the two other sites ( Fig. 2A), especially under HN, where it ranged from 0.14 kg DM m -2 to 0.77 kg DM m -2 . Under HN, for most of the years, GY was similar at RR and CF, but under LN on average it was 6% higher at CF than at RR. At the three sites, GY was more variable under HN than under LN. On average GY was 88% higher under HN than under LN. As for GY, GPC was more variable at AV (ranging from 6.2% to 12.2% under LN and from 10.5% to 16.5% under HN) than at the two other sites. As for GY, the model predicted significant site×N interactions for GPC (Fig. 2B). Although it is usually observed that GPC is higher under elevated temperature or water deficit because carbon assimilation is more affected by these environmental constraints than N assimilation (Triboi and Triboi-Blondel, 2002;Aguirrezábal et al., 2014), simulated GPC under HN was lower at AV than at the two other sites despite AV being a more stressful environment. This result can be explained by significant water deficit during the stem extension period at AV, which reduced the expansion of the canopy, and therefore limited N storage in vegetative tissues during the pre-anthesis period. Severe water deficits at this site also reduced soil N supply to the plant due to their negative effects on organic N mineralization. Under LN, GPC was higher at AV than at the two other sites, because of the strong limitation of crop growth and N accumulation, independently of the climate. For both N treatments, CF had higher GPC than RR. This was expected since CF had a higher probability of water deficit and supraoptimal temperature during grain filling.

Screening of the model input parameters
The influence of all 75 parameters of the SiriusQuality2 model (Table 1) was first assessed using the Morris screening method (Morris, 1991). The perturbation of the parameters resulted in -14 d to +17 d differences in anthesis date (Supplementary Table S1 at JXB online). GY (-45% to +16%) and GPC (-20% to +35%) also showed large variations in response to parameter variations. The responses of eight key model outputs were analysed in detail. Grain protein deviation (GPD), namely the residual from the regression of GPC on GY, was also calculated (Monaghan et al., 2001). Several authors have used GPD to identify genotypes with higher GPC than would be expected from their GY (e.g. Monaghan et al., 2001;Oury and Godin, 2007;Bogard et al., 2010). For each site-N-year combination, GPD was calculated by regressing GPC against GY for the 760 parameter combinations ('genotypes') generated in the Morris analysis. Under HN at AV the regression was significant in only 75% of the years and the median value of r 2 was only 0.05 (Supplementary Table S2). For the other site-N combinations, the correlations were highly significant in >95% of the years, and the ranges and distributions of the slope [median value -15.6% kg -1 (DM) m 2 ] and r 2 (median value 0.53) were in close agreement with that calculated from experimental data (e.g. Oury et al., 2003;Bogard et al., 2010).
The 'overall' effect (μ*) and the non-linear and/or interaction effects (σ) of the parameters were well correlated for all model outputs, as illustrated in Fig. 3 for GY, GPC, and GPD at CF under LN. The influence of the parameters on the model outputs showed a high level of uncertainty associated with the year, and at a given site the ranking of the parameters depended strongly on the years.
The median values of the rescaled (see the Materials and methods) 'overall' effect ( Fig. 4) and of the non-linear and/ or additive effects ( Supplementary Fig. S2 at JXB online) are summarized on heat maps, which give a grand-view of parameter×site×N interactions. Anthesis date was most influenced, in order of importance, by the phyllochron (P) and the absolute minimum possible leaf number (MinL), the duration of the period between the appearance of the flag leaf ligule and anthesis (PFLLAnth), and the response of vernalization rate to temperature (VBEE; Fig. 4). Except for the maximum temperature for vernalization (MaxTvern), which had no significant overall effect on anthesis date, all the phenology parameters had significant overall and non-linear and/or interaction effects on canopy expansion and biomass and N accumulation, grain N yield, post-anthesis N uptake, and GPD ( Fig. 4; Supplementary Fig. S2). As expected, the influence of the phenology parameters on anthesis date was independent of N supply; however, their effect on the other model outputs were highly dependent on N supply. In all site-N combinations, all the parameters associated with phenology had an effect of the same sign on anthesis date and crop DM at maturity, but of opposite sign on anthesis date and GY. In other words, early flowering crops had higher GY associated with higher DM accumulation after anthesis and DM harvest index. Three parameters associated with leaf layer expansion (AreaPL, NLL) and light interception (Kl) had some indirect effects on anthesis date due to non-linearity or interactions between parameters ( Supplementary Fig.  S2). This result could be explained by their effect on canopy temperature.
At CF and RR under HN, the green area index (GAI) at anthesis was essentially influenced in order of importance by the number (NLL) and the potential size (AreaPL and aSheath) of the leaves produced after floral initiation, and the potential ratio of the flag leaf to penultimate leaf size (RatioFLPL). Under these growing conditions, the potential phyllochronic duration between the end of expansion and the beginning of senescence of the leaves produced after floral initiation had some non-linear and/or interaction effects on GAI ( Supplementary Fig. S2 at JXB online). Under dry conditions and non-limiting N supply NLL, P, AreaPL, MinL, and the FTSW threshold for the response of leaf expansion to water deficit (LowerFSTWexp) were the most influential parameters (in order of importance) for GAI at anthesis ( Fig. 4; Supplementary Fig. S2). For most of the model outputs and parameters, the sensitivity indices were significantly higher under HN than under LN. The only significant exceptions were GAI at anthesis and total crop DM at maturity which were more influenced by the parameter related to N allocation under LN than under HN. Under LN, canopy expansion was primarily determined by N supply; therefore, it was most influenced by the N requirement of growing leaf tissues (CritSLN). It was also significantly influenced by RUE and Kl. At all site-N combinations, RUE had the highest influence on GY. Under HN, at CF and RR the potential duration of grain filling (Dgf) was the second most influential parameter. In contrast, at AV, Dgf had no significant overall effect on GY. At the three sites under LN, the phyllochron was the second most influential parameter with respect to GY. At CF and RR, the potential number of culm leaves (NLL) had a significant overall effect on GY under both N treatments, while the potential size of the culm leaves (AreaPL) and Kl had significant overall effects on GY only under HN. NLL and GY were negatively related under all site-N combinations. At RR under LN, the relative rate of accumulation of grain structural DM during the endosperm cell division period (Kcd) and the duration of the endosperm cell division period (Dcd) were significantly associated with GY uncertainty. The parameters determining the response of RUE to the ratio of diffuse to direct light (SlopeFR), leaf N (TauSLN), temperature (Tmax and Topt), and water deficit (LowerFTSWrue and UpperFTSWrue) had much less effect on GY than the parameters mentioned above. At AV, the response of leaf senescence to water deficit (LowerFTSWsen and MaxDSF) and the capacity of the crop to extract soil water (MaxRWU) had a high overall effect on GY, especially under HN. None of the parameters associated with DM or N partitioning had a significant effect on GY.
As regards GY, at all site-N combinations, RUE was the parameter influencing GPC the most. The effect of the other parameters on GPC depended on both the site and the N supply (Fig. 3). The phyllochron was the second most important parameter for GPC at AV (under both N), while at the two other sites it was the potential duration of grain filling (Dgf) under HN and the number of culm leaves (NLL) under LN. At all site-N combinations, Kl and Dgf were among the five most influential parameters. Under HN at the three sites, the N storage capacity of the leaves (MaxSLN) and NLL had the strongest overall positive effect on GPD. This was obtained despite a negative effect of these two parameters on post-flowering N uptake. Under LN, RUE had a strong negative effect on GPD. Kl also had significant overall negative     . The median of μ* was rescaled to [0, 1] across the sites and N treatments so for a given output they can be compared across the sites and N treatments. Negative values (blue colour) indicate that the parameter negatively influences the corresponding model output, and vice versa. The parameters were grouped according to the submodel to which they belong.
effects on GPD at the three sites. All the parameters with a significant influence on GY or GPC had opposite effects on these two variables. However, the phyllochron, at the three sites and both N treatments, the parameters related to the response of leaf senescence to water deficit (UpperFTSWsen and MaxDSF), and the capacity of the crop to extract soil water (MaxRWU), at AV and to a lesser extent at CF, had an effect of the same sign on both GY and GPD.

Quantitative sensitivity analysis of the influential parameters
The Morris design is an efficient screening method to identify influential model parameters. However, it does not allow quantification of the interactions between the model parameters, and the Morris sensitivity indexes are only qualitative. Therefore, based on the Morris analysis, 41 influential parameters were identified (Table 4) and their individual first-order (S i ) and total (S Ti , first-order plus interactions of any order with all the other parameters) effects were calculated using the E-FAST method (Saltelli et al., 1999). Although the number of parameter combinations generated in the E-FAST analysis was considerably higher than that generated in the Morris analysis, the range of variations in anthesis date, GY, and GPC was similar for the two methods (Supplementary Table S3 at JXB online). In general, parameter ranking for the two methods was in good agreement ( Supplementary Fig. S3). The major difference was that in the Morris method the phenology parameters had greater effects on DM and N accumulation, GPC, and GPD under HN than under LN, while the opposite was observed with the E-FAST method (Figs 4, 5).
The analysis of the results of the E-FAST method was focused on the most influential parameters. The parameters accounting for 90% of the sum of S Ti for GY, GPD, or GPC in at least 50% of the years and for at least one site and N treatment were selected. For each site-N treatment-model output combination, 6-17 parameters were selected, representing a total of 32 different parameters. Averaged across years these 32 parameters accounted for 48-87% of S i for GY, GPC, and GPD, while the interactions between all of the parameters accounted for 8-56% of S Ti . For GY and GPC, the contribution of the interactions was 2-fold higher under LN (averaging 34%) than under HN (averaging 16%; Supplementary  Fig. S4 at JXB online). At the different sites and for both N treatments, RUE and Dgf accounted for 22-71% of the sum of S i with respect to GY, and GPC, except at AV where Dgf had no significant first-order effect on GY ( Supplementary  Fig. S4). At AV, maxRWU contributed 12% and 8% of the sum of S i for GY under HN and LN, respectively. These results suggest that there are a few very influential parameters and a large majority of non-influential parameters. However, since the interactions between the parameters accounted for a large part of the total variance, the relative contribution of the different parameters to the variance of the model outputs cannot be quantified only from their first-order effect.
The S Ti results showed that under HN the phenology parameters had much less influence on GY than estimated through the Morris screening analysis (Fig. 6A). At CF, GY was mainly influenced by the duration of the period between the appearance of the flag leaf ligule and anthesis (PFLLAnth), while at RR it was mainly influenced by the phyllochron. The ranking of the other parameters was similar to that observed in the screening analysis, and this quantitative analysis confirmed the primary role of RUE in determining GY under HN. As shown by the screening analysis, under LN, GY was determined by a large number of parameters (14-17 depending on the site) with equivalent effects (Fig. 6A). As for GY, the ranking of the parameters with respect to GPC and GPD was similar for the two methods (Fig. 6B, C). At AV under HN, the S Ti value of the three most influential parameters (RUE, P, and MaxStemN) with respect to GPC and GPD showed a high dependency on the year, as indicated by the spread of the error bars.

Discussion
Crop simulation models based on physiologically sound mechanisms have the potential to accommodate the uncertainty associated with gene and environment context dependencies (Campos et al., 2004;Sinclair et al., 2004; Hammer, 2005;Hammer et al., 2006). The analysis of system dynamics and properties through simulation experiments using such models provides a holistic approach to study soilplant-atmosphere systems through the properties captured in the model. The in-depth simulation analysis carried out here contributes to the identification of traits with stable effects in a given population of environments needed to implement fully physiological trait-based breeding, for which phenotyping solutions are being developed (Lopes et al., 2014;Pask et al., 2014).
The results of the E-FAST analysis are quantitative for the 41 selected parameters analysed. However, the results might be taken as qualitative when considering the 75 parameters of the model, since possible interactions between the selected parameters and the remaining parameters that were unchanged in the analysis are not accounted for. However, since the 34 parameters that were unchanged in the E-FAST analysis were selected based on a screening Morris analysis, it is possible to be confident that the E-FAST results can be taken as quantitative when extended to the original (75-parameter) model.
In good agreement with the pleiotropic nature of the genetic and physiological antagonisms between GY and GPC (Aguirrezábal et al., 2014), all the parameters had antagonist effects on these two traits. The present results show that GPC can be increased by increasing the N storage capacity of the leaves and to a lesser extent that of the stem. In most cases, and in good agreement with reported phenotypic correlations (Bogard et al., 2010), a higher GPD was associated with a higher post-anthesis N uptake.
Overall, RUE and Dgf were the two parameters that influenced GY the most. Large genetic variability for RUE has been reported for wheat (García et al., 2014), and genetic improvement in GY in the UK after 1983 has been positively associated with higher RUE (Shearman et al., 2005). Genetic analyses have found genetic variations for the duration of grain filling of up to 40% for both bread (Charmet et al., 2005) and durum (Akkaya et al., 2006) wheat. Therefore, the present study supports earlier suggestions (Evans and Fischer, 1999;Reynolds et al., 2012) that RUE and Dgf could contribute to improve GY in wheat. However, according to the results, Dgf had no significant effect on GY in Mediterranean environments with terminal water deficit. To a lesser extent, the effect of RUE on GY was also larger and more stable (less interannual variability) in wet than in dry environment.
Changing the number of culm leaves (NLL) independently of the final leaf number had a major influence on GY. With a constant phyllochron, increasing NLL could be an approach to increase the duration of the stem extension period, which has been suggested to increase ear dry mass at anthesis and grain number per square metre and thus yield potential (Slafer et al., 2001;Foulkes et al., 2011). However, in the present study, NLL was negatively associated with GY (i.e. increasing NLL decreased GY), and the influence of NLL on GY was higher under HN than under LN. This result is explained by the fact that even with 20 g N m -2 of N fertilizer in HN, the higher early growth of simulated genotypes with higher NLL value depleted soil N, and the growth of the upper culm leaves was then limited by soil N supply. A higher number of culm leaves thus resulted in a reduction in the amount of light capture during the stem extension period which had a negative effect on ear growth and thus on grain number. The production of biomass during the grain-filling period was also reduced with higher NLL values.
Therefore, although increasing the number of culm leaves might be a good strategy to improve GY, it would require increasing the N supply and probably the number of fertilizer applications. A more precise monitoring of crop N needs would also be required so that leaf growth is not limited by N and applied N is not lost by leaching. This analysis of the effect of NLL is thus a good example where the complex response of the crop can hardly be understood through experimentation because of the complex interactions that occur. It is also an example where agronomic solutions should be developed along with genetic solutions in order to take full advantage of the genetic improvements.
Physiological constraints should also be considered when trying to manipulate plant development. Indeed changing the number of culm leaves independently of the final leaf number and therefore without affecting the anthesis date might be challenging. The leaf stage at the occurrence of the double ridge and terminal spikelet stage is tightly related to the final leaf number in wheat (Slafer et al., 1994;Giunta et al., 2001;Jamieson et al., 2007). However, variations of about one leaf in this relationship have been observed between winter and spring near isogenic lines with the same final leaf number (Brooking and Jamieson, 2002). This result suggests that the range of variation tested in the present study (±0.9 leaf) could be achieved by genetics without changing the final leaf number and therefore the anthesis date.
Parameters related to phenology had a significant influence on almost all of the investigated model outputs and, as expected, their influence on GY was higher in dryer environments. In general, early flowering genotypes had higher GY because of a higher post-anthesis biomass production and higher harvest index. This trait has played a major role in the adaption of wheat to the Mediterranean climate (Richards, 1991). A recent analysis of genetic improvement in the CIMMYT spring wheat breeding programme has also shown a significant association of GY improvement with fewer days to anthesis (Lopes et al., 2012). The importance of flowering time has been highlighted in many genetic analyses where co-locations of quantitative trait loci (QTLs) for flowering time and GY or GPC were found (e.g. Habash et al., 2007;Laperche et al., 2007). In particular, Bogard et al. (2011) showed that the phenotypic correlation between monocarpic leaf senescence and GY and GPD is largely explained by anthesis dates, as major QTLs for these traits co-locate with major flowering time QTLs. Consequently, the benefits of using populations with a narrow range of phenology to identify chromosomic regions associated with adaptive traits have been demonstrated (e.g. Olivares-Villegas et al., 2007).
The present results show that phyllochron could be a valuable trait to improve GY. Among the phenology parameters, phyllochron appeared as the most influential. Phyllochron is a key determinant of early vigour, especially in growing conditions where yield is mainly determined by pre-anthesis biomass production (Rebetzke and Richards, 1999), which can in turn affect tillering through the carbohydrate supply/demand balance (Bos and Neuteboom, 1998;Rebolledo et al., 2012). In the analysis reported here, under LN, phyllochron was the most important parameter determining GY at AV, and the second most important at the other two sites. The importance of phyllochron under LN was mainly due to its strong interactions with other parameters. Interestingly, phyllochron was the only parameter that had an effect of the same sign on both GY and GPD at the three sites. Significant genetic variations exist for phyllochron in wheat (He et al., 2012), and specific QTLs were found in durum wheat (Sanna et al., 2014), rice (Miyamoto et al., 2004;Morita et al., 2005), and barley (Borràs-Gelonch et al., 2010).
If physiological traits are to be used in crop breeding programmes, there is a need to develop protocols for their efficient phenotyping. For phyllochron, Jamieson and Munro (1999) have proposed a simplified field-based protocol that could be implemented in a breeding programme. In contrast, determining the number of extended internodes and therefore the number of culm leaves might be more difficult to quantify with the necessary precision (Weightman et al., 1997), and it might require determination of the size of the leaves associated with each phytomer under near optimal growth conditions (Fournier et al., 2005), which in overwintering conditions would require at least two plant samplings. Regarding RUE, high-throughput protocols in the greenhouse are currently being developed and will provide a relevant alternative to field protocols that use punctual measurements of canopy light interception (Moreau et al., 2012;García et al., 2014). Finally, for Dgf, visual assessment of physiological maturity might be a valuable alternative to ear sampling (Hanft and Wych, 1982;Talbert et al., 2001).
In the Morris analysis, μ* and σ were closely correlated. Therefore, a parameter that is important in the model is usually also involved in non-linear or interaction effects. In other words, the 'overall' importance of a model parameter is primarily determined by the non-linear response of the model or by its interactions with other model parameters. The interactions were quantified in the E-FAST analysis: 8-56% of the variance in GY, GPC, and GPD were accounted for by the interactions between the parameters. Therefore, the effect of a parameter cannot be accurately quantified based on its firstorder effect. The effect of the parameters on GY was more balanced under LN than under HN. The higher trait×trait interactions found in the present study under LN compared with HN may explain, at least in part, the lower genetic progress in GY under LN than under HN (Le Gouis, 2011). This result has important implications for physiological breeding as the introgression of traits whose value depends on interactions with other traits in an improved genotype will be more difficult than in the case of traits with lower interactions. This also calls for the definition of ideotypes that combine several traits (Semenov and Stratonovitch, 2013;Martre et al., 2014a).
The interactions between parameters found in this study depend on the hypothesis on which the model is founded regarding the interactions between physiological processes and forcing (environmental) factors for particular processes. Recent studies have shown that multimodel ensembles give more accurate predictions of GY (Asseng et al., 2013;Bassu et al., 2014) and of processes leading to yield than any individual model; that is, ensemble estimators have less compensation of errors than individual models and therefore causal relationships are more likely (Martre et al., 2014b). The difficulty in the use of multimodel ensembles to identify putative traits that could be used for genetic improvement would be to achieve the correspondence between the parameters of individual models.
In summary, an in-depth simulation exercise was conducted using the wheat simulation model SiriusQuality2 to quantify the value of putative morpho-physiological traits for breeding and their environmental stability. To the best of the authors' knowledge, it is the first study analysing the uncertainty of the whole set of parameters of a process-based crop simulation model. The results show that the overall effect of the parameters was dominated by their interactions with other parameters. Under non-limiting N supply, a few influential parameters with respect to GY were identified (e.g. RUE, Dgf, or P). However, under LN, >10 parameters showed similar effects on GY and GPC. All parameters had opposite effects on GY and GPC, but leaf and stem N storage capacity appeared as good candidate traits to shift the negative relationship between GY and GPC, confirming the results of an earlier analysis (Martre et al., 2007).

Supplementary data
Supplementary data are available at JXB online. Figure S1. Long-term weather at the three sites considered in this study. Figure S2. Heat map of the median values of the rescaled standard deviation of the distribution of elementary effect from the Morris screening analysis. Figure S3. Comparison of parameter rankings obtained for the Morris analysis and for the E-FAST analyses. Figure S4. Radar plots of the median values of the E-FAST first-order sensitivity index for the 32 most influential parameters with respect to grain yield, grain protein concentration, and grain protein deviation. Table S1. Ranges of differences in anthesis date and percentage changes in grain yield and grain protein concentration explored in the Morris screening analysis. Table S2. Statistics of the linear regression between simulated grain protein concentration and grain yield for the 760 parameters combinations of the Morris screening analysis. Table S3. Ranges of differences in anthesis date and percentage changes in grain yield and grain protein concentration explored in the E-FAST analysis.