Modeling crop breeding for global food security during climate change

The 2009 Declaration of the World Summit on Food Security requires a sustainable increase in global food production of 70% by 2050. However, the rate of increase in global crop grain yield, which rose during the late 20th century, is now declining (Alexandratos & Bruinsma, 2012), and most studies predict a decrease in grain yield of 10%–25% in the late 21st century as a result of higher global temperatures (IPCC, 2014). Grain crops such as wheat, maize, and rice produce more than 60% of the Received: 6 July 2018 | Revised: 27 September 2018 | Accepted: 2 October 2018 DOI: 10.1002/fes3.157


| INTRODUCTION
The 2009 Declaration of the World Summit on Food Security requires a sustainable increase in global food production of 70% by 2050. However, the rate of increase in global crop grain yield, which rose during the late 20th century, is now declining (Alexandratos & Bruinsma, 2012), and most studies predict a decrease in grain yield of 10%-25% in the late 21st century as a result of higher global temperatures (IPCC, 2014). Grain crops such as wheat, maize, and rice produce more than 60% of the stress tolerance (HST) simultaneously as the climate warms (Battisti & Naylor, 2009).
We stochastically model a rapid-cycle spring wheat (Triticum aestivum L.) breeding program based on best linear unbiased prediction (BLUP) of breeding value with pedigree relationships (known as pedigree BLUP, or ABLUP) and compare truncation selection with optimal contributions selection (OCS) during 60 years of global warming. Complex polygenic variation exists for HST in wheat landraces and contemporary lines (Atlin et al., 2017;Farooq et al., 2011), and is assumed to be present in the founder population. Model conditions for the founder population include genetically diverse contemporary spring wheat lines with moderate to low heritability for grain yield, disease resistance, stem strength, and HST during anthesis. HST is selected to match expected increases in average temperature during anthesis from 2017 to 2077, as predicted by global warming models (IPCC, 2014), while simultaneously selecting for an economic index composed of grain yield, stem strength, flowering time, and disease resistance.
Optimal contributions selection was developed in animal breeding to maximize genetic gain for an economic index under the constraint of the maximum permissible target inbreeding rate (Woolliams, Berg, Dagnachew, & Meuwissen, 2015), and was recently applied to self-pollinating grain crops (Cowling et al., 2017). OCS produced sustainable and superior long-term genetic improvement in an economic index including grain yield compared with truncation selection (Cowling et al., 2017). Here, we compare OCS vs. truncation selection to improve both HST and grain yield simultaneously in a model wheat breeding program, and thereby avoid the threat of reductions in grain yield caused by global warming during the 21st century.
The value of HST for grain yield depends on the prevailing climate, with little benefit from excess tolerance, but probable disaster resulting from insufficient tolerance. We thus use strict constraints to ensure that HST levels match the prevailing requirements predicted by global warming models, while targeting maximum long-term response in the economic index and its components.
Grain yield of spring wheat is expected to decrease by 10% per 1°C increase above the critical temperature threshold of 25°C during anthesis, reaching zero yield at the limiting temperature of 35°C (Deryng, Conway, Ramankutty, Price, & Warren, 2014). This is in line with other studies which show that yield of wheat is reduced as growing season temperature increases (Lobell & Field, 2007), with a failure temperature during anthesis of 34°C (Hatfield et al., 2011) or 35°C (Prasad & Djanaguiraman, 2014). We use this information to model HST in wheat at an effective temperature of 30°C (HST 30 ) ( Figure 1). Effective temperature is the temperature half-way between the mean and the highest temperature each day for four consecutive days during anthesis (Deryng et al., 2014). Our founder wheat population is assumed to have a F I G U R E 1 A model of heat stress tolerance (HST 30 ) during anthesis and yield penalty in wheat, following Deryng et al. (2014). (a) Yield penalty begins when the threshold effective temperature during anthesis exceeds 25°C in 2017 (blue line) or 29°C in 2077 (red line), following successful selection for HST 30 score of +4 between 2017 and 2077. Grain yield decreases by 10% for every 1°C increase in effective temperature above the threshold temperature, such that a 50% yield penalty occurs at 30°C in 2017 or 34°C in 2077 following successful selection for HST 30 score of +4. (b) Normal curve of true breeding values for HST 30 in 2017 (blue curve), and in 2077 (red curve) after successful selection for a target increase in HST 30 score of +4 between 2017 and 2077 mean HST 30 score of 30 in 2017, which must increase by 2077 to match global warming models (IPCC, 2014). In our model, if selection for HST 30 fails to meet the target in each year, grain yields are penalized by 10% per unit below target HST 30 , following the yield loss assumptions of Deryng et al. (2014) (Figure 1, Supporting information Table S1).
Plant breeders of inbred crops often use independent culling, where a minimum performance is set for each trait (Bernardo, 2010). In this process, individuals which survive culling for the first trait are selected for the next trait, and so on, until only a few selected individuals are available for mating to begin the next cycle. Independent culling was shown to be less efficient than index selection for economic traits in animal breeding (Hazel & Lush, 1942) and in self-pollinating crops (Pesek & Baker, 1969).
During the past 30-40 years, index selection based on BLUP breeding values has been adopted widely in animals and perennial tree crops, but not in self-pollinating grain crops (Bauer & Léon, 2008). High accuracy of BLUP breeding values was obtained with 2-year cycles of recurrent selection in a self-pollinating crop (Cowling et al., 2015), and long-term genetic gain in an economic index composed of several low heritability traits was greater with OCS compared to truncation selection (Cowling et al., 2017).
We build on this experience to model wheat breeding from 2017 to 2077 to improve an economic index composed of grain yield and three other economic traits (Cowling et al., 2017), while precisely controlling HST 30 to match prevailing climatic conditions. We have chosen an ABLUP breeding method based on all pedigree relationship information across cycles, including selfing, which results in high accuracy of predicted breeding values (Cowling et al., 2015) at relatively low cost, and is suitable for implementation in developing countries. Three levels of selection of HST 30 are modeled: (a) no selection, (b) selection for an increase in HST 30 of +2 units by 2077, which matches the "average" climate change model predictions for +2°C increase in land temperatures in 2077 (IPCC, 2014), and (iii) selection for an increase of HST 30 of +4 units by 2077, which matches the "high" model prediction of +4°C increase in land temperatures in 2077 (IPCC, 2014). For each HST 30 target, we evaluate two methods of selection (OCS and truncation) and two intensities of selection (high and moderate), for a total of 12 scenarios. Achieving HST 30 targets in each year was a primary constraint, and the economic index was improved in the face of that constraint.

| Simulations
The starting values for mean, phenotypic standard deviation, narrow-sense heritability, genetic correlations among traits, and economic weights for grain yield, resistance to a hypothetical disease, stem strength, and flowering time in the model wheat breeding program were the same as reported previously (Cowling et al., 2017), except the economic weight for flowering time was adjusted from −0.020 to −0.025, and a new trait for heat stress tolerance (HST 30 ) was added (Supporting information Table S2). The founder population of genetically diverse wheat had a mean yield of 1.5 t/ha in 2017, which is close to the projected average yield of rainfed wheat in developing countries in 2017 (1.8 t/ha) (Alexandratos & Bruinsma, 2012).
The recurrent selection scheme was based on 2-year cycles in the population improvement phase of crop breeding (Supporting information Figure S1). Records from S x -derived S x+1 plots were used to predict breeding values of S x individuals, following the theory that the self-family mean provides an improved estimate of the breeding value of the parent for crossing (Walsh & Lynch, 2018). Records were obtained on both cross progeny and selfs of parent plants, because this increased the accuracy of estimated breeding value (EBV) due to inclusion of self-relatives in the analysis (Cowling et al., 2015(Cowling et al., , 2017. When a genotype was selected for crossing, remnant self-progeny seeds were used in crossing (Cowling et al., 2017). We follow the nomenclature in crop breeding where "F" denotes progeny of crosses among homozygous parents, and "S" denotes progeny of crosses among heterozygous parents.
Simulations were carried out using the PopSim module of Genup (http://bkinghor.une.edu.au/genup.htm, accessed 3 March 2018), which was developed and used in the context of the infinitesimal model (Webb et al., 2012). PopSim was modified to include OCS as an option for mate selection and mate allocation decisions at each breeding cycle (Kremer, Newman, Wilson, & Kinghorn, 2010), and to handle bisexuality and selfing as required for self-pollinating crops (Cowling et al., 2017).

| Preparation of founder lines
Cycle 0 was used to prepare the founder populations (Supporting information Figure S1). In the moderate selection intensity scenario, 100 non-inbred individuals were generated by the simulation software on the basis of the starting values (Supporting information Table S2), and the founder population of 100 inbred individuals was formed by single seed descent for 6 generations. Fifty random pair-wise matings were made among these founder inbred lines, and their offspring were intermated again (25 matings) to generate 32 heterozygous segregating S 0 progeny per mating (800 S 0 progeny), plus four selfs per parent plant (200 F 2 progeny), to give a total of 1,000 progeny to begin cycle 1. BLUP-derived EBVs of the 1,000 progeny were used to construct an economic index for use in selection of parents to begin cycle 1 (Cowling et al., 2017).
In the moderate selection intensity scenario, 40 noninbred individuals were generated by the simulation software on the basis of the starting values (Supporting information Table S2), and the founder population of 40 inbred individuals was formed by single seed descent for six generations. Twenty random pair-wise matings were made among these founder inbred lines and their offspring were intermated again (10 matings) to generate 80 heterozygous segregating S 0 progeny per mating (total 800 S 0 progeny), plus 10 selfs per parent plant (200 F 2 progeny), to give a total of 1,000 progeny to begin cycle 1. BLUP-derived EBVs of the 1,000 progeny were used to construct an economic index for use in selection of parents to begin cycle 1 (Cowling et al., 2017).

| Economic index
The emphasis placed on each trait was calculated using the desired gains approach (Brascamp, 1984), implemented using the program Desire (http://bkinghor.une.edu.au/desire.htm, accessed 3 March 2018). Economic weightings were adopted for calculation of BLUP-based economic index values: where e is a vector of implied economic weights, and EBV i,j is a vector of estimated breeding values for individual i calculated by BLUP.
BLUPindex i was calculated from EBVs based on multipletrait BLUP analysis on phenotypes and pedigree information generated in the simulations. Simulation was based on the starting values and genetic parameters for traits (Supporting information Table S2).

| Optimal contributions selection
The selection and mate allocation method used previously (Cowling et al., 2017;Kinghorn, 2011) involves a function whose key components relate to genetic gain and genetic diversity. The optimization of this function results in OCS. OCS was based on the breeding program implementation platform "Matesel" (Kinghorn & Kinghorn, 2018). The practical implementation of this method is based on an evolutionary algorithm, with constraints easily invoked to ensure practical relevance and precise control of response in HST 30 . Matesel dictates which individuals to select and the actual mating allocations and/or selfings to be made.

| Four selection scenarios
From the beginning of cycle 1, we evaluated four scenarios based on two methods of selection (OCS and truncation) and two intensities of selection (high and moderate). The model was run 10 times for each scenario, and the mean and standard deviation reported for each cycle (Supporting information Table S3). The cost of each breeding scenario was approximately the same, since the same total number of genotypes was phenotyped in each scenario in each cycle.

| Truncation with moderate selection pressure ("truncation-moderate")
Progeny from the previous cycle were ranked on economic index, and the top-ranked 100 progeny out of 1,000 (10% selection proportion) were used in 50 random pair-wise matings with 16 S 0 progeny per mating (800 S 0 progeny) and two selfs per parent plant (200 selfs) for a total of 1,000 progeny. HST 30 was the priority trait for independent culling (see below).

| Truncation with high selection pressure ("truncation-high")
Progeny from the previous cycle were ranked on economic index, and the top-ranked 40 progeny out of 1,000 (4% selection proportion) were used in 20 random pair-wise matings with 40 S 0 progeny per mating (800 S 0 progeny) and five selfs per parent plant (200 selfs) for a total of 1,000 progeny. HST 30 was the priority trait for independent culling (see below).

| OCS and moderate selection pressure ("OCS-moderate")
Economic index of progeny and their relationships from the previous cycle were used to optimize a design for 50 matings (the same number of matings as in truncation-moderate) based on a balance strategy of target 45 degrees, where 0 degrees gives full emphasis to short-term genetic gain in index and 90 degrees gives full emphasis to minimizing parental coancestry, and a maximum of five matings per selection (Cowling et al., 2017), with 16 S 0 progeny per mating and two selfs per parent plant for a total of 1,000 progeny. This process was slightly modified for selection on HST 30 (see below).

| OCS and high selection pressure ("OCS-high")
Economic index of progeny and their relationships from the previous cycle were used to optimize a design for 20 matings (the same number of matings as in truncation-high) based on a balance strategy of target 45 degrees and a maximum of five matings per selection (Cowling et al., 2017), with 40 S 0 progeny per mating and five selfs per parent plant for a total of 1,000 progeny. This process was slightly modified for selection on HST 30 (see below).

| Selection for HST 30 score
HST 30 score was simulated as follows: for two days before and two days after anthesis, at least 30 plants of each genotype of the founder population were exposed to an effective temperature of 30°C in a controlled environment while keeping soil moisture at field capacity (Supporting information  Table S1). Effective temperature is the temperature half-way between the mean and the highest temperature each day for 4 days during anthesis (Deryng et al., 2014). In our simulation, an effective temperature of 30°C was achieved for 4 days during anthesis in a growth room which ranged from 21°C minimum, 27°C mean, and 33°C maximum, which represents a heat wave during flowering of wheat. This effective temperature will cause 50% yield loss in contemporary wheat varieties (Deryng et al., 2014), and contemporary wheat varieties are therefore allocated an average HST 30 score of 30.
The founder wheat population starts with a mean HST 30 score of 30 with a standard deviation of true breeding values of 1 (Figure 1). The starting narrow-sense heritability of HST 30 score is 0.3, and we assume HST 30 is not genetically correlated with any other traits in the selection index, assuming these traits are expressed in a non-heat stressed environment (Supporting information Table S2). We also assume HST 30 is not affected by genotype × environment interactions because it is assessed under controlled environment conditions. Under our simulated environment, a HST 30 score of 30 is recorded when 50% of the tested plants set viable seed at an effective temperature of 30°C, and the score ranges from 25 when no plants set seed to 35 when all plants set viable seed (Supporting information Table S1).
Three strategies for selection of HST 30 were implemented (a) no selection (ΔHST 30 score = 0), (b) moderate selection (ΔHST 30 score increases by +2 over 30 cycles), and (c) high selection (ΔHST 30 score increases by +4 over 30 cycles). For ΔHST 30 score +2, the target population mean HST 30 score was increased by 0.0667 in each cycle for a mean HST 30 score of 32 at the end of cycle 30. This increase in HST 30 allowed for the predicted increase in global average surface temperatures of +1.8°C by 2077 (average of all models) (IPCC, 2014) and we use +2°C since temperatures will be higher on land than on sea (IPCC, 2014). For ΔHST 30 score +4, the target population mean HST 30 score was increased by 0.1333 in each cycle for a mean HST 30 score of 34 at the end of cycle 30. This increase in HST 30 allowed for the predicted increase in global average surface temperatures of +3°C by 2077 (highest model) (IPCC, 2014) and we use +4°C since temperatures will be higher on land than on sea (IPCC, 2014).
HST 30 was the priority trait for independent culling. EBVs for HST 30 were sorted into descending order, and the lowest scoring candidate was discarded. This continued until the mean HST 30 in the remaining candidates at least equaled the target mean HST 30 for the next cycle. The remaining candidates then competed on index using independent culling as described above.
Under OCS, the overriding constraint was to achieve the predicted HST 30 score in the progeny generation to meet the target HST 30 score set for the year of their production, while economic index selection was optimized in the face of that constraint.

| Heat-adjusted grain yield
Based on Deryng et al. (2014), we calculate the heat-adjusted grain yield in each cycle as the predicted grain yield less 10% for every degree Celsius by which the achieved HST 30 falls below the target HST 30 in selected parents in the previous cycle (Figure 1). Future research may demonstrate a non-linear relationship between temperature and the effects of heat stress; in that case, an updated relationship between HST 30 and heat-adjusted grain yield (Figure 1) can be accommodated in the model.

| RESULTS
Under truncation selection, there was rapid early gain in index, but this slowed over 60 years as population coancestry increased and genetic diversity decreased (Supporting information Figures S2-S4). Truncation-high failed to reach the goal of +2 or +4 units in HST 30 in 2077 (Figure 2), despite the priority selection for HST 30 . In the absence of global warming, the economic index ( Figure 3) and grain yield (Figure 4) reached an early plateau and heat-adjusted grain yield began falling before 2077 due to a +2°C ( Figure 5) or +4°C ( Figure 6) increase in global temperatures. The plateau and fall in heat-adjusted grain yield before 2077 resulted from the failure of truncation selection to achieve target HST 30 and economic index, simultaneously, especially at high selection intensity.
Without selection for HST 30 over the next 60 years, under all scenarios there were major losses in heat-adjusted grain yield before 2077 due to a +2°C ( Figure 5) or +4°C ( Figure 6) increase in global temperatures, despite continued investment in wheat breeding.
In contrast, selection for HST 30 under OCS-moderate kept pace with target HST 30 levels (Figure 2), economic index nearly tripled its starting value in 2017 (Figure 3), and heat-adjusted grain yield in 2077 was more than double the founder population yield in 2017 after an increase in global temperatures of +2°C or +4°C (Figures 5 and 6). OCS-moderate achieved an average rate of increase in grain yield over 60 years of 1.8% per year, which is twice the rate of gain in wheat production expected in the 21st century (Alexandratos & Bruinsma, 2012), and double historical rates of improvement in yield of wheat (<1% per year) (Brancourt-Hulmel et al., 2003;Robertson, Kirkegaard, Rebetzke, Llewellyn, & Wark, 2016).
There was a small "cost" of increasing HST 30 by +2 or +4 units by 2077 in the absence of global warming, observed as a slightly reduced economic index and grain yield compared with no selection for HST 30 (Figures 3 and 4), but major potential benefits when HST 30 was improved sufficiently to prevent catastrophic yield declines in the presence of global warming (Figures 5 and 6).

| DISCUSSION
Breeding for HST in the world's major food crops (wheat, maize, rice) is very important for improving future global food security (Atlin et al., 2017). Rapid cycles of crop breeding with frequent phenotyping for HST and adoption of new technologies were considered essential for crop improvement during global warming (Atlin et al., 2017). Our simulation of wheat breeding shows that rapid cycles of recurrent selection with OCS and moderate selection intensity is the best strategy in the long run to improve all traits in the economic index. In contrast, truncation selection failed to reach targets for HST and heat-adjusted grain yield began falling before 2077. A potentially serious disruption to global food security may be averted by adopting OCS with moderate selection intensity on an economic index, with priority selection for HST.
In self-pollinating grain crops, most breeders practice some form of independent culling (Bernardo, 2010). In our BLUP-based simulations based on the infinitesimal model, high intensity selection with truncation for HST and an economic index ("truncation-high") caused a rapid loss of genetic diversity in the population as population coancestry approached unity (Supporting information Figures S2-S4). Under high selection intensity, genetic drift caused a reduction in genetic gain in the long term, and reduced the viability of investment in crop breeding. OCS reduced the negative impact of drift, and OCS-moderate always achieved a better result in the long term than OCS-high.
Our focus is on the population improvement phase of crop breeding (Supporting information Figure S1), which is the most logical place to apply pedigree BLUP or genomic selection (Hickey, Chiurugwi, Mackay, & Powell, 2017), and selection is based on S 0 -derived S 1 family selection with 2year selection cycles. Cowling and Li (2018) modeled independent culling on near-homozygous S 5 -derived wheat lines following rapid single seed descent in 3-year cycles from 2017 to 2077, for the same traits, and compared independent culling to selection on economic index and OCS. All methods achieved target HST 30 of +4 (the first trait in order of independent culling) by 2077, but OCS on economic index was superior overall, and especially for grain yield which was last in the order of independent culling (Cowling & Li, 2018). This confirms the conclusions of Hazel and Lush (1942) and Pesek and Baker (1969) that index selection out-performs independent culling in self-pollinating crops.
In the famous 100-year University of Illinois maize experiment, under moderate selection intensity (approximately 20%) for higher oil or protein content, genetic gain continued continuously over 100 years of annual cycles of recurrent selection (Dudley & Lambert, 2004). The response largely followed predictions from the infinitesimal model, with potentially new (small effect) mutations contributing to the long-term gain in oil or protein content (Walsh, 2004). In contrast, modern crop breeding programs are normally based on high selection intensity for multiple traits by independent culling, and are therefore very vulnerable to the detrimental effects of genetic drift. Our models show that BLUP breeding with OCS based on an economic index improves outcomes from crop breeding for multiple traits, especially when HST must be added to meet the demands of global warming (Figures 2-6).
While we do not explicitly accommodate genotype x environment x management (GxExM) interaction in our models, we make allowance for this by using low starting values for narrow-sense heritability (h 2 = 0.3) for grain yield, disease resistance and other traits (Supporting information Table  S2). The relative ranking of the four breeding scenarios will not change if strong GxExM exists in the target region of the breeding program, or if patterns of GxExM change over time-OCS will out-perform truncation selection. If adverse genetic correlations are discovered between HST and other traits, or if new research suggests a different relationship between HST 30 and yield penalty (Figure 1), this will not change the relative merit of OCS compared to truncation selection. Similarly, the effects of climate change may require that economic weightings (Supporting information   Table S2) are adjusted in future, for example, to select earlier flowering types in a warmer climate, or to reduce weighting on disease resistance which is less economically important in a drier climate. It is also possible that HST 30 may be negatively correlated with other economic traits. The strategic implementation of OCS in Matesel can be adjusted for these changes as knowledge improves over time.
Two approaches to genomic selection have been proposed recently in wheat: (a) combining existing pedigree and genomic information in "single-step" genomic prediction (Ashraf et al., 2016), and (b) "parallel cycles" of recurrent selection based on genomic predicted breeding values (Gaynor et al., 2017). With parallel cycles, trained genomic markers are available from the previous phenotypic cycle, and selections are returned to the next phenotypic cycle. The accuracy of genomic EBVs decreases rapidly in each cycle until genomic markers are retrained (Gaynor et al., 2017). With single-step genomic prediction (Ashraf et al., 2016), all relationship information is included (pedigree and genomic relationships), so that selection accuracy increases in each cycle.
Genomic selection may accelerate genetic progress through shorter crossing cycles and higher accuracy of selection in individuals without records (Gaynor et al., 2017;Hickey et al., 2017), but efficient genomic selection depends on accurate phenotyping on training populations that are closely related to the commercial breeding program (Jonas & de Koning, 2016).
We deliberately based our model on spring wheat in a low-yielding rainfed environment, typical of many developing countries, because these environments offer the greatest opportunity for substantial increases in global food production (Tester & Langridge, 2010) and this is where food security is under greatest threat due to climate change. The cost of OCS or truncation selection is similar (1,000 progeny tested in each cycle), but the outcomes of each breeding program are strikingly different from 2017 to 2077 during a period of global warming (Figures 5 and 6). In our modeling, the wheat breeding option of OCS based on an economic index will achieve significantly higher rates of genetic improvement in grain yield, HST, and other economic traits during 60 years of global warming than the option based on truncation selection. OCS with moderate selection intensity was the only method to achieve a 70% increase in wheat yield by 2050 (from 1.50 to 2.55 t/ha) during global warming, as motivated by the 2009 World Summit on Food Security ( Figure 6).
ABLUP breeding with OCS is readily implemented into existing breeding programs, and is a low-cost approach with high potential value in developing countries, where it may improve grain yields during climate change. ABLUP can be readily "upgraded" to GBLUP or single-step HBLUP based on "single-step" genomic prediction as defined in (Ashraf et al., 2016). OCS will be important to optimize genetic progress with genomic selection and with new alleles from gene editing (Jenko et al., 2015;Yin, Gao, & Qui, 2017). OCS should help to reverse the trend toward lower heat tolerance in recent US wheat varieties (Tack, Barkley, & Nalley, 2015), and to reverse the trend toward narrow genetic diversity in many crop breeding programs (Cowling, 2013;Tack, Lingenfelser, & Jagadish, 2017). BLUP breeding with OCS will provide greater opportunity to select for HST and to adapt crops to climate change (Tack et al., 2015(Tack et al., , 2017, by maximizing the rate of genetic gain for a given "acceptable" rate of population inbreeding.