Climate change has increased genotype-environment interactions in wheat breeding


 The International Maize and Wheat Improvement Center (CIMMYT) develops and distributes annually elite wheat lines as international trials worldwide to assess their performance in different environments and utilization by partners for use in breeding or release as varieties. However, as elsewhere, the collaborator test sites where trials are evaluated have experienced climate change, with implications for how adapted wheat genotypes are bred. Using a standard quantitative genetic model and archived datasets for four global spring wheat trials, we show that the genotype-environment-interaction (GEI) has increased by up to 500% over recent decades. Notably crossover has increased over time, a critical indicator of changes in the ranking of cultivar performance in different environments. Climatic factors explain over 70% of the year-to-year variability in GEI and crossover interactions for yield. Examining yield responses of all genotypes in all trial environments from 1985 to 2017 reveals that climate change has increased GEI by ~ 49% and ranking change by ~38%. Genetic improvement of wheat targeted to high-yielding environments has exacerbated this increase, but the performance of new wheat germplasm developed to withstand heat and drought stress is more adapted and stable, offsetting the increase in ranking changes due to the warmer climate.

Wheat yield must be improved more rapidly to ensure food security for a growing human population amid production constraints due to diminishing resources, land loss, and climate change 1 . Different methods forecast global wheat yields will decline in response to worldwide warming if crop adaptation is neglected, particularly in warmer regions [2][3][4] . This is despite continuous genetic gains, improved agronomic practices 5,6 , and supportive policies 7 . Significantly, the impact of climate warming on crop breeding itself, vital to future agriculture, has yet to be scrutinized.
Genotype-environmental interactions (GEI) are the differential responses of genotypes to environments 8 , and they affect every aspect of decision making in crop breeding. For example, a given superior cultivar may perform poorly in an environment to which it is not adapted, which creates uncertainty for breeders when translating the relative performance of numerous cultivars into long-term strategies 9 . GEI may impact the overall genetic gains and efficiency of many breeding programs by influencing how well target growing environments are identified, and how resources are allocated to different breeding projects. A clear understanding of the temporal and spatial pattern of GEI is indeed prerequisite in planning an effective breeding program 10 . It is particularly useful at the pre-breeding stage to identify useful traits and haplotypes expressing yield stability and adaptation to boost yield potential and broaden the genetic bases of the crop 11 .
However, obtaining this level of information is challenging due to the limitations of current breeding trial design, relatively small area coverage, discontinuous testing periods, and a variety of confounding factors.
Wheat breeders collaborating in the CIMMYT International Wheat Improvement Network (IWIN) test approximately 1,000 new, well-adapted, disease-resistant wheat genotypes each year at approximately 700 field sites in over 90 countries 12,13 . IWIN serves primarily to distribute improved germplasm globally, but the data returned from these trials provide a valuable resource to assess GEI and its long-term trends 14 . To date, IWIN has collected over 10 million raw phenotypic data points, while delivering germplasm that is estimated to be worth several billion dollars annually in extra productivity to hundreds of millions of farmers in less developed countries 15 . Using historical IWIN datasets, here we investigate the changes in GEI over time for four different spring wheat trials ( Fig. 1  evaluates spring bread wheat that is specifically bred to maintain yields under dry and droughtprone conditions 19 . To investigate whether GEI has changed over time, we first used the standard quantitative genetic model (Equation 1) to measure the GEI for each year and nursery. Changes in GEI over time were measured by the sum of squares (SS) of GEI, which increased significantly (p<0.05) and substantially (73%~485%) for all the trials (Fig. 2a). HTWYT showed the most rapid increase in GEI, at a ratio of 14.9% per year, corresponding to a four fold increase in total from 1992 to 2017. The slowest increase in GEI was for IDYN, which nevertheless increased at a rate of 2.1% per year and 72.8% in total from 1980 to 2014. Many factors may have contributed to these increases, including improved germplasm and possibly increased genetic diversity, more frequent extreme weather, and environmental heterogeneity across sites. For example, the number of testing lines increased from 30 to 50 over the years, and the number of test sites increased more than six fold for some nurseries (Supplementary Table S1). Temperatures in the three stages of the wheatgrowing season (vegetative, reproductive, and grain filling) increased significantly over the years, accompanied by erratic precipitation patterns and changes in other climatic variables, such as vapor pressure deficit ( Supplementary Fig. S1). There were also changes in the spatial variation of some climatic variables, such as minimum temperature, humidity, which reflects the enhanced spatial heterogeneity of the wheat-growing environments ( Supplementary Fig. S2).
Different types of GEI can be grouped into two categories -crossover interaction that represents the difference in genotype rankings, and non-crossover interaction that only covers the difference in the level of expression of genetic differentiation between environments 21 . Crossover GEI is what concerns breeders most, as it potentially jeopardizes the selection of the better genotypes, with the likelihood of smaller genetic gains the more cross-overs occur 22 . We can count the number of crossovers for any possible pair of lines using the simple yet powerful Finlay-Wilkinson (FW) regression as expressed in equation (2) 23 . Fig. 2b shows the relative number of crossovers that have occurred for the four trials each year, taking into account that different numbers of lines were tested for some years. All nurseries experienced steady and fast increases in crossover occurrence, suggesting that ranking changes have become more frequent since the four trials began. These increases are markedly higher with the Elite trials of ESWYT and IDYN both of which evaluated high-yielding wheat. In contrast, for HTWYT and SAWYT, the increase in the crossover interaction is relatively small, possibly because the breeding goals of these trials are specifically for heat and drought tolerance. There is an apparent jump in the crossover number around 1994/1995, which corresponds with more testing lines being used during that period (Supplementary Table S1). Genetic variation in the yield-climate response is likely responsible for the GEI and ranking changes. To test this hypothesis, we estimated the relationship between yield variation and climate for each line by fitting the yield variation (yijk -(k + Gik + Ejk)) with six climatic variables split across three stages of the growing season and their quadratic terms. Replacing the GEI component of equation (1) by the yield variation-climate relationship gives a new yield prediction equation (3) in which climate is the only driver of GEI. Results indicate that the regression model between yield variation and climate can well predict the temporal variability in both GEI and number of crossovers. However, the prediction ability of the regression model on GEI/crossover depends on the number of climatic variables included, with R 2 between real and simulated values ranging from 0.48 to 0.87 for GEI, and from 0.63 to 0.97 for the number of crossovers (Supplementary Table   S2). We chose a model by compromising between modeled bias on GEI/crossover and predicted accuracy on the final yield by equation (3)

(Supplementary Methods). A comparison between real
and simulated yields suggests that replacing the GEI components by the yield-climate relationships in equation (3) can well predict the observed yields, with an R 2 of 0.84 averaged between training and testing datasets ( Supplementary Fig. S3). Inferred yield variations from the yield-climate relationship explain 71.4% of temporal variability in GEI and 84.5% of variability in crossovers ( Fig. 3). Furthermore, a pairwise comparison between real and predicted values indicates inferred yield variations capture 56.8% of the real GEI averaged across years, with the highest value for HTWYT (69.1%) and the lowest for ESWYT (42.0%). For crossovers, modeled yield variations correctly predict the occurrence of ranking changes in over 70% of line pairs, ranging from 73.7% for ESWYT to 83.8% for HTWYT (Supplementary Table S2). Theoretically, yields of any line at any site can be predicted from the genotype-climate relationship and the site's climate, so we applied this to all the line-location combinations to eliminate confounding effects in estimating the GEI. By doing so, we were able to untangle the contributions of climate on GEI by using the mean of GEI/crossover across any subset of genotypes with respect to years. Similarly, by splitting the lines into groups going from old to new genotypes, we could separate the effects of breeding progress by using the mean of GEI/crossover across different climatic environments with respect to line groups. For each trial, there is a maximum of 27 line groups with 50 lines each representing breeding progress from 1985 to 2017.
To observe the trends, the modeled values of the SS of both GEI and number of crossovers were plotted against year and line group respectively for each trial (Fig. 4), with the changes in GEI and crossovers converted into percentage relative to the mean of the first 5-years or oldest line group as appropriate.
We found a detectable effect of climate change on both GEI and crossover, with differences apparent among trials. Most trials exhibit increasing trends for both GEI and crossovers, suggesting climate change generally reduces the performance stability of genotypes across environments. Three trials -ESWYT, IDYN, and SAWYT -show increasing GEI and crossover at mean rates of 1.01% and 0.34%, respectively, corresponding to increases of 24% and 11% in total from 1985 to 2017 (Fig. 4a, 4b). Among these three nurseries, only IDYN demonstrates significant (P<0.05) changes for both GEI and crossover, and indeed shows the largest increase in genotypic ranking changes (38%). HTWYT is the only nursery experiencing decreasing trends for both GEI (23%) and crossover (27%), probably because it includes many heat-tolerant genotypes.
Not all line groups respond similarly to climate change, even in the same nursery. These differences indicate how shifts in relative performance of some varieties may lead to deviation from the planned course towards breeding goals. On a positive note, the information provides breeders with options for adapting improved germplasm for warming climates by taking advantage of the varietal differences revealed.
Climate change affects GEI/crossover through three mechanisms. First, the warming trends across the test sites accelerate wheat growth, reducing mean wheat yield across the environments 2 .
The spatial analysis of crossover for reported yields reveals a normal distribution of crossover possibility across environments, with the highest chance of crossover occurring in environments with a yield around 2t/ha (Supplementary Fig. S4). Current yield averaged across all test environments is 4.1t/ha, with warmer climate tend to give lower mean yields in most environments 3 , potentially increasing the occurrence of crossover. Secondly, the warmer climate likely increases the variations of several environmental factors across sites, such as vapor pressure deficit, directly resulting in larger yield variation and an increase in GEI. Thirdly, the increased frequency of extreme events under climate change results in unexpected yield penalties for some sites and years, driving up the GEI and crossover. For instance, for most trials there have been more years with extreme SS of GEI since 2000 than before ( Supplementary Fig. S5), which partly relates to the spatial variation in precipitation. When taking equal account of extreme points, estimated increases in GEI/crossover due to climate change are higher than current estimates for all the trials (Supplementary Fig. S6), suggesting that extreme weather has had a considerable effect on GEI. Genotypes tested in each trial are selected from CIMMYT international breeding programs that contribute to over 80% varieties in Asia, Africa, and Latin America 24 , so we can assume they are representative of progress in wheat breeding. Compared to the impacts caused by climate change, the effects of genetic improvement since the 1980s have various effects on GEI or crossover depending on the trial. Averaging across the four trials, the SS of GEI decreases at a ratio of 0.48%, and the number of crossovers declines by 0.02% (Fig 4c, 4d). Targeted breeding goals, including selection in simulated environments, different distribution of test sites, and distinct testing periods could be responsible for the varying influence of genetic improvement among trials.
Of the four trials, the two processing stress-tolerant selection (HYWTY and SAWAY) have exhibited substantial decreases in both GEI and crossover, with mean ratios of -3.1% and -0.7%, respectively. These large decreases in GEI indicate an improved and stable genotypic adaptation of wheat cultivars for a wider range of environment and can be largely attributed to recent global scale efforts to develop drought and heat tolerant wheat in coping with climate change 25,26 .
Considering the trials that include high-yielding wheat in optimal management and environments, ESWYT has experienced a decrease in GEI (-2.3%), but an increase in altered genotype rankings (0.6%), while IDYN is the only trial showing increasing trends for both GEI and crossover due to genetic improvement, with ratios of 6.5% and 0.8%, respectively. The higher values observed for IDYN may be because durum wheat is cultivated in relatively restricted regions, limiting the development of genotypes suited to other types of environment. It was also observed that breeding effects vary across climatic environments, particularly for the change in the number of crossovers. This is because of the large climatic variation in some years, implying that the efficiency of breeding efforts is dependent on the testing environments in which the crop is grown.
Statistical, the estimated trends are not significant (P<0.1) in many cases (Fig. 4). This is because of the large variations in the GEI/crossover predictions caused by increasingly extreme weather and also the yield variation-climate model chosen for this analysis. Model sensitivity analysis indicates that estimated trends vary when choosing different yield variation-climate models. A more complex model including detailed climate factors leads to higher prediction accuracy, with the drawback that large variations in predicted GEI/crossovers obscure the temporal trends of effects. Alternative models using a varying number of climate factors were found to have different explanatory powers but gave similar results as described in the Supplementary information ( Supplementary Fig. S7). Another consideration is that observation errors and missing data for both yield and phenologies are inevitable in such a large dataset that relies on records from breeders worldwide who operate with different resources. For example, a large portion of data points lack recorded phenologies, particularly for heading. In many cases, we used a surrogate of location mean or a phenology prediction method to estimate the dates for computing the climatic variables, potentially increasing the uncertainty of the results (Supplementary Method).
Despite the limitations, our results are sufficiently meaningful and robust to reach two main conclusions. First, climate change since the 1980s has increased GEI in global wheat breeding.
This finding confirms a widely held belief that a higher GEI would be expected due to projected warmer climate. Importantly, our study quantifies the effects on the global scale. The increase in GEI, and particularly altered genotype rankings, has a number of implications for breeding, including less efficient genotype selection 27 , the discarding of more potentially useful lines 28 , increased environmental heterogeneity, modification of mega-environment classification 26 , and even reduced genetic gains. We find climate change from 1985 to 2017 has increased GEI by up to 49% depending on the breeding program, highlighting the increasingly difficult task for wheat breeders to accelerate the genetic gains under an erratic climate.
Secondly, in spite of climate warming, there have been increases in grain yield over the years due to improved germplasm. However, different global breeding trials differ in how they have responded with respect to GEI, either accelerating or offsetting observed ranking changes caused by climate change. The different crossover trends partly explain observed distinctions in ranking change (Fig. 2b) and the real genetic gains of wheat yield (Supplementary Fig. S8) between the high-yielding and stress breeding programs. Genetic development in the high-yielding breeding programs has jointly contributed to the increase in ranking change, potentially undermining the capacity of high-yielding wheat to adapt to climate change. Conversely, cultivars in abiotic stress breeding programs, particularly the heat-tolerant wheat, have substantially improved the yield performance stability across broad environments, consequently increasing the capacity of adaptation to the warmer climate as expected. These results indicate that breeding efforts to produce heat-tolerant wheat can not only benefit current climate change adaptation, but also offer an opportunity to obtain larger and faster genetic gains in the future. For the datasets, we extracted the yield and phenology (sowing, heading, and maturity dates) for each line, environment, and year. Datapoints with at least two phenologies (sowing and maturity) and yield were kept.
Any replicates among the datapoints were averaged for the line-location-year combination. Before analysis, we conducted a data quality control to remove likely observation errors by discarding lines deemed to be phenological outliers for each location. The lines whose phenology was beyond 1.5 times the interquartile were thus treated as outliers and removed.
Yield data were paired with reconstructed daily weather data (Supplementary Method). The daily weather dataset covers the period from 1985 to 2018, and includes maximum and minimum temperatures, shortwave radiation, precipitation, vapor pressure deficit, relative humidity, and wind speed. The dataset is based on the ERA5 gridded dataset, a global reanalysis dataset at 30 km resolution, with consistent data sequences from 1980 onwards in hourly to 3-hourly intervals. ERA5 was interpolated for each nursery site then quality-controlled against the nearest observation by Meteoblue, a specialist company (Switzerland).
Reported phenology was used to define the growing period, and compute the mean or accumulated climatic variables for the three growth stages as relevant. We used the method described in ref 14 to define the three growth stages of wheat: vegetative is from sowing to 300 growing degree-days (GDD) before heading; reproductive is from 300 GDD before heading to 100 GDD after heading; and grain-filling is from 100 GDD after heading to harvest.
The standard quantitative genetic model. We used a standard quantitative genetic model to investigate the GEI 29 .
where yijk is the yield response for genotype i in environment j and year k, is the grand mean, G and E are the genotypic and environmental effects, respectively, and GEijk is the interaction. Observed GEI was directly computed by two-way ANOVA. We did not separate error from GEI and considered the sum of squares of GEijk as the value of GEI. The crossover between any pair of lines was measured by using the Finlay-Wilkinson (FW) model with regression of individual genotype performance against environmental means 23 .
where bik is the genotypic slope on environmental means such that each genotype has an intercept Gik and the slope bik, with a crossover (equivalent to a ranking change) if the regression lines of any two genotypes interact within the environmental ranges. We computed the total number of crossovers divided by the number of genotypes for each trial and year.
The modified quantitative genetic model. Wheat GEI is caused by the phenotypic variation in yield when a genotype is exposed to different environments 30 . In most sites, particularly those for ESWYT and IDYN, agronomic management is typically optimized, so climate is the main factor causing such variation. Therefore, we replaced the GEijk component of equation (1) by a genotypic yield-climate response as described in equation (3) to simulate the yield across genotypes, environments, and years.
× is an extracted yield-climate relationship for the genotype i, which was estimated by fitting the yield variation − ( + + ) and corresponding climate vector Wijk. The climate vector included maximum and minimum temperature, shortwave radiation, vapor pressure deficit, relative humidity, and wind speed, split into three growing stages -vegetation (sowing to heading), reproduction (heading to flowering), and grain filling (flowering to maturity). The quadratic terms for each climatic variable were also included to take into account the non-linear effects. After fitting , the sum of squares (SSk) indicating the GEI at any year k can be directly estimated with equation (4) We also estimated the number of crossovers for modeled yields with the FW model, but instead regressing modeled yield variation plus genotype variation × + against environmental variation . As we intended to investigate the temporal trends in GEI/crossovers and there is a structural difference in germplasm, environment, and management between trials, all the regressions were run separately for each trial and year, e.g. ESWYT in 2001 only, ESWYT in 2002 only and so on.
Untangling the contributions of climate and breeding to temporal changes in GEI. Genotype and environmental effects on GEI can be separated by statistical models such as the additive main effect and multiplicative model (AMMI) 31 . However, these models rely on the trial design and are challenging when investigating the temporal variability and eliminating confounding factors such as site/line changes. We therefore extended the model (3) to include all sites and lines to untangle the contribution of climate trend and breeding progress on the changes in GEI/crossover over time. For this, we developed the model (5)  = + + + + × where is the grand mean of yield for each nursery estimated from the reported yield; Yk is the yearly variation of the yield from the grand mean; Gi and Ej are the genotypic and environmental variations from the grand mean, respectively, which are the deviations from the genotypic or environmental means, estimating from the reported yield for specific years but averaged over years if there were multiple year observations; represents the fitting parameters estimated from equation (3); and Wjk indicates the climatic variables in trial j and year k, computed from the daily weather data based on identical trial phenology. For each test site of the trial, we used the median of sowing, heading, and maturity dates from 1985 to 2017 to define the growth stages. For sites that lack observed phenology, such as heading, the date was estimated from its growing-season mean temperatures by taking into account the known degree-day requirements of the genotypes and sowing date (Supplementary Method).
For each nursery, we pooled all genotypes into line groups according to their testing years, each containing 50 lines. We applied equation (5) to all line groups (the maximum of 27 was for ESWYT) with the climatic variables from 1985 to 2017. For each line group and year, we repeated the calculation of the sum of squares of GEI and the number of crossovers with the above approach. Effect of climate change on GEI/crossover was investigated by a trend analysis between GEI/crossover plotted against years, while for breeding progress, the effect was measured by the changing trend of the GEI/crossover ploted against line groups.