# Genetic control of plasticity of oil yield for combined abiotic stresses using a joint approach of crop modeling and genome-wide association

Brigitte Mangin, Pierre Casadebaig, Eléna Cadic, Nicolas Blanchet, Marie-Claude Boniface, Sebastien Carrere, Jerome Gouzy, Ludovic Legrand, Baptiste Mayjonade, Nicolas Pouilly, et al.

<figure><img src="image_1.png" /></figure>

## To cite this version:

Brigitte Mangin, Pierre Casadebaig, Eléna Cadic, Nicolas Blanchet, Marie-Claude Boniface, et al. Genetic control of plasticity of oil yield for combined abiotic stresses using a joint approach of crop modeling and genome-wide association. Plant, Cell and Environment, 2017, 40 (10), pp.2276-2291. 10.1111/pce.12961. hal-01608380

HAL Id: hal-01608380

https://hal.science/hal-01608380

Submitted on 26 May 2020

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Genetic control of plasticity of oil yield for combined abiotic stresses using a joint approach of crop modeling and genome-wide association

Brigitte Mangin¹a, Pierre Casadebaig²a, Eléna Cadic<fnref n="1" />, Nicolas Blanchet<fnref n="1" />, Marie-Claude Boniface<fnref n="1" />, Sébastien Carrère<fnref n="1" />, Jérôme Gouzy<fnref n="1" />, Ludovic Legrand<fnref n="1" />, Baptiste Mayjonade<fnref n="1" />, Nicolas Pouilly<fnref n="1" />, Thierry André<fnref n="3" />, Marie Coque<fnref n="4" />, Joël Piquemal<fnref n="5" />, Marion Laporte<fnref n="6" />, Patrick Vincourt<fnref n="1" />, Stéphane Muños<fnref n="1" />, Nicolas B. Langlade<fnref n="1" />

## Corresponding author

Nicolas Langlade (email: nicoas.lnaglade@inra.fr)

## Addresses

This article has been accepted for publication and undergone full peer review but has not been through the copyediting, typesetting, pagination and proofreading process which may lead to differences between this version and the Version of Record. Please cite this article as doi: 10.1111/pce.12961

<footnote>aThese authors contributed equally to this work</footnote>

<footnote n="1">¹ LIPM, Université de Toulouse, INRA, CNRS, Castanet-Tolosan, France</footnote>

<footnote n="2">² AGIR, Université de Toulouse, INRA, INPT, INP-EI PURPAN, Castanet-Tolosan, France</footnote>

<footnote n="3">³ SOLTIS, Domaine de Sandreau, Mondonville F-31700 Blagnac, France</footnote>

<footnote n="4">⁴ BIOGEMMA, Domaine de Sandreau, Mondonville F-31700 Blagnac, France</footnote>

<footnote n="5">⁵ SYNGENTA SEEDS, 12 chemin Hobit, F-31790 Saint Sauveur, France</footnote>

<footnote n="6">⁶ RAGT 2n, Site de Bourran, F-12000 Rodez, France</footnote>
# Abstract

Understanding the genetic basis of phenotypic plasticity is crucial for predicting and managing climate change effects on wild plants and crops. Here, we combined crop modeling and quantitative genetics to study the genetic control of oil yield plasticity for multiple abiotic stresses in sunflower.

First we developed stress indicators to characterize 14 environments for three abiotic stresses (cold, drought and nitrogen) using the SUNFLO crop model and phenotypic variations of three commercial varieties. The computed plant stress indicators better explain yield variation than descriptors at the climatic or crop levels. In those environments, we observed oil yield of 317 sunflower hybrids and regressed it with three selected stress indicators. The slopes of cold stress norm reaction were used as plasticity phenotypes in the following genome-wide association study.

Among the 65,534 tested SNP, we identified nine QTL controlling oil yield plasticity to cold stress. Associated SNP are localized in genes previously shown to be involved in cold stress responses: oligopeptide transporters, LTP, cystatin, alternative oxidase, or root development. This novel approach opens new perspectives to identify genomic regions involved in genotype-by-environment interaction of a complex traits to multiple stresses in realistic natural or agronomical conditions.

## Brief Summary

First, we developed a novel method to describe the abiotic stresses perceived by the plant using the dynamic ecophysiological model included in the SUNFLO crop model. This allowed us characterize 17 environments corresponding to a large range of cold, drought and nitrogen stresses situations.

Secondly, we characterized, in the 17 environments, the oil yield plasticity to those three abiotic stresses in a sunflower core-collection, defined as the linear response of oil yield to the modeled stresses. This allowed us to perform a genome-wide association study on these oil yield plasticity phenotypes. In this article, we describe the GWAS only for cold stress, for which we identified 9 regions bearing candidate genes. We restricted the publication of our results to cold because drought and nitrogen stress plasticities were genetically correlated and we believe the GWAS results on those phenotypes must be interpreted with caution. Nevertheless, the analysis on cold stress identified very promising candidate genes and represents a proof of concept of the general approach. All together, our study represents a novel combination of crop modeling and quantitative genetics that addresses an important question for biology.
# Introduction

Adaptation to climate change requires new crop varieties adapted to new management options. Adaptation of agriculture is a key factor to lessen the impact of climate change (Lobell et al., 2008). The crop exposition to unfavorable growing periods can be partially controlled by adapting crop management, i.e by shifting planting dates or choosing a cultivar with an adequate phenology (Acosta-Gallegos and White, 1995). But such adaptations also have side-effects and while flowering can be successfully desynchronized from the period of occurrence of water deficit, crop emergence would be more exposed to cold stress. In this case, adaptation to climate change should also include the development of new crop varieties (Rosenzweig et al., 1994), with new or improved properties such as tolerance to cold or other abiotic stresses. Moreover, because the multiplicity of cultivation conditions (soils, climatic uncertainty), a single genotype can be exposed to random unfavorable growing conditions within its cultivation area, ultimately impacting the expected crop performance. To ensure a stable performance under uncertain conditions, newly developed genotypes not only need to be tolerant, i.e. adapted to a single type of environment (specialisation) but also plastic, i.e. to be able to adapt in most of growing conditions encountered in the targeted cultivation area (Sambatti and Caylor, 2007).

Phenotypic plasticity is a key process for crop productivity under climate change. One way plants will respond to these changes is through environmentally induced shifts in phenotype (phenotypic plasticity) (Nicotra et al., 2010). While the process of phenotypic plasticity is mostly studied on natural systems, its implications on crop productivity under climate change are interesting for plant breeding (DeWitt and Langerhans, 2004; Sadras et al., 2009). Empirically relationships between plant traits and environmental variables (norms of reaction) are known to vary within species in nature as well as in crop species, are considered heritable traits themselves subject to natural or artificial selection (Via and Lande, 1985; Sambatti and Caylor, 2007). However, crop growth in a fluctuating environment generates complex and dynamic interactions between plant and environment, under the control of cultural practices.

It is necessary to unravel and measure abiotic stress levels before assessing plasticity in plant traits. In these conditions, assessing plasticity in plant traits is limited by our capacity to unravel those interactions and estimate abiotic stress levels at the plant scale. Actually, each growth condition creates a unique combination of those stress levels, with possible identical combinations in different growth conditions. For example, crops growing in continental climates might be exposed to both cold (during emergence) and heat (during flowering) stresses; with high temperatures driving a strong evaporative demand and water deficit, which also limits plant nitrogen uptake from the transpired water stream (Kiani et al., 2016). Accordingly, in a given cultivation area, stresses are not independent (Vile et al., 2012) and need to be characterized and modelled prior studying their impact on plant traits.

Crop simulation and modeling can help to characterize environment from the crop point of view. Because the environment is the largest component of the phenotypic variability of most plant traits and of course crop yield, its quantitative characterization is of major importance (Lake et al., 2016). Crop simulation models are based on mathematical equations representing the crop growth and development as a function of environment (climate, soil and management). Such tools can give access to plant-level state variables, such as time-series of several abiotic stresses, in large range of growing conditions otherwise difficult to characterize with sensors. This methodology
was recently implemented and allowed to identify major types of water deficit patterns for rainfed wheat in the Australian target population of environments (Chenu et al., 2013) or for coupled thermal and water stress patterns for chickpea in Australian National Variety Trials (Lake et al., 2016). In sunflower, a crop model was developed by Casadebaig et al. (2011) and takes into account water, cold, heat, and nitrogen stresses to estimate their impact on grain yield and oil content.

Genetic studies of genotype-by-environment interaction (GxE) Although plasticity has long been recognized as an interesting trait in ecological and crop science since the pioneering work of Bradshaw in the 1950s, the identification of genetic variation involved in plasticity is scarce and recent. An approach to deal with phenotypic variation due to environmental effects is to develop multi-environmental QTL analysis to identify QTLxE interactions as reviewed in Van Eeuwijk et al. (2010). Although this approach was widely used for single or multiple traits, this can provide a greater sensitivity for environment-dependent QTL, inform on the stability of the QTL effect, but not on the environmental factor(s) the QTL is responding to. To overcome this issue, several studies compared two conditions, varying a single environmental factor and were able to study the genetic control of the plasticity to a specific stress for a particular trait (McKay et al., 2008; El-Soda et al., 2015). Several genes involved in GxE have been cloned in plants (reviewed by Des Marais et al., 2013) by using a combination of fine QTL mapping and candidate gene approaches. Most genes are involved in flowering time control (FT, PPD1, FLC, FRI, PHYC, CO). Examples concerning abiotic stress responses are still rare: CBF2 for cold stress (Alonso-Blanco et al., 2005), P5CS1 for osmotic stress (Kesari et al., 2012), SUB1A for submersion in rice (Fukao et al., 2011) and RAS1 for saline stress tolerance (Ren et al., 2010). No gene showing a GxE interaction for drought stress nor for any yield-related trait were cloned or fine-mapped yet. This certainly lies respectively in the difficulty to characterize drought stress in natural conditions and in the complex genetic architecture of yield-related traits, which implies small effect mutations. Indeed successful examples were obtained with easy-to-phenotype traits (flowering time) and easy-to-setup environmental factors controlled one at a time.

This reductionist vision of the environment prove to be efficient for gene cloning. However, more complex and realistic approaches are needed to understand plant and crop responses to environmental factors and therefore breed for trait plasticity ultimately providing stable crops. In the context of the phenomics era (Großkinsky et al., 2015), with a greater description of the environment, and with the genomic tools accessible on many species, these approaches should flourish. However, to our knowledge, no study tackled yet the identification of the genetic control of plant plasticity to combined environmental factors in nature or field trials (Mahalingam, 2015). Thanks to prior crop modeling, this approach is now amenable, and shall take advantage of a precise description of the plant stress factors, of the statistical power of multi-environment trials and of the realistic nature of the measured traits.

In this study, we developed a novel approach that combined crop modeling and quantitative genetics to identify the genetic basis of oil yield plasticity in sunflower. Our methodology consisted in estimating four abiotic stress indicators in a range of environments using a crop model and in selecting the indicators explaining the best the grain yield plasticity of commercial varieties. The plant-level stress estimations on these varieties were used to characterize each experimental sites for the three selected abiotic stresses. We could then estimate the plasticity of oil yield for each stress and for every hybrid of a sunflower diversity panel cultivated in every site. Therefore, as a demonstration on cold stress, we could successfully perform a genome-wide association study to identify genomic regions putatively involved in oil yield plasticity to cold.
# Material and methods

## Plant material

Association mapping was carried out on a panel of 317 inbred lines from INRA and sunflower breeding companies. This panel was a subset of the core collection of 384 inbred lines of Cadic et al. (2013) chosen for its diversity from an initial set of 752 inbred lines (Coque et al., 2008). It was comprised of both elite lines, parents of commercial hybrids, and lines with introgressions from several wild Helianthus accessions, including H. annuus, H. argophyllus and H. petiolaris.

Oil yield was observed on testcross progeny obtained by crossing panel lines with testers according to their status (maintainers of cytoplasmic male sterility [B-lines] or fertility restorers [R-lines]), as described in Cadic et al. (2013). R-lines were crossed with the two CMS PET1 counterparts of B-line proprietary testers (FS71501 or AT0521) while the B-lines were crossed with two R-line testers: 83HR4gms and SOLR001M. 83HR4gms is derived from the 83HR4 line and was converted to female by the introduction of a genetic male sterility, SOLR001M is a proprietary line carrying PEF1 cytoplasmic male sterility (Serieys, 1984; Crouzillat et al., 1991) which it maintains, although it is a restorer for classical PET1 cytoplasm (Table 1).

## Description of the multi-environment trial (MET)

From 2008 to 2010, eight locations located in the center and Southwest of France were planted with the testcross progeny. In six location, trials were conducted with and without irrigation, providing a total of 17 location x treatment x year combinations, (designated as environments). The panel lines were evaluated on the same tester in each environment (Table 1). Each experiment was an Augmented-Design (Federer, 1961) formed of blocks, with 24 or 30 entries replicated in two sub-blocks. Each sub-block was randomized separately and contained two to four control hybrids.

The climatic variability on experimental locations was summarized by computing the mean air temperature, the sum of water inputs (rainfall and irrigation) and the climatic water deficit (difference between precipitations, P and potential evapotranspiration, PET) on the cropping period (see supplementary Figure S1 and Table S1). No environment had a climatic water deficit (ranging from -177 to -458 mm), meaning that the climatic evaporative demand was always above the water supply (even accounting for irrigation). Rainfall on the cropping period ranged from a low 36 mm to 368.5 mm and the average amount of irrigation was 74 mm. Trials were performed on various soils depth leading to a soil water capacity (SWC) from 112 mm to 240 mm. Mean nitrogen fertilization was 60 kg/ha (eq. mineral nitrogen).

Among the 17 environments, we discarded 3 environments (AI09_I, AI09_NI, CO09_NI). The first two were outliers for the observed yield and the SUNFLO model failed at simulating yield phenotypes close to the observed one for the controls. The last one did not exhibit genotypic effect for the panel phenotypes, the estimated genotypic variance was judged non-significantly different from zero by a Z-ratio test in the linear mixed model used to correct for micro-environment effects and to predict the genotypic value of the yield panel lines.
## Intra environment phenotypic data analysis

Within each environment, the oil yield was adjusted for micro-environment effects using ASReml-R (Butler *et al.*, 2009), as described in Cadic *et al.* (2013). A linear mixed model (named naïve in Cadic *et al.*, 2012) with a random effect for the genotypic value of the panel lines, including blocks and sub-blocks as fixed effects was compared to two spatial models. The spatial models included (a) random effects of row and column or (b) a first-order autoregressive process in the residuals to take into account autocorrelation between neighbour plots. The three models were compared using the Akaike criterion (AIC) and the best linear unbiased predictors (BLUP) of the genotypic values were extracted from the best model according to AIC, for the next step of analysis.

## Estimation and choice of abiotic stresses in each environment using simulation

SUNFLO is a process-based model for the sunflower crop which was developed to simulate the grain yield and oil concentration as a function of time, environment (soil and climate), management practice and genetic diversity (Debaeke *et al.*, 2010; Casadebaig *et al.*, 2011; Lecoeur *et al.*, 2011). The model simulates the main soil and plant processes: root growth, soil water and nitrogen content, plant transpiration and nitrogen uptake, leaf expansion and senescence and biomass accumulation, as a function of main environmental constraints (temperature, radiation, water and nitrogen deficit).

This model is based on a conceptual framework initially proposed by Monteith (1977) and now shared by a large family of crop models (e.g. Jones *et al.*, 2003; Keating *et al.*, 2003; Brisson *et al.*, 2003). In this framework, the daily crop dry biomass ($DM_t$) is calculated as an ordinary difference equation (Equation 1) function of incident photosynthetically active radiation (PAR, MJ m⁻²), light interception efficiency (1-exp⁻$kLAI_t$) and radiation use efficiency (RUE, g MJ⁻¹, Monteith (1994)). The light interception efficiency is based on Beer-Lambert's law as a function of leaf area index ($LAI_t$) and light extinction coefficient ($k$).

$$DM_t = DM_{t-1} + RUE_t (1 - \exp^{-kLAI_t}) PAR_t \quad (1)$$

Broad scale processes of this framework, the dynamics of LAI, photosynthesis (RUE) and biomass allocation to grains were split into finer processes (e.g plant phenologic development, leaf expansion and senescence, response functions to environmental stresses) to reveal genotypic specificity and to allow the emergence of GxE interactions. Globally, the SUNFLO crop model has about 50 equations and 64 parameters (43 plant-related traits and 21 environment-related). When evaluated on the presented MET dataset, the SUNFLO model predicted accurately the performance of control hybrids across environments: the root of mean square error (RMSE) was 0.3 t ha⁻¹, relative RMSE was 9.6 %, bias was -0.14 t ha⁻¹.

Using the SUNFLO model, we computed two indicators (continuous and discrete) per type of considered abiotic stresses to characterize the different environments (Table 2). Each indicator was integrated over three key periods: vegetative stage (veg), flowering period (flo) and grain filling period (fil). We also considered the sum over two periods and during the whole cropping period, for a total of seven time periods per indicator. All these indicators corresponded to the mean stress felt by the control hybrids during the above periods.
For water stress, the fraction of transpirable soil water (FTSW) represents yield limitation through water deficit (integration of 1 minus FTSW); ETR is the conditional sum of days, if the ratio of the real evapotranspiration (ET) to potential evapotranspiration (PET) was less than 0.6 (threshold for photosynthesis limitation). For cold stress, LTs is the conditional sum of days if mean air temperature was below 20°C and LTi represent low temperatures impact on photosynthesis (integration of 1 minus equation 2). Heat stress indicators were computed following the same logic, albeit representing high temperatures impact on photosynthesis. Equations (2) and (3) are used in the crop model to define the radiation use efficiency (RUE) response to temperature (Villalobos et al., 1996). For nitrogen deficit, NAB is the amount of absorbed nitrogen in the considered cropping period and NNI is the sum of 1 minus nitrogen nutrition index, which indicates crop nitrogen deficit (Lemaire and Meynard, 1997; Debaeke et al., 2012).

$$LTRUE = \begin{cases} 0 & \text{if } T_m < T_b \\ \frac{T_m}{T_{ol} - T_b} - \frac{T_b}{T_{ol} - T_b} & \text{if } T_b < T_m < T_{ol} \\ 1 & \text{else} \end{cases} \quad (2)$$

$$HTRUE = \begin{cases} 1 & \text{if } T_m < T_{ou} \\ \frac{T_m}{T_{ou} - T_c} - \frac{T_c}{T_{ou} - T_c} & \text{if } T_{ou} < T_m < T_c \\ 0 & \text{else} \end{cases} \quad (3)$$

where $T_m$ denotes the mean daily air temperature (°C), $T_b$=4.8 is the base temperature (°C), $T_{ol}$=20 is the optimal lower temperature (°C), $T_{ou}$=28 is the optimal upper temperature (°C) and $T_c$=37 is the critical temperature (°C).

Three sunflower varieties (Melody, Pacific, Pegasol) which were used as controls within the environments had their phenotypic characteristics previously included in the SUNFLO model (Supplementary Table S2). The soil characteristics, the crop management and climate data were collected to allow simulation for each environment. One soil characteristic, the soil depth, is difficult to observe and has a strong impact on the simulated yield. Instead of using the approximated value given by the breeders and the farmers, we adjusted this parameter by minimizing the empirical mean square error between the observed and the simulated yields.

Model choice to select the stress indicators was made with the AIC using the native R function lm. AIC of each model was computed for all three control varieties and model choice was made on the mean over the three controls. Compared linear models, which fitted the grain yield with the indicators, were limited to combinations of only one indicator per type of stresses. All models integrating one, two, three or four stresses were compared; in total 50,624 models were computed and compared.

The R function lm was also used to compute the p-value of the Fisher test for each indicator including in the best model.
# Genetic study

## Estimation of plasticities of oil yield to water, nitrogen and cold stresses in the diversity panel

In order to get plasticity phenotypes that reflect the responses of the panel lines to the different abiotic factors, we adjusted the BLUP phenotype of each panel line with the following linear model:

$$Y_{ij} = a_i + b_i DS_j + c_i CS_j + d_i NS_j + \varepsilon_{ij} \quad (4)$$

where $Y_{ij}$ is the BLUP of the phenotype for the $i$th genotype, $a_i$ the potential phenotype in an environment with average stresses, $DS_j$ the water stress indicator in the $j$th environment calculated as the mean over the 3 controls, $b_i$ the slope linked to the water stress, $CS_j$ the cold stress indicator in the $j$th environment calculated as the mean over the 3 controls, $c_i$ the slope linked to the cold stress, $NS_j$ the nitrogen stress indicator in the $j$th environment calculated as the mean over the 3 controls, $d_i$ the slope linked to the nitrogen stress and $\varepsilon_{ij}$ the residual variance. The covariates $DS$, $CS$, and $NS$ were centered before computing the regression model. The three estimated covariate coefficients $\hat{b}_i$, $\hat{c}_i$, and $\hat{d}_i$ of this model are the plasticity phenotypes of interest, i.e. the genotype slopes in response to water, cold and nitrogen stresses respectively.

To compare model (Error! Reference source not found.) to a more simple regression model, we also fitted the BLUP phenotype of each panel line using the cold stress indicator as a single regressor.

Before computing the above linear models, the oil yield missing data were imputed using the missMDA R package (Josse et al., 2012; Josse and Husson, 2016). All recorded traits (24 traits for a total of 198 traits x environments) were used to impute missing oil yield values. However, panel lines with less than six observations per oil yield trait over the environments were discarded.

## Combined stresses in a single multi-stress index of plasticity

In order to compare the panel lines for the stability against multiple stresses together, we defined a multi-stress plasticity index accounting for the variance-covariance of the stress slopes, equals to:

$$(\hat{b}_i, \hat{c}_i, \hat{d}_i) V^{-1} (\hat{b}_i, \hat{c}_i, \hat{d}_i)^t \quad (5)$$

where $V$ denotes the (3 by 3) variance-covariance matrix.

To visualize the strategy of the panel lines versus the three stresses, we drew a star representation using the star function of the R package graphics. The three plasticity phenotypes were first taken in absolute value and scaled to 0-1 in order to have three comparable values with values close to 0 noted the stability and values close to 1 noted the instability. Then they were normalized by panel line to sum to 1 in order to get comparable values for all panel lines. The star representation was applied on these scaled and normalized values that can be interpreted as the percentage of stability dedicated to each stress.
## Genotyping and map building.

A set of 197,914 SNP were used to produce an AXIOM® genotyping 96-array (Affymetrix, Santa Clara, CA, USA). These SNP were selected from either genomic re-sequencing or transcriptomic experiments. An additional set of 6,800 non-polymorphic sequences were added as controls. Combined with internal technical controls, the AXIOM® genotyping 96-array was designed with a total of 445,876 probesets.

Genomic DNA from the 317 panel lines and two recombinant inbred lines (RIL) populations INEDI and FUxPAZ2 obtained from the cross between XRQ and PSC8 lines (180 RIL) and from the cross between FU and PAZ2 lines (87 RIL) respectively, were genotyped with the AXIOM® array. All hybridization experiments were performed by Affymetrix and the genotypic data were obtained with the GTC software (Affymetrix). From the 197,914 SNP, 35,562 were polymorphic between XRQ and PSC8 and 28,529 between FU and PAZ2.

We used CarthaGène v1.3 (Givry et al., 2005) to build the genetic maps of the INEDI population and the FUxPAZ2 population separately. For the INEDI population, we added to the set of AXIOM® SNP the markers previously mapped by Cadic et al. (2013) which allowed completing this map and assigning AXIOM® markers to appropriate linkage group (LG). We built a consensus map with common markers of both previous maps using Biomercator v4.0 (Sosnowski et al., 2012) and we projected the specific markers of the two previous maps on this consensus map. Unmapped AXIOM® SNP were placed by BLAST analysis on the RHA280xRHA801 genetic map based on 454 sequencing (Kane et al., 2011) and finally projected on the INEDI and FUxPAZ2 consensus map. The remaining AXIOM® SNP were located by computing the linkage disequilibrium (LD) measurements proposed by Mangin et al. (2012). They were mapped to the same position of the mapped SNP in maximum LD. For this, we used as LD statistics the maximum of the $r_V^2$ and $r_{VS}^2$ measurements that correct for relatedness and for both structure and relatedness, respectively.

## Association tests

Association mapping was based on a set of 65,534 SNP with MAF > 0.05. Similarly to our previous work (Cadic et al., 2013), two association models were performed using EMMA (Kang et al., 2008) based on the Yu et al. model (2006). Both models included a correction for the genomic relatedness using the alike-in-state (AIS) kinship estimated with EMMA version v1.1.2 R package (Villanova et al., 2011) using all the above SNP. The population structure modelled as the restorer or maintainer status of the panel lines was added in the second model leading to the following model:

$$Sl_i = \sum_{c} X_{ic} \alpha_c + M_{il} \Theta_l + u_i + e_i \quad (6)$$

$Sl_i$ is the slope for ith line, $X_{ic}$ is the line status, $\alpha_c$ is the effect of the line status c, $M_{il}$ is the genotype of the ith line at locus l, $\Theta_l$ is the effect of locus l. $\alpha_c$ and $\Theta_l$ are considered to be fixed effects. $u_i$ is the random polygenic effect modelling genetic relatedness with $Var(u)=\sigma_u^2 K_{ais}$ where $K_{ais}$ is an AIS matrix and $Var(e)=\sigma_e^2 I$ where I denotes the identity matrix. Multiple testing correction was achieved using an approximate effective number of
tests ($M_{eff}$) based on the eigen-values of the SNP correlation matrix as proposed by Li and Ji (2005). We computed $M_{eff}$ by using blocks of 250 markers and assuming that the blocks were independent.

In order to conduct a multi-loci analysis, the genotypic data was imputed using Beagle v4.0 (Browning and Browning, 2007) and MLMM (Segura et al., 2012) was used for this purpose. We stopped the forward approach of MLMM when the variance of the polygenic term was non significant, which is a little different compared to the initial forward procedure of MLMM that stops when the estimator of this polygenic variance is equal to 0. The significance of the polygenic variance component was judged with a log-likelihood ratio test and a risk of 1%. This log-likelihood ratio test compared the model with and without the polygenic effect. Two times the difference between the log-likelihood of the models is known to follow asymptotically a mixture between a Dirac at 0 and a Chi-square distribution with 1 degree of freedom (Self and Liang, 1987). When the forward approach stopped, a SNP was judged associated if its Bonferonni corrected p-value, using $M_{eff}$ as the number of independent tests, was inferior to the chosen type I error. The associated SNP were put all together with the polygenic term to create the final multi-loci linear model. The associated SNP effect, their reevaluated p-values were computed with the base R function lm in this final model (R Core Team, 2014).

## Functional annotation of associated SNP

Context sequences of associated SNP were compared to the sunflower genome (line XRQ) sequenced using the PacBio (Pacific Biosciences, Menlo Park, CA, USA) technology (https://www.heliagene.org/HanXRQ-SUNRISE/) by blast analysis. In case of ambiguous positioning on the genome, we retained the chromosomal position in accordance to genetic map location of the marker. Similarly, we positioned associated SNP on the reference transcriptome (https://www.heliagene.org/HaT13l/).

Accepted for publication [This article is protected by copyright. All rights reserved.
# Results

## Estimation of abiotic stresses in field environment.

The characterization of abiotic stresses at the plant level through modelling and simulation explained observed yield variability better than climate-based indicators (Figure 1). For example, the correlation between yield observed in field experiments and water deficit computed from (1) climate data only (P - PET), (2) climate, soil, and management data (P - PET + SWC + irrigation) and (3) simulated plant data (ETR, defined in Table 2) was gradually stronger (respectively -0.01, -0.65 and -0.86). The correlation between yield and water deficit index increased in strength with the "proximity" of the regressor to the plant, revealing the expected negative impact of water deficit on crop yield.

The correlation between different type of abiotic stresses (Figure 2), indicated that cold (LT) and heat stresses (HT) were naturally the most highly negatively correlated combination. The nitrogen indicator NAB was not correlated to other indicators except the negative correlation with the water stress indicator (ETR, particularly during grain filling). The other nitrogen indicator NNI, which was not correlated to NAB showed both a positive correlation with heat stress indicators and a negative correlation with cold indicators. The water stress indicator FTSW during grain filling was also positively correlated to heat stress indicators (HT) and nitrogen stress (NNI), and negatively correlated to cold stress indicators. Within a single type of stress, correlations between indicators computed during different cropping periods were generally positive and high (Figure 2), ranging 0.40 for water deficit (FTSW) to 0.98 for heat stress (HTi), with nitrogen and cold indicators in-between.

<figure><img src="image_3.png" /></figure>
## Model selection to estimate the best combination of abiotic stresses.

The model selected as having the best AIC on average over the control genotypes was a three indicator model, including a cold stress indicator during the vegetative stage (LTi<sub>veg</sub>), a water stress indicator during the vegetative and the flowering period (FTSW<sub>veg+flo</sub>), and a nitrogen indicator during the whole growth period (NAB). Table 3 presents the p-value of the Fisher test for each control in this best model and the proportion of variance explained by this best model for each control.

Each stress indicator had a significant effect on the grain yield for each control genotype. The highest impact was due to cold stress for Melody (p-value = 2.39 10<sup>-4</sup>) and the smallest was due to water stress for Pacific (p-value = 2.91 10<sup>-2</sup>) as indicated in Table 3. All together, these indicators explained very well the grain yield variability of each control since the percentages of explained variance were 90%, 93% and 94% for Pegasol, Pacific and Melody respectively.

Although cold stress has a strong impact on yield (-0.35 q h<sup>-1</sup> d<sup>-1</sup>), the relatively low number of day of cold stress in the MET (8.7 on average, variation in the MET shown on Figure 3) reduces its impact, whereas water stress impacts the most yield due to the high number of days of stress (33.9 on average) even if sunflower is relatively tolerant (-0.14 q h<sup>-1</sup> d<sup>-1</sup>).

## Plasticity in the diversity panel.

As we described the actual cold, nitrogen and water stresses felt by sunflower in 14 environments, we also measured oil yield in a diversity panel in those environments. This allowed us to compute plasticities of oil yield to abiotic stresses for each line in the panel as the slopes of regression of oil yield to individual stress indicators. These three slopes represented how plastic was the response of a line faced to a given stress and are therefore referred as plasticity phenotypes later. Minimum, mean, maximum, and variance of these plasticities as well as their correlation are presented in Table 4 and their distribution histograms are in supplementary Figure S2. Cold stress plasticity appears genetically independent to water and nitrogen stress plasticities. On the contrary, nitrogen and water stress plasticity are highly correlated (Pearson correlation of 0.62) suggesting common genetic control of these traits in our panel.

To illustrate the importance of multi-stress indicator modelling, we compared regression of oil yield in a single stress indicator model (including only cold stress indicator) and in the multi-stress indicator model in Figure 4. Oil yields of the most sensitive and tolerant lines were plotted against the cold stress indicator. Point clouds are closer to the regression lines in the multi-stress indicator model indicating a better characterization of cold stress impact in this modelling approach.
## A multi-stress plasticity index.

Following the estimation of the three single stress plasticities, we were interested in calculating a multi-stress plasticity index to describe the general abiotic tolerance of every line. This index is a weighted sum of the three abiotic stress plasticities taking into account the correlation between them, specifically the one between drought and nitrogen. This index allowed us to rank the panel lines and to describe the different strategies observed in the panel to tolerate combined stresses. This is illustrated on Figure 5 that shows the multi-stress plasticity index of every panel lines against its mean oil yield in the MET. Stable panel lines, with a small multi-stress plasticity index, as well as panel lines belonging to the border of the point cloud were highlighted by a star representation (a triangle in our case) illustrating their abiotic stress tolerance strategy. First, we observed that unstable panel lines were generally more sensitive to cold stress compared to stable lines. This was confirmed by a significant correlation between the multi-stress plasticity index and the percentage of plasticity due to the cold stress (p-value of 6 *10⁻⁴) although this was not observed for nitrogen or water stress (p-value of 0.28 and 0.10 respectively). Second, we could not identify a consensus strategy for stable panel lines. As examples: (i) the three similarly stable panel lines (multi-stress index around 45) having high oil yield (around 15) developed tolerance to all stresses but showed different plasticity patterns, (ii) the two most stable panel lines had opposite profile of stability with the most stable being sensitive to cold and the second to water.

## Genetic map.

In order to position the genomic regions controlling the abiotic stress plasticity, we constructed a genetic map with the markers genotyped on the diversity panel using three RIL populations. The final genetic map was composed of 89,979 markers positioned in 4,782 genetic positions for a total distance of 1,398.5 cM. On this map, 4,094 markers were used to build the consensus map between the INEDI (XRQxPSC8) and FUxPAZ2 RIL populations. Among the other markers, 27,663 were mapped using the INEDI population alone, 13,807 were mapped using the FUxPAZ2 population alone, 29,586 were located thanks to a genotyping by sequencing map on the RIL population RHA801xRHA280 (Kane et al., 2011) and the remaining 14,829 markers were placed by linkage disequilibrium. Details of these maps can be found in the supplementary Table S3.

## Association study of oil yield plasticity to cold stress.

For the genetic analysis, we focused, as a proof of principle of the approach, on the plasticity to cold stress as it is the most impacting stress and appeared genetically independent from the others. Among the 65,534 association tests, the effective number of tests (Li and Ji, 2005) was estimated around 14,000. Using this effective number of tests, we kept SNP associated with a Bonferroni corrected p-value of 7.1 *10⁻⁶ for a family wise type I error of 10%.

Only 2 SNP were detected by a marker by marker association analysis using EMMA. The most significant SNP was located at the end of LG 5 in a QTL named LG05.64, and was detected with a model including the maintainer
or restorer status as a fixed structure effect. The second is located at the center of LG 17 and was detected without the line status effect.

We completed this study by a forward approach of the multi-loci association analysis (MLMM) for both models (with or without the maintainer/restorer status). Both models stopped with 6 SNP among which four were judged associated. None of them was in common between the two models but one SNP corresponded to the previously identified SNP in LG05.64. The most associated SNP was located at the end of the LG13.

In total 9 QTL could be identified using MLMM and EMMA procedures with or without population structures: two located on LG 5 and 10 and one on LG 9, 13, 14, 16 and 17. Their phenotypic effects on cold stress plasticity varies from 10 to 21% of the average plasticity in the panel (Table 5).

## Genes located in QTL controlling oil yield plasticity to cold stress.

We were interested in genes containing associated SNP to link the genetic identification to molecular and physiological processes putatively involved in cold tolerance. All associated SNP were located within coding sequences as expected from the AXIOM® genotyping array design. Functional annotation of the corresponding genes pointed out homologues of: NPF3.1, LTP, CYS6, NPF5.3, GMII, RPD1, PPX1, HAOX2 and IAR4 (from the most to the least significantly associated, as shown in Table 5). Strikingly, two close homologues of oligopeptide transporters (NPF3.1 and NPF5.3) are present in associated QTL on LG 5 and 14 and two homologues of genes involved in root development (RPD1 and IAR4) on LG 9 and 17. In addition, homologues of a lipid transfer protein, a cystatin, an alpha-mannosidase, a protein phosphatase and an aldolase are also in QTL controlling oil yield plasticity to cold stress.

Accepted
# Discussion

In this work, we developed a novel method to characterize the abiotic stress levels on different environments by using crop modelling and simulation and subsequently exploited it to identify genetic control of stress plasticity. We implemented this environment characterisation method on 17 locations from a multi-environment trial for sunflower and four abiotic stresses (water, nitrogen, cold and heat). The SUNFLO model (Casadebaig et al., 2011) was used to simulate stress patterns dynamics for three varieties used as controls in each location and integrated indicators were computed from these data considering different crop phenological stages and physiological processes (eight stress indicators over seven periods). We used a model selection approach to select the best linear model among combinations of stress indicators used as regressors for yield. Water, nitrogen and cold stresses were retained as the most explicative abiotic stresses for yield variability in this MET. Using reaction norms as conceptual reference, we computed abiotic stress plasticities as the slopes of the linear regression of oil yield on selected stress indicators. We then conducted an association study with a panel of 317 lines genotyped for nearly 65,000 markers on oil yield plasticity for cold stress as its was the most impacting stress (per time unit) and was not correlated to other abiotic stress plasticities.

## Crop modelling helped to analyze abiotic stress patterns and to explain their impact on yield.

In a location, stress indicators can be climate-based (precipitations minus evapotranspiration), crop-based (simple water balance including soil water capacity and irrigation) or plant-based (simulated dynamic water balance). We observed that the observed grain yield was best explained by plant-based stress indicators because the interactions between climate, leaf area dynamics, plant stomatal conductance (isohydric vs anisohydric behaviours, e.g. Casadebaig et al. (2008) for sunflower) and management practices can be partly reproduced by the crop model algorithm. This is for example illustrated by environment rankings, where some irrigated locations (CO08_I, CO09_I) still display a high level of water stress while a rainfed one (CA10_NI, VE09_NI) show reduced water deficit.

Abiotic stresses also do not have the same impact on crop physiology according to their timing of occurrence during the crop cycle (Table 3). Among the seven possible combinations between main crop phases (i.e. vegetative, flowering, grain filling), we indeed observed that the relevance of these timings was specific to the type of abiotic stress. For cold stress, the detection of early crop growth (vegetative period) was expected because this stress occurrence is strongly determined by the climate (low temperature during crop installation). However, in continental climates, where sunflower is mainly grown, we can also observe cold temperatures at the end of the crop cycle. For water stress, where interactions between crop growth and climate variability are more important, vegetative and flowering periods were identified, which is consistent with numerous previous reports on sunflower (Blanchet et al., 1990; Cabelguenne et al., 1999). Regarding nitrogen stress, the importance of this process over the whole crop cycle was highlighted. Indeed, recent reports indicated that post-flowering nitrogen absorption could also be significant (Andrianasolo et al., 2016). Remarkably, heat stress was not identified as a major contribution to yield variability in the MET. Actually, depreciative effect of high temperatures on photosynthesis were caused by temperatures that were almost never reached in our experimental conditions. According to the current parameterization of the simulation model, heat stress indicator (Error! Reference
source not found.) was only significant in one location (CO09) and null (9/17 locations) or weak in the others (Figures 2 and S1).

Genetic control of plasticity of oil yield to cold stress

The most and fourth most associated SNP pointed to two genes highly homologous to NPF3.1 and NPF5.3 that are both oligopeptide transporters. These two independent association signals on chromosomes 5 and 14 strongly suggest a role of oligopeptide transport in tolerance to cold stress observed in our experimental conditions i.e. when young plants are exposed to chilling. In plants, these transporters are key players in nitrogen nutrition and therefore plantlet growth. The importance of oligopeptide transport to tolerate cold is corroborated by the demonstrated molecular adaptation of this transporter family in antarctic icefish (*Chionodraco hamatus*) adapted to sub-zero temperatures (Maffia *et al.*, 2003; Rizzello *et al.*, 2013). The role in N nutrition of these transporters in animal and plants indicates that nutrient transport can be a limiting factor at low temperature that likely limits remobilization of seed stocks and/or absorption and transport of N from roots to aerial organs.

The second most associated SNP is located in a putative lipid transfer protein (LTP) (QTL LG16.48). Many LTP have been reported to be transcriptionally induced by freezing (reviewed in Liu *et al.*, 2015) and over-expression of LTP3 provided freezing tolerance in *A. thaliana* (Guo *et al.*, 2013). This action could be due to membrane stabilization as demonstrated by Hincha *et al.* (2001) in preventing chloroplastic damages induced by freezing, or through its role in seed lipid mobilization during germination and seedling growth (Pagnussat *et al.*, 2015) as shown in sunflower for another LTP (Pagnussat *et al.*, 2009). Another candidate genes is the cystatin CYS6 homologue located in QTL LG13.72. Several homologues of cystatin were shown to be induced during cold exposure in barley (Gaddour *et al.*, 2001), maize (Massonneau *et al.*, 2005), wheat (Talanova *et al.*, 2012), and increase, when over-expressed, cold tolerance in Arabidopsis (Zhang *et al.*, 2008). Based on Prins *et al.* (2008), the sunflower cystatin could provide a better regulation of Rubisco turnover in chloroplasts in cold conditions. The alternative oxidase HAOX2 homologue found in QTL LG10.34 constitutes another good candidate gene for cold tolerance. The Alternative Oxidase Pathway (AOP) has been described in many plants to be involved in cold stress response as a biochemical protection against overproduction of Reactive Oxygen Species due to the cold inhibition of the electron transport chain in mitochondria (Feng *et al.*, 2008). Furthermore, the AOP was shown to participate in differential cold-sensitivity between two maize genotypes (Ribas-Carbo *et al.*, 2000).

Interestingly, two genes (homologues to RPD1, IAR4) were found (QTL LG09.27, LG17.49) and their Arabidopsis counterparts share similar features: both are involved in root development and both mutants show temperature-sensitive phenotypes (at 20°C and 28°C) (Konishi and Sugiyama, 2006; Quint *et al.*, 2009). This suggests that root setting could also be temperature- and genotype-dependent in sunflower. All together, the functional annotation of QTL associated to cold stress plasticity of oil yield identified several candidate genes and physiological processes. Most of them were already described to be involved in cold stress tolerance in other plants which supports our study and indicates a probably short genetic distance between associated SNP and causal mutations.

To complete our understanding of how these processes that act during the early growth phase of sunflower, impact the final seed yield, further studies putting in relation dynamic measurements of plant growth rate and cell physiology would be enlightening.
# Phenotypic plasticity and tradeoff for potential yield

In the studied multi-environment network, mean oil yield was positively correlated with a high sensitivity to environmental stresses, indicating that a global gain in performance was generally associated to a higher yield instability (Figure 5). This is globally in accordance with previous claims that an increase in phenotypic plasticity allowed to achieve better yield stability across environments but at the expense of greater performance in low stress conditions (Sadras *et al.*, 2009). However the presence of both high-yielding and stress-tolerant genotypes suggests this general observation can be genetically by-passed and leaves room for more efficient sunflower varieties.

Stability of complex traits such as yield or fitness depends on plasticity of numerous intermediate traits (likely physiological and developmental processes) that are yet unknown. In our approach, we studied directly the plasticity of the complex trait with the idea of stabilizing it. The molecular and physiological processes pointed out by the genetic analysis allow us to identify some candidate processes: oligopeptide transport, root development, ROS scavenging, chloroplast and mitochondrial physiology (i.e. intermediate traits). In this kind of approach, we can wonder whether we could detect key regulators such as transcription factors (none in our case). Indeed, genetic variation in those would likely have trade-off effects on various physiological processes. Then, they would impact intermediate traits significantly but with possibly opposite effects and at the end, no significant impact on the resulting complex trait.

On the breeding strategy point of view, the considered environmental stresses (temperature, nitrogen and water) do not necessarily coexist in the French target population of environments: e.g. south-western production regions are exposed to drought and heat stress while the temperatures in northern regions are low enough and necessitate short crop cycle cultivars. This lack of spatial superposition of environmental constraints allow to exploit the differential sensitivities in the studied genetic material and adapt cultivar choice or breeding according to local growing environment. On the other hand, breeding for phenotypic plasticity of intermediate traits would potentially result in resilience to increasingly unpredictable environments (Nicotra *et al.*, 2010).

# Conclusion

Improving crop performance in low-input cropping systems requires a coordinated improvement of genotypes and agronomical practices (Sadras and Denison, 2016). In these growth-limiting conditions, abiotic stresses occur in combined and dynamic patterns. Therefore, disentangling those using models allows to understand their specific impacts on complex traits (such as yield) and the genetic factors potentially reducing those.

In our study, precise characterization of water, cold, heat and nitrogen stresses allowed accurate identification of nine QTL and underlying genes controlling stress plasticity. This joint approach between crop modelling and quantitative genetics also permitted estimations of allelic variation in natural conditions.

Such inter-disciplinary approach should be useful to conduct different breeding strategies and adapt crop to climate variability through local adaptation: maximizing performance in a given environment type, or through global adaptation: maximizing yield stability over different environment types.
Acknowledgement

This work benefited from the GENOPLANTE program “HP1” (2001–2004), the “SUNYFUEL” project, financially supported by the French National Research Agency (ANR-07-GPLA-0022, 2008–2011), the “OLEOSOL” project (2009–2012) with the financial support from the Midi Pyrénées Region, the European Fund for Regional Development (EFRD), and the French Fund for Competitiveness Clusters (FUI), and the “SUNRISE” project of the French National Research Agency (ANR-11-BTBR-0005, 2012-2019).

Accepted Article
References

Acosta-Gallegos J, White JW. 1995. Phenological plasticity as an adaptation by common bean to rainfed environments. *Crop Science* **35**, 199–204.

Alonso-Blanco C, Gomez-Mena C, Llorente F, Koornneef M, Salinas J, Martínez-Zapater JM. 2005. Genetic and molecular analyses of natural variation indicate CBF2 as a candidate gene for underlying a freezing tolerance quantitative trait locus in arabidopsis. *Plant Physiology* **139**, 1304–1312.

Andrianasolo FN, Champolivier L, Debaeke P, Maury P. 2016. Source and sink indicators for determining nitrogen, plant density and genotype effects on oil and protein contents in sunflower achenes. *Field Crops Research*.

Blanchet R, Texier V, Gelfi N, Viguier P. 1990. Articulation des divers processus d'adaptation à la sécheresse et comportements globaux du tournesol. In: CETIOM, ed. *Le tournesol et l'eau : Adaptation à la sécheresse, réponse à l'irrigation*. Paris, France, 45-55.

Brisson N, Gary C, Justes E et al. 2003. An overview of the crop model STICS. *European Journal of Agronomy* **18**, 309-332.

Browning SR, Browning BL. 2007. Rapid and accurate haplotype phasing and missing-data inference for whole-genome association studies by use of localized haplotype clustering. *The American Journal of Human Genetics* **81**, 1084-1097.

Butler D, Cullis B, Gilmour A, Gogel B. 2009. ASReml-r reference manual. Queensland Department of Primary Industries, Queensland, Australia.

Cabelguenne M, Debaeke P, Bouniols A. 1999. EPICphase, a version of the EPIC model simulating the effects of water and nitrogen stress on biomass and yield, taking account of developmental stages: Validation on maize, sunflower, sorghum, soybean and winter wheat. *Agricultural Systems* **60**, 175-196.

Cadic E, Coque M, Vear F et al. 2013. Combined linkage and association mapping of flowering time in sunflower (helianthus annuus L.). *Theoretical and Applied Genetics* **126**, 1337-1356.

Casadebaig P, Debaeke P, Lecoeur J. 2008. Thresholds for leaf expansion and transpiration response to soil water deficit in a range of sunflower genotypes. *European Journal of Agronomy* **28**, 646-654.

Casadebaig P, Guilioni L, Lecoeur J, Christophe A, Champolivier L, Debaeke P. 2011. SUNFLO, a model to simulate genotype-specific performance of the sunflower crop in contrasting environments. *Agricultural and Forest Meteorology* **151**, 163-178.

Chenu K, Deihmfard R, Chapman SC. 2013. Large-scale characterization of drought pattern: A continent-wide modelling approach applied to the Australian wheatbelt-spatial and temporal trends. *New Phytologist* **198**, 801-820.

Coque M, Mesnildrey S, Romestant M, Grezes-Besset B, Vear F, Langlade N, Vincourt P. 2008. Sunflower line core collections for association studies and phenomics. *Proc. 17th Int Sunflower Conf, Cordoba (Spain)*, 725-8.
Crouzillat D, Canal L, Perrault A, Ledoigt G, Vear F, Serieys H. 1991. Cytoplasmic male sterility in sunflower: Comparison of molecular biology and genetic studies. *Plant molecular biology* **16**, 415-426.

Debaeke P, Casadebaig P, Haquin B, Mestries E, Palleau J-P, Salvi F. 2010. Simulation de la réponse variétale du tournesol à l'environnement à l'aide du modèle SUNFLO. *Oléagineux, Corps Gras, Lipides* **17**, 143-51.

Debaeke P, Oosterom E van, Justes E, Champolivier L, Merrien A, Aguirrezabal L, González-Dugo V, Massignam A, Montemurro F. 2012. A species-specific critical nitrogen dilution curve for sunflower (helianthus annuus L.). *Field Crops Research* **136**, 76-84.

Des Marais DL, Hernandez KM, Juenger TE. 2013. Genotype-by-environment interaction and plasticity: Exploring genomic responses of plants to the abiotic environment. *Annual Review of Ecology, Evolution, and Systematics* **44**, 5-29.

DeWitt T, Langerhans R. 2004. Integrated solutions to environmental heterogeneity: Theory of multimoment reaction norms. *Phenotypic Plasticity: Functional and Conceptual Approaches*, 98-111.

El-Soda M, Kruijer W, Malosetti M, Koornneef M, Aarts MGM. 2015. Quantitative trait loci and candidate genes underlying genotype by environment interaction in the response of arabidopsis thaliana to drought. *Plant, Cell & Environment* **38**, 585-599.

Feng H, Li X, Duan J, Li H, Liang H. 2008. Chilling tolerance of wheat seedlings is related to an enhanced alternative respiratory pathway. *Crop Science* **48**, 2381.

Fukao T, Yeung E, Bailey-Serres J. 2011. The submergence tolerance regulator SUB1A mediates crosstalk between submergence and drought tolerance in rice. *The Plant Cell* **23**, 412-427.

r Gaddour K, Vicente-Carbajosa J, Lara P, Isabel-Lamoneda I, Díaz I, Carbonero P. 2001. A constitutive cystatin-encoding gene from barley (icy) responds differentially to abiotic stimuli. *Plant Molecular Biology* **45**, 599-608.

Givry S de, Bouchez M, Chabrier P, Milan D, Schiex T. 2005. CAR(H)(T)AGene: Multipopulation integrated genetic and radiation hybrid mapping. *Bioinformatics* **21**, 1703-1704.

Großkinsky DK, Svensgaard J, Christensen S, Roitsch T. 2015. Plant phenomics and the need for physiological phenotyping across scales to narrow the genotype-to-phenotype knowledge gap. *Journal of Experimental Botany* **66**, 5429-5440.

Guo L, Yang H, Zhang X, Yang S. 2013. Lipid transfer protein 3 as a target of MYB96 mediates freezing and drought stress in arabidopsis. *Journal of Experimental Botany* **64**, 1755-1767.

Jones JW, Hoogenboom G, Porter C, Boote K, Batchelor W, Hunt L, Wilkens P, Singh U, Gijsman A, Ritchie J. 2003. The dSSAT cropping system model. *European Journal of Agronomy* **18**, 235-265.

Josse J, Husson F. 2016. missMDA: A package for handling missing values in multivariate data analysis. *Journal of Statistical Software* **70**, 1-31.
Josse J, Chavent M, Liquet B, Husson F. 2012. Handling missing values with regularized iterative multiple correspondence analysis. *Journal of Classification* **29**, 91–116.

Kane NC, Gill N, Lindhauer MG *et al*. 2011. Progress towards a reference genome for sunflower. *Botany* **89**, 429–437.

Kang HM, Zaitlen NA, Wade CM, Kirby A, Heckerman D, Daly MJ, Eskin E. 2008. Efficient control of population structure in model organism association mapping. *Genetics* **178**, 1709–1723.

Keating BA, Carberry PS, Hammer GL, Probert ME, Robertson MJ, Holzworth D, Huth NI, Hargreaves JNG, Meinke H, Hochman Z. 2003. An overview of APSIM, a model designed for farming systems simulation. *European Journal of Agronomy* **18**, 267–288.

Kesari R, Lasky JR, Villamor JG, Marais DLD, Chen Y-JC, Liu T-W, Lin W, Juenger TE, Verslues PE. 2012. Intron-mediated alternative splicing of arabidopsis p5CS1 and its association with natural variation in proline and climate adaptation. *Proceedings of the National Academy of Sciences* **109**, 9197–9202.

Kiani M, Gheysari M, Mostafazadeh-Fard B, Majidi MM, Karchani K, Hoogenboom G. 2016. Effect of the interaction of water and nitrogen on sunflower under drip irrigation in an arid region. *Agricultural Water Management* **171**, 162–172.

Konishi M, Sugiyama M. 2006. A novel plant-specific family gene, ROOT PRIMORDIUM DEFECTIVE 1, is required for the maintenance of active cell proliferation. *Plant Physiology* **140**, 591–602.

Lake L, Chenu K, Sadras V. 2016. Patterns of water stress and temperature for australian chickpea production. *Crop and Pasture Science*.

Lecoeur J, Poiré-Lassus R, Christophe A, Pallas B, Casadebaig P, Debaeke P, Vear F, Guilioni L. 2011. Quantifying physiological determinants of genetic variation for yield potential in sunflower. SUNFLO: a model-based analysis. *Functional Plant Biology* **38**, 246–259.

Lemaire G, Meynard J. 1997. Use of the nitrogen nutrition index for the analysis of agronomical data. In: *Diagnosis of the Nitrogen Status in Crops*. Ed. G. Lemaire. Springer-Verlag, Berlin, 45–55.

Li J, Ji L. 2005. Adjusting multiple testing in multilocus analyses using the eigenvalues of a correlation matrix. *Heredity* **95**, 221–227.

Liu F, Zhang X, Lu C, Zeng X, Li Y, Fu D, Wu G. 2015. Non-specific lipid transfer proteins in plants: Presenting new advances and an integrated functional analysis. *Journal of Experimental Botany* **66**, 5663–5681.

Lobell DB, Burke MB, Tebaldi C, Mastrandrea MD, Falcon WP, Naylor RL. 2008. Prioritizing climate change adaptation needs for food security in 2030. *Science* **319**, 607–610.

Maffia M, Rizzello A, Acierno R, Verri T, Rollo M, Danieli A, Döring F, Daniel H, Storelli C. 2003. Characterisation of intestinal peptide transporter of the antarctic haemoglobinless teleost chionodraco hamatus. *Journal of Experimental Biology* **206**, 705–714.

Mahalingam R. 2015. Consideration of combined stress: A crucial paradigm for improving multiple stress tolerance in plants. In: Mahalingam R, ed. *Combined stresses in plants*. Springer International Publishing, 1–25.
Mangin B, Siberchicot A, Nicolas S, Doligez A, This P, Cierco-Ayrolles C. 2012. Novel measures of linkage disequilibrium that correct the bias due to population structure and relatedness. *Heredity* **108**, 285–291.

Massonneau A, Condamine P, Wisniewski JP, Zivy M, Rogowsky PM. 2005. Maize cystatins respond to developmental cues, cold stress and drought. *Biochimica Et Biophysica Acta-Gene Structure and Expression* **1729**, 186–199.

McKay J, Richards J, Nemali K, Sen S, Mitchell-Olds T, Boles S, Stahl E, Wayne T, Juenger T. 2008. Genetics of drought adaptation in arabidopsis thaliana II. QTL analysis of a new mapping population, kAS-1 x tSU-1. *Evolution* **62**, 3014.

Monteith JL. 1977. Climate and the Efficiency of Crop Production in Britain. *Philosophical Transactions of the Royal Society of London. Series B, Biological Sciences* **281**, 277–294.

Monteith JL. 1994. Validity of the correlation between intercepted radiation and biomass. *Agricultural and Forest Meteorology* **68**, 213–220.

Nicotra AB, Atkin OK, Bonser SP et al. 2010. Plant phenotypic plasticity in a changing climate. *Trends in Plant Science* **15**, 684–692.

Pagnussat LA, Lombardo C, Regente M, Pinedo M, Martin M, Canal L de la. 2009. Unexpected localization of a lipid transfer protein in germinating sunflower seeds. *Journal of Plant Physiology* **166**, 797–806.

Pagnussat LA, Oyarburo N, Cimmino C, Pinedo ML, Canal L de la. 2015. On the role of a lipid-transfer protein. arabidopsis ltp3 mutant is compromised in germination and seedling growth. *Plant Signaling & Behavior* **10**, e1105417.

Prins A, Van Heerden PDR, Olmos E, Kunert KJ, Foyer CH. 2008. Cysteine proteinases regulate chloroplast protein content and composition in tobacco leaves: A model for dynamic interactions with ribulose-1,5-bisphosphate carboxylase/oxygenase (rubisco) vesicular bodies. *Journal of Experimental Botany* **59**, 1935–1950.

Quint M, Barkawi LS, Fan K-T, Cohen JD, Gray WM. 2009. Arabidopsis IAR4 modulates auxin response by regulating auxin homeostasis. *Plant Physiology* **150**, 748–758.

**R Core Team**. 2014. *R: A language and environment for statistical computing*. Vienna, Austria: R Foundation for Statistical Computing.

Ren Z, Zheng Z, Chinnusamy V, Zhu J, Cui X, Iida K, Zhu J-K. 2010. RAS1, a quantitative trait locus for salt tolerance and ABA sensitivity in arabidopsis. *Proceedings of the National Academy of Sciences* **107**, 5669–5674.

Ribas-Carbo M, Aroca R, González-Meler MA, Irigoyen JJ, Sánchez-Díaz M. 2000. The electron partitioning between the cytochrome and alternative respiratory pathways during chilling recovery in two cultivars of maize differing in chilling sensitivity. *Plant Physiology* **122**, 199–204.

Rizzello A, Romano A, Kottra G, Acierno R, Storelli C, Verri T, Daniel H, Maffia M. 2013. Protein cold adaptation strategy via a unique seven-amino acid domain in the icefish (chionodraco hamatus) PEPFS71501 transporter. *Proceedings of the National Academy of Sciences* **110**, 7068–7073.
Rosenzweig C, Parry ML, others. 1994. Potential impact of climate change on world food supply. *Nature* **367**, 133–138.

Sadras VO, Denison RF. 2016. Neither crop genetics nor crop management can be optimised. *Field Crops Research*.

Sadras VO, Reynolds MP, De la Vega A, Petrie PR, Robinson R. 2009. Phenotypic plasticity of yield and phenology in wheat, sunflower and grapevine. *Field Crops Research* **110**, 242–250.

Sambatti J, Caylor KK. 2007. When is breeding for drought tolerance optimal if drought is random? *New Phytologist* **175**, 70–80.

Segura V, Vilhjálmsson BJ, Platt A, Korte A, Seren Ü, Long Q, Nordborg M. 2012. An efficient multi-locus mixed-model approach for genome-wide association studies in structured populations. *Nature Genetics* **44**, 825–830.

Self SG, Liang K-Y. 1987. Asymptotic properties of maximum likelihood estimators and likelihood ratio tests under nonstandard conditions. *Journal of the American Statistical Association* **82**, 605–610.

Serieys H. 1984. Wild helianthus species, a potential source of androsterilities. *Second eucarpia meeting on the sunflower*. Leningrad.10–14.

Sosnowski O, Charcosset A, Joets J. 2012. BioMercator v3: An upgrade of genetic map compilation and quantitative trait loci meta-analysis algorithms. *Bioinformatics* **28**, 2082–2083.

Talanova VV, Titov AF, Topchieva LV, Frolova SA. 2012. Effects of abscisic acid treatment on the expression of cysteine proteinase gene and enzyme inhibitor during wheat cold adaptation. *Russian Journal of Plant Physiology* **59**, 581–585.

Van Eeuwijk F, Bink M, Chenu K, Chapman S. 2010. Detection and use of qTL for complex traits in multiple environments. *Current Opinion in Plant Biology* **13**, 193–205.

Via S, Lande R. 1985. Genotype-environment interaction and the evolution of phenotypic plasticity. *Evolution*, 505–522.

Vile D, Pervent M, Belluau M, Vasseur F, Bresson J, Muller B, Granier C, Simonneau T. 2012. Arabidopsis growth under prolonged high temperature and water deficit: Independent or interactive effects? *Plant, Cell & Environment*.

Villalobos F, Hall A, Ritchie J, Orgaz F. 1996. OILCROP-SUN: A development, growth and yield model of the sunflower crop. *Agronomy Journal* **88**, 403–415.

Villanova L, Smith-Miles K, Hyndman. RJ. 2011. EMMA: Evolutionary Model-based Multiresponse Approach.

Zhang X, Liu S, Takano T. 2008. Two cysteine proteinase inhibitors from arabidopsis thaliana, AtCYSa and AtCYSb, increasing the salt, drought, oxidation and cold tolerance. *Plant Molecular Biology* **68**, 131–143.
<table><caption>**Table 1** Details on location, treatment, year, testers and observed hybrids for the 17 environments. <sup>a</sup> The locations are designated as follows: AI: Aigrefeuille (Center West), CA: Castelnaudary (South West), CO: Cornebarrieu (South West), GA: Gaillac (South West), VE: Verdun (South West), LO: Loudun (Center West), SE: Segoufielle (South West), CHA: Chateauroux (Center). <sup>b</sup> I for irrigation, NI for non-irrigation.</caption><thead><tr><th>Environment</th><th>Location<sup>a</sup></th><th>Treatment<sup>b</sup></th><th>Year</th><th>Tester for B-line</th><th>Tester for R-line</th><th>Number of hybrids</th></tr></thead><tbody><tr><td>AI08_I</td><td>AI</td><td>I</td><td>2008</td><td>83HR4gms</td><td>FS71501</td><td>192</td></tr><tr><td>AI08_NI</td><td>AI</td><td>NI</td><td>2008</td><td>83HR4gms</td><td>FS71501</td><td>193</td></tr><tr><td>CO09_I</td><td>CO</td><td>I</td><td>2009</td><td>83HR4gms</td><td>FS71501</td><td>278</td></tr><tr><td>CO09_NI</td><td>CO</td><td>NI</td><td>2009</td><td>83HR4gms</td><td>FS71501</td><td>278</td></tr><tr><td>GA09_I</td><td>GA</td><td>I</td><td>2009</td><td>83HR4gms</td><td>FS71501</td><td>275</td></tr><tr><td>GA09_NI</td><td>GA</td><td>NI</td><td>2009</td><td>83HR4gms</td><td>FS71501</td><td>274</td></tr><tr><td>LO10_NI</td><td>LO</td><td>NI</td><td>2010</td><td>83HR4gms</td><td>FS71501</td><td>284</td></tr><tr><td>VE10_I</td><td>VE</td><td>I</td><td>2010</td><td>83HR4gms</td><td>FS71501</td><td>289</td></tr><tr><td>AI09_I</td><td>AI</td><td>I</td><td>2009</td><td>SOLR001M</td><td>AT0521</td><td>280</td></tr><tr><td>AI09_NI</td><td>AI</td><td>NI</td><td>2009</td><td>SOLR001M</td><td>AT0521</td><td>280</td></tr><tr><td>VE09_I</td><td>VE</td><td>I</td><td>2009</td><td>SOLR001M</td><td>AT0521</td><td>273</td></tr><tr><td>VE09_NI</td><td>VE</td><td>NI</td><td>2009</td><td>SOLR001M</td><td>AT0521</td><td>273</td></tr><tr><td>CA10_NI</td><td>C1</td><td>NI</td><td>2010</td><td>SOLR001M</td><td>AT0521</td><td>306</td></tr><tr><td>CO08_I</td><td>CO</td><td>I</td><td>2008</td><td>SOLR001M</td><td>AT0521</td><td>249</td></tr><tr><td>CO08_NI</td><td>CO</td><td>NI</td><td>2008</td><td>SOLR001M</td><td>AT0521</td><td>249</td></tr><tr><td>SE10_NI</td><td>SE</td><td>NI</td><td>2010</td><td>SOLR001M</td><td>AT0521</td><td>285</td></tr><tr><td>CHA10_I</td><td>CHA</td><td>NI</td><td>2010</td><td>SOLR001M</td><td>AT0521</td><td>306</td></tr></tbody></table>
<table><caption>**Table 2 Description of abiotic stress indicators simulated by the crop model.** The SUNFLO crop model was used to simulate the interactions between plant growth and available environmental resources. The evolution of resource level and abiotic constraints during the crop cycle was summarized via computing eight stress indicators during seven cropping periods: vegetative, flowering, grain filling, and their combination. $1_{[x]}$ equals 1 if x is true and 0 else.</caption><thead><tr><th>Stress</th><th>Symbol</th><th>Description</th><th>Unit</th><th>Formula</th></tr></thead><tbody><tr><td>temperature</td><td>HTi</td><td>high temperature (continuous)</td><td>-</td><td>∫ 1-<i>HTRUE</i> dt</td></tr><tr><td>temperature</td><td>HTs</td><td>high temperature (discrete)</td><td>d</td><td>Σ 1<sub>[T<sub>m</sub> > 28]</sub></td></tr><tr><td>temperature</td><td>LTi</td><td>low temperature (continuous)</td><td>-</td><td>∫ 1-<i>LTRUE</i> dt</td></tr><tr><td>temperature</td><td>LTs</td><td>low temperature (discrete)</td><td>d</td><td>Σ 1<sub>[T<sub>m</sub> < 20]</sub></td></tr><tr><td>water</td><td>FTSW</td><td>Edaphic water deficit (continuous)</td><td>-</td><td>∫ 1-<i>FTSW</i> dt</td></tr><tr><td>water</td><td>ETR</td><td>Edaphic water deficit (discrete)</td><td>d</td><td>Σ 1<sub>[ET/PET < 0.6]</sub></td></tr><tr><td>nitrogen</td><td>NAB</td><td>Absorbed nitrogen</td><td>kg ha<sup>-1</sup></td><td>∫ NAB dt</td></tr><tr><td>nitrogen</td><td>NNI</td><td>Nitrogen deficit (continuous)</td><td>-</td><td>∫ 1-<i>NNI</i> dt</td></tr></tbody></table>
<table><caption>**Table 3** Results of the regression between grain yield and abiotic environmental indicators: cold, nitrogen, and water stresses. The plasticity of the three control genotypes (i.e., commercial varieties) are presented as their slopes and p-values of the Fisher test in the best-fit model, with the proportion of explained variance for each control genotype. <sup>a</sup> LTi_veg: integration of low temperatures during the vegetative stage, <sup>b</sup> NAB: absorbed nitrogen during the whole growth period, <sup>c</sup> FTSW_veg+flo: fraction of the transpirable soil water during the vegetative and the flowering periods. The best explanatory linear model was determined with AIC and was selected from all regression models with one to four linear regressors and one indicator per abiotic stress.</caption><thead><tr><th></th><th></th><th>Unit</th><th>Melody</th><th>Pacific</th><th>Pegasol</th></tr></thead><tbody><tr><td>Yield</td><td></td><td>q ha<sup>-1</sup></td><td>32.41</td><td>29.18</td><td>31.30</td></tr><tr><td>Cold stress<sup>a</sup></td><td>slope</td><td>q ha<sup>-1</sup> d<sup>-1</sup></td><td>-0.36</td><td>-0.36</td><td>-0.32</td></tr><tr><td>Cold stress</td><td>p-value</td><td></td><td>2.39 10<sup>-4</sup></td><td>1.97 10<sup>-3</sup></td><td>1.95 10<sup>-3</sup></td></tr><tr><td>Nitrogen stress<sup>b</sup></td><td>slope</td><td>q ha<sup>-1</sup> kg<sup>-1</sup></td><td>0.27</td><td>0.22</td><td>0.25</td></tr><tr><td>Nitrogen stress</td><td>p-value</td><td></td><td>1.16 10<sup>-3</sup></td><td>7.66 10<sup>-3</sup></td><td>6.6 10<sup>-4</sup></td></tr><tr><td>Water stress<sup>c</sup></td><td>slope</td><td>q ha<sup>-1</sup> d<sup>-1</sup></td><td>-0.16</td><td>-0.12</td><td>-0.13</td></tr><tr><td>Water stress</td><td>p-value</td><td></td><td>3.27 10<sup>-3</sup></td><td>2.91 10<sup>-2</sup></td><td>8.21 10<sup>-3</sup></td></tr><tr><td>Explained variance</td><td></td><td></td><td>0.93</td><td>0.94</td><td>0.90</td></tr></tbody></table>
Accepted Article

<table><caption>**Table 4** Variation and correlation of oil yield plasticity for different abiotic stresses in the diversity panel. Minimum, mean, maximum, and variance of the plasticity phenotypes and their correlations. The plasticity phenotypes are calculated as the slopes in a linear model including the three stress indicators that best characterized the environments. <sup>a</sup> cold stress plasticity (q h<sup>-1</sup> d<sup>-1</sup>), <sup>b</sup> nitrogen stress plasticity (q h<sup>-1</sup> kg<sup>-1</sup>), <sup>c</sup> water stress plasticity (q h<sup>-1</sup> d<sup>-1</sup>).</caption><thead><tr><th rowspan="2"></th><th rowspan="2">Min.</th><th rowspan="2">Mean</th><th rowspan="2">Max.</th><th rowspan="2">Var.</th><th colspan="2">Cor.</th></tr><tr><th>Nitrogen</th><th>Water</th></tr></thead><tbody><tr><td>Cold stress<sup>a</sup></td><td>-0.37</td><td>-0.23</td><td>-0.09</td><td>2.5 10<sup>-3</sup></td><td>0.12</td><td>0.03</td></tr><tr><td>Nitrogen stress<sup>b</sup></td><td>0.01</td><td>0.11</td><td>0.21</td><td>8.1 10<sup>-4</sup></td><td></td><td>0.62</td></tr><tr><td>Water stress<sup>c</sup></td><td>-0.25</td><td>-0.08</td><td>0.16</td><td>3.7 10<sup>-3</sup></td><td></td><td></td></tr></tbody></table>
<table><caption>**Table 5 Positions and estimated effects of SNP and genes associated with oil yield plasticity to cold stress.** Genetic distance determined on a consensus map built from three RIL populations (see details in Materials and Methods). Physical position was determined via BLASTing against the genome of line XRQ (early access to version HanXRQv1.1). Association tests were performed using the MLMM procedure unless otherwise stated and either with or without the maintainer/restorer status structure (noted as B/R status). Effects of SNP were estimated using a linear model including all associated SNP. Their relative effect compared with the average plasticity of the diversity panel is indicated. Sunflower gene names correspond to the genome HannXRQ v1.1 (early access). The closest *Arabidopsis thaliana* homologue is indicated with its TAIR code, name, and a brief functional description.</caption><thead><tr><th>SNP name</th><th>QTL</th><th>LG</th><th>Genetic distance (cM)</th><th>Physical position (bp)</th><th>Association test p-value</th><th>Association model</th><th>Effect on plasticity</th><th>Effect percentage</th></tr></thead><tbody><tr><td>AX-84511295</td><td>LG05.17</td><td>5</td><td>16.6</td><td>14,284,297</td><td>1.54 10<sup>-8</sup></td><td>Without B/R status</td><td>2.80 10<sup>-2</sup></td><td>12%</td></tr><tr><td rowspan="3">AX-84248033</td><td rowspan="3">LG05.64</td><td rowspan="3">5</td><td rowspan="3">64.4</td><td rowspan="3">206,893,338</td><td>5.70 10<sup>-6</sup></td><td>EMMA with B/R status</td><td>3.15 10<sup>-2</sup></td><td>13%</td></tr><tr><td>3.64 10<sup>-6</sup></td><td>MLMM with B/R status</td><td>3.34 10<sup>-2</sup></td><td>14%</td></tr><tr><td>6.76 10<sup>-7</sup></td><td>With B/R status</td><td>-2.35 10<sup>-2</sup></td><td>-10%</td></tr><tr><td>AX-84488969</td><td>LG09.27</td><td>9</td><td>27.4</td><td>203,648,151</td><td>3.87 10<sup>-6</sup></td><td>Without B/R status</td><td>-5.06 10<sup>-2</sup></td><td>-21%</td></tr><tr><td>AX-84358846</td><td>LG10.34</td><td>10</td><td>34.1</td><td>220,548,129</td><td>5.11 10<sup>-7</sup></td><td>With B/R status</td><td>-3.14 10<sup>-2</sup></td><td>-13%</td></tr><tr><td>AX-84316040</td><td>LG10.45</td><td>10</td><td>44.8</td><td>93,526,595</td><td></td><td></td><td></td><td></td></tr><tr><td></td><td></td><td></td><td></td><td>or 217,685,071</td><td></td><td></td><td></td><td></td></tr><tr><td>AX-84586324</td><td>LG13.72</td><td>13</td><td>72.4</td><td>190,712,390</td><td>1.17 10<sup>-7</sup></td><td>Without B/R status</td><td>-2.94 10<sup>-2</sup></td><td>-12%</td></tr><tr><td>AX-84436021</td><td>LG14.27</td><td>14</td><td>27.0</td><td>168,814,054</td><td>2.70 10<sup>-7</sup></td><td>With B/R status</td><td>3.61 10<sup>-2</sup></td><td>15%</td></tr><tr><td>AX-84337313</td><td>LG16.48</td><td>16</td><td>47.7</td><td>129,668,002</td><td>6.65 10<sup>-8</sup></td><td>Without B/R status</td><td>-2.72 10<sup>-2</sup></td><td>-12%</td></tr><tr><td>AX-84507515</td><td>LG17.49</td><td>17</td><td>49.1</td><td>153,030,442</td><td>6.81 10<sup>-6</sup></td><td>EMMA without B/R status</td><td>3.55 10<sup>-2</sup></td><td>15%</td></tr></tbody></table>

(Table 5 continued...)

<table><thead><tr><th>SNP name</th><th>Sunflower gene</th><th>TAIR code</th><th>Description</th><th>Arabidopsis name</th></tr></thead><tbody><tr><td>AX-84511295</td><td>HannXRQ_Chr05g0131961</td><td>AT1G68570</td><td>Oligopeptide transporter</td><td>NPF3.1</td></tr><tr><td>AX-84248033</td><td>HannXRQ_Chr05g0159711</td><td>AT4G26720</td><td>Protein phosphatase</td><td>PPX1</td></tr><tr><td>AX-84488969</td><td>HannXRQ_Chr09g0274811</td><td>AT4G33495</td><td>Root development, Temperature sensitive</td><td>RPD1</td></tr><tr><td>AX-84358846</td><td>HannXRQ_Chr10g0311881</td><td>AT3G14150</td><td>ROS protection</td><td>HAOX2</td></tr><tr><td>AX-84316040</td><td>HannXRQ_Chr10g0292461</td><td>AT5G14950</td><td>N-glycan mannose hyper-osmotic salinity stress</td><td>GMII</td></tr><tr><td>AX-84586324</td><td>HannXRQ_Chr10g0311021</td><td>AT3G12490</td><td>Cystatin, protease inhibitor, Tolerance to abiotic stresses incl. cold</td><td>CYS6</td></tr><tr><td>AX-84436021</td><td>HannXRQ_Chr14g0460351</td><td>AT5G46040</td><td>Oligopeptide transporter</td><td>NPF5.3</td></tr><tr><td>AX-84337313</td><td>HannXRQ_Chr16g0519611</td><td>AT3G07450</td><td>Lipid Transfer Protein</td><td></td></tr><tr><td>AX-84507515</td><td>HannXRQ_Chr04g0119921</td><td>AT1G24180</td><td>Root development, Auxin</td><td>IAR4</td></tr></tbody></table>
Accepted Article

<figure><img src="image_4.png" /><figcaption>**Figure 1** Relation between observed grain yield and several abiotic stress indicators.</figcaption></figure>

Panels A, B, and C display the regression line between grain yield and water stress indicators, computed at different levels: using climatic data only (panel A, precipitation - potential evapotranspiration, Pearson correlation (r) of -0.01), using both climatic and crop data (panel B, precipitation - potential evapotranspiration + irrigation + soil water capacity, r = -0.65, p-value = 9.6 10⁻⁶) and using simulated plant data (panel C, evapotranspiration ratio, r = -0.86, p-value = 3.3 10⁻¹²). Panel D displays grain yield predicted as a function of combined abiotic stresses (linear model of water, nitrogen, and cold stress indicators, r = 0.91, root mean square error of 0.21 t ha⁻¹).
<figure><img src="image_5.png" /></figure>

**Figure 2** Heat map of the correlations between stress indicators. The three selected stress indicators are indicated in blue, yellow and green for water, cold and nitrogen respectively. For water stress, the fraction of transpirable soil water (FTSW) represent yield limitation through water deficit (integration of 1 minus FTSW); ETR is the conditional sum of days, if the ratio of the real evapotranspiration (ET) to potential evapotranspiration (PET) was less than 0.6 (threshold for photosynthesis limitation). For cold stress, LTs is the conditional sum of days if mean air temperature was below 20°C and LTi represent low temperatures impact on photosynthesis (integration of 1 minus equation 2). Heat stress indicators were computed following the same logic, albeit representing high temperatures impact on photosynthesis. Equations (2) and (3) are used in the crop model to define the radiation use efficiency (RUE) response to temperature (Villalobos et al., 1996). For nitrogen deficit, NAB is the amount of absorbed nitrogen in the considered cropping period and NNI is the sum of 1 minus nitrogen nutrition index, which indicates crop nitrogen deficit (Lemaire and Meynard, 1997; Debaeke et al., 2012).

<figure><img src="image_6.png" /></figure>
<figure><img src="image_7.png" /><figcaption>**Figure 3** Range of variation in abiotic stresses in the multi-environment trial. Water (A), cold (B), and nitrogen stress (C) indicators, computed by the SUNFLO crop model for each environment, and averaged over the three control genotypes. The SUNFLO crop model was used to compute stress indicators for three control commercial hybrids (Melody, Pacific, Pegasol) according to the developmental stage of the crop: vegetative growth (green), flowering (yellow), and seed filling (brown).</figcaption></figure>
<figure><img src="image_8.png" /></figure>

<figure><img src="image_9.png" /><figcaption>**Figure 4** Regression of oil yield against the cold stress indicator for the most tolerant (green triangles) and sensitive (red squares) panel lines showing the ability of multi-stress modelling to better characterize the environment. A) Regression in a single (cold) stress indicator model. B) Regression in a multi-stress (cold, nitrogen, water) indicator model.</figcaption></figure>

<figure><img src="image_10.png" /></figure>
Accepted Article

<figure><img src="image_11.png" /><figcaption>**Figure 5** Relation between the multi-stress plasticity index and the average oil yield showing differences in stress tolerance strategies in the diversity panel of sunflower lines. The three-branch star represents the strategy of the most extreme panel lines in response to combined abiotic stresses. The length of star branches represents the relative plasticity against the corresponding stress: longer equals more sensitive.</figcaption></figure>