Separating heat stress from moisture stress: analyzing yield response to high temperature in irrigated maize

Several recent studies have indicated that high air temperatures are limiting maize (Zea mays L.) yields in the US Corn Belt and project significant yield losses with expected increases in growing season temperatures. Further work has suggested that high air temperatures are indicative of high evaporative demand, and that decreases in maize yields which correlate to high temperatures and vapor pressure deficits (VPD) likely reflect underlying soil moisture limitations. It remains unclear whether direct high temperature impacts on yields, independent of moisture stress, can be observed under current temperature regimes. Given that projected high temperature and moisture may not co-vary the same way as they have historically, quantitative analyzes of direct temperature impacts are critical for accurate yield projections and targeted mitigation strategies under shifting temperature regimes. To evaluate yield response to above optimum temperatures independent of soil moisture stress, we analyzed climate impacts on irrigated maize yields obtained from the National Corn Growers Association (NCGA) corn yield contests for Nebraska, Kansas and Missouri. In irrigated maize, we found no evidence of a direct negative impact on yield by daytime air temperature, calculated canopy temperature, or VPD when analyzed seasonally. Solar radiation was the primary yield-limiting climate variable. Our analyses suggested that elevated night temperature impacted yield by increasing rates of phenological development. High temperatures during grain-fill significantly interacted with yields, but this effect was often beneficial and included evidence of acquired thermo-tolerance. Furthermore, genetics and management—information uniquely available in the NCGA contest data—explained more yield variability than climate, and significantly modified crop response to climate. Thermo-acclimation, improved genetics and changes to management practices have the potential to partially or completely offset temperature-related yield losses in irrigated maize.


Introduction
Maize (Zea mays L.) is the leading cereal crop produced worldwide. Maize yields have steadily increased over the past century due to breeding and improved management (Duvick 2005). Demand for maize grain is expected to expand by 50% in the developing world by 2050 (Msangi and Rosegrant 2011) while, over the same period, mean global temperature is expected to increase by up to approximately 2°C (IPCC 2014). Statistical studies have indicated that daily maximum temperature greater than approximately 30°C limit maize yields Roberts 2009, Lobell et al 2011). These studies suggest that temperature mediated limitations on maize yields will increase with projected increases in growing season temperatures (Schlenker and Roberts 2009, Lobell et al 2011, Lobell et al 2013, Deryng et al 2014. High temperatures and high vapor pressure deficits (VPD)-an important driver of evaporative losses-are often associated with drought and low soil moisture, and VPD has been found to be a significant explanatory variable of maize yields (Roberts et al 2012, Lobell et al 2013. These studies either examined rainfed maize systems only or did not differentiate between irrigated and rainfed maize so that it is not possible to determine if there is a direct effect of high temperatures on maize yields independent of soil moisture stress. Because specific correlation structures between temperature and moisture will likely not be conserved under climate change (Trenberth et al 2014, Williams et al 2014, accurate yield projections made using future climate scenarios require separately quantifying the mechanisms by which high temperature and soil moisture stress impact yield. One direct approach to do this is to analyze irrigated maize where soil moisture stress is minimized. Prior analysis of irrigated maize versus rainfed maize (Shaw et al 2014) found that exposure to high temperatures (>29°C) reduced some rainfed maize yields but not irrigated maize yields. This analysis (Shaw et al 2014) examined only a limited number of sites and used county-level yield data with no information on cultivar, planting date and specific management practices (e.g., planting rate, planting date, tillage, previous crop).
Recent field and laboratory studies have reported disruption of physiological and developmental processes in irrigated maize exposed to high temperatures, indicating that high temperatures may impact maize yields by pathways independent of associated moisture stress (Schoper et al 1987, Dupuis and Dumas 1990, Wilhelm et al 1999, Cicchino et al 2010, Rattalino et al 2011. These studies indicate threshold temperatures for damage of approximately 32°C-34°C. However statistical analyzes of reported irrigated maize yields in the US Corn belt over multiple climate years found little impact of high temperature on yield even when the reported temperature thresholds for damage were exceeded (Shaw et al 2014). These results, and a survey of high temperature thresholds for maize (Sanchez et al 2014), indicate that there is still uncertainty regarding these thresholds. In addition, both crop genetics (Duncan and Hesketh 1968, Duvick 2005 and crop management practices (Winstanley andChangnon 1999, Challinor et al 2014) may cause shifts or modulation of apparent high temperature thresholds. These data were often not available in previous statistical analyzes of regional high temperature impacts on maize, limiting the conclusions of those studies and constraining statistical analyzes of the relevance of such parameters over broad spatial scales.
There were two objectives of this analysis. The first was to analyze irrigated maize yield data for evidence of temperature-induced yield loss at high temperature thresholds reported in experimental trials. The second was to test if genetics and crop management offer a means to adapt to high temperatures. For our analysis, we obtained irrigated maize yield data from the National Corn Growers Association (NCGA) 'National Corn Yield Contest' from 2005 to 2012 for three states in the western US Corn Belt (Nebraska (NE), Kansas (KS), Missouri (MO)). These states have some of the highest growing season temperatures and VPDs in the US Corn Belt. The NCGA data set included grain yields (bushels acre −1 at 15.5% moisture, converted to Mt ha −1 ) from nearly 2000 entries (Methods). Critical management data (planting date, planting rate, cultivar brand and name, tillage, and previous crop) not often available in other studies were also included (Methods). These additional data allowed us to statistically control for management impacts on yield variance. In addition, with specific cultivar information available based on cultivar brand, we were able to map crop phenological growth stages over time for each yield contest entry whereas other studies have mainly used state-wide averages for characterizing crop phenology. By being able to isolate specific crop growth stages for each yield contest entry, we were able to analyze whether high temperatures coincided with specific developmental processes in maize.

Methods
Maize yield and climate station data Irrigated maize yield data from 2005 to 2012 were obtained from the NCGA National Corn Yield Contest (http://www.ncga.com/for-farmers/national-cornyield-contest) for NE, KS, and MO. Yield contest data included county location, maize yield, farm ID, planting date, planting rate, previous crop, tillage practice, and cultivar name (supporting data (SD) Methods figure S1). Contest grain yields are reported in bushels acre −1 (15.5% moisture) and were con verted to Mt ha −1 in our analyzes. Grain yields were measured with a calibrated yield monitor from a 0.5 ha section of a field that contest rules specify must be at least 4 ha in area planted to the reported cultivar. The NCGA has strict rules governing the contests, and verifies that all operations are done following these rules. All contest yields are verified by the NCGA. The high reported contest yields indicate near-optimal management in all contest entries. These near-optimal management conditions, together with the extensive data on management among entries, allowed for yield response due to climate to be analyzed relative to yield response due to management.
Controlling for temperature impacts on crop phenology Yield contest entries were associated with hourly climate data from 66 state mesonet stations associated with the Automated Weather Data Network of the High Plains Regional Climate Center and Missouri Mesonet. Since the reported entry location was only by county and no finer measure of location was provided, it was necessary to identify a weather station in close proximity to, or within, each county with yield contest entries. To do this, the average growing season growing degree day (GDD) accumulation for each year was calculated for each NCGA yield contest county using the nearest weather station from a list of stations within 200 GDD of the associated county mean (SD Methods, figure S2). This associated weather station was also used as a proxy for location in the mixed effects models described below.
In the US Corn Belt, maize phenological development is primarily a function of thermal-time (counted as GDD) (Bonhomme et al 1994). Planting date and product information on GDD to silking and GDD to physiological maturity for the cultivars in the NCGA data set were used to determine the cultivar-specific GDD requirement for each of the three crop growth stages: 'early growth' (EG), 'reproductive period' (Sen), and 'grain-fill' (G) (SD Methods, figure S3). These GDD data, combined with planting date and weather station associated with each entry, allowed us to calculate climate variables over both the entire growing season and by growth stage for each cultivar. Controlling for temperature impacts on crop phenology was important in this study for several reasons. First, knowing when a heat stress event occurred during phenological development provides evidence regarding the mechanism by which heat stress impacts yield. Second, there is a broad range of maturity classes and, therefore, thermal-time requirements for development among the cultivars in the NCGA yield contest. Given this broad range, it would be expected that yield contest entry-specific evaluation of growth development stage timing and duration would improve accuracy when determining whether high temperature events coincided with key growth development phases. Finally we were interested in identifying high temperature impacts on yield independent of temperature impacts on phenology. Capturing phenological impacts through climate variable processing allowed us to do this.
The temperature factors or indices examined in our study included 'killing degree days (KDD), cumulative degree-days (°C d) where daily maximum temperature exceeded 29°C (Butler and Huybers 2013), cumulative degree-hours over 30°C (TT30) and 34°C (TT34) (Cicchino et al 2010), and an additional index, average night temperature (AvgNT), which was the mean temperature between the hours of 7 pm and 7 am. We also modeled the impact of average maximum daily VPD and cumulative incident solar radiation (RAD) on irrigated maize yields. As irrigated canopy temperature and air temperature may differ, maximum daily canopy temperature was calculated using an empirically derived relationship between air temperature, VPD, and RAD from AmeriFlux data over irrigated maize in Meade, NE (Payero and Irmak 2006). Maximum daily canopy temperature was then used to calculate a 'canopy' KDD (Canopy T KDD) (SD Methods, figure S4).
Temperature indices, RAD and VPD were then calculated for each entry-specific thermal-time defined growing season and crop growth stage (EG, Sen, and G) durations. As we calculated all of our climate variables over thermal-time defined intervals, phenological impacts on cumulative intercepted radiation (RAD) were accounted for in our data processing (SD Methods figure S5), allowing us to look for other high temperature impacts on yields.

Statistical analyses
To prevent model scaling issues and allow for comparison of coefficients, RAD, AvgNT, yield, planting rate, planting date (as day of year), and cultivar GDD to maturity were standardized (z-values). Because of the significance of a 'zero' value of high temperature threshold indices in the context of this study, KDD, Canopy T KDD, TT30, and TT34 were log-transformed with zeros retained as zeros (McCune et al 2002). VPD was not transformed.
Yield and climate data were analyzed using a series of nested linear mixed-effects models (lme4 package in R) (McCune et al 2002). Candidate models (table 1) were designed to test hypotheses about specific mechanisms of maize heat stress response from the literature, while controlling for yield variance due to management, and quantifying remaining inter-farm, regional, inter-annual and genetic yield variance.

Base model
The 'base' model (SD Results table S1) contains fixed effects for planting date, planting rate, planting rate crossed with planting date, and cultivar GDD to maturity. All subsequent models contained these terms as well. In addition to management fixed effects, year of entry was also evaluated in the base model as a continuous variable to detect a linear trend in yields, but was found to be insignificant (p = 0.40).

Seasonal model
The seasonal models follow from recent statistical/ econometric analyzes Roberts 2009, Lobell et al 2013) that determine a seasonal coefficient for yield response to a climate variable (KDD or VPD) that is then used to project climate impacts on maize yields under different future climate scenarios. Yield response to climate over the entire growing season (S) was tested by inputting individual seasonal climate indices into the base model as a fixed effect (resulting in a series of 'season' models; one for each individual climate variable) and testing the significance by likelihood ratio testing relative to the 'base' (serving as the reduced) model. The Akaike Information Criterion (AIC) was calculated for the base model and each where the AIC from the radiation model is AIC min ) are reported in table 2 (Zuur et al 2009).

Crop growth stage model
The crop growth stage model was designed to test for possible high temperature impacts that have been reported for specific growth stages (Schoper et al 1987, Dupuis and Dumas 1990, Wilhelm et al 1999, Cicchino et al 2010, Rattalino et al 2011. To do this, yield response to high temperature was analyzed independently during EG, Sen, and G (SD Methods figure S3). In this analysis, we controlled for both increased phenological development associated with high temperature, and soil moisture stress. This was done by calculating all climate variables, including solar radiation, over entry-specific thermal-growth defined EG, Sen, and G and by examining only irrigated yield contest entries.
Since heat stress tolerance, and therefore heat stress response, may represent a genotype by environment (G x E) interaction (Raudenbush andBryk 2002, Sung et al 2003), we introduced terms to test whether exposure to high temperatures during the prior crop growth stage(s) (EG and Sen) modified yield response to high temperature during Sen and/or G.

Management model
The management model was designed to evaluate whether management modified yield response to these specific temperature response mechanisms. Crossed effects between planting rate, planting date, cultivar GDD to maturity, and climate variables calculated by crop growth stage (EG, Sen, and G) were introduced.
The 'null' model (random effects only, table 1) was fit using the restricted maximum likelihood procedure. All random effects were evaluated using likelihood ratio Table 1. Diagram of nesting and structure of linear mixed effects models used in analysis.

*
The Null and Base models nest within all other models. The Seasonal models are compared to the base model (likelihood ratio testing) and each other (using AIC). The Crop growth stage model and Management model nest with each other, but not the Seasonal models. Formula symbols follow notation used in the lme4 package in R (Bates et al 2014 lme4: linear mixed-effects models using Eigen and S4. ArXiv e-print; submitted to J. Stat. Soft., arXiv: 1406.5823). Fixed effects terms listed in bold were determined to be significant following backwards elimination with likelihood ratio testing (Ramirez-Villegas et al 2015).
testing (where a full model was compared to a reduced model which did not contain the variable of interest). All other models were fit using the maximum likelihood procedure, and fixed effects were selected by backwards elimination using likelihood ratio testing (Zuur et al 2009). Information on eliminated variables can be found in SD Results (tables S2 and S3). Variance inflation factor (VIF) was tested on each term in each model to prevent errors from multi-collinearity in predictor variables; no VIFs greater than three were observed. All crop growth stage tiered models ('Base', 'Crop growth stage' and 'Management') were compared using AIC values to prevent over-parameterization. The amount of inter-annual and location-based yield variance explained by the inclusion of significant climate variables in the model was evaluated by calculating the percent reduction in the year-level and location-level variance components of the crop growth stage model from those corresponding components in the base model (Raudenbush and Bryk 2002). Fixedeffects level yield variance was calculated used a hybrid of two methods (Raudenbush andBryk 2002, Nakagawa andSchielzeth 2013). Fixed effects were partitioned into groups based on interaction terms. These groups of fixed effects were inserted into the null model one group at a time, and the fixed effects level variance (Nakagawa and Schielzeth 2013) was determined for each test model.

Seasonal analyzes
Growing season daytime temperature stress indices (KDD, Canopy T KDD, TT30 and TT34) and VPD were not significantly related to maize yields (table 2). This was not due to unusually low growing season temperatures. The study period included a broad range of seasonal average maximum temperatures and seasonal RAD, as well as significant inter-annual variability in irrigated maize yields (figures 1(a) and (b)). For example, the average KDD in 2012 across all farms was approximately 350, well within the range that has previously been reported to reduce maize yields (Butler and Huybers 2013). As RAD was calculated by entry-specific thermal-time defined intervals, high growing season RAD could either be indicative of a long growing season resulting from below average temperatures (as in 2009), or overall high growing season RAD (as in 2012) (figures 1(a) and (b)).
Furthermore, in agreement with recent literature (Butler and Huybers 2013), we found little spatial structure to the yield values even though high temperatures, measured as degree days with a daily maximum temperature in excess of 30°C, followed a distinct latitudinal gradient (figures 2(a) and (b)).
While growing season daytime temperature indices and VPD were not related to yields, we did find that the coefficients of RAD (p-value = 5.4 × 10 -6 ), and AvgNT (p-value = 2.5 × 10 −4 ) were both significant (table 2). Effect sizes of RAD and AvgNT are approximately equal in magnitude (estimates of 0.12 and −0.12, respectively). Cumulative incident solar radiation and AvgNT are correlated in our data set (p-value of a bivariate fit <2 × 10 −16 , R 2 = 0.44) (figure 3). Because of the GDD model (34°C/8°C) we used to calculate the calendar time for crop development over the entire growing season and by individual growth stages, increasing night temperatures >8°C were modeled to increase phenological development and, therefore, reduce the calendar time for the season and for each growth stage. Since RAD was calculated as a cumulative value based on the calendar time duration of the entire season and of each growth stage as determined by the calculated thermal-time, higher night temperatures were associated with reduced RAD over the growing season ( figure 3).
AvgNT did not significantly modify yield response to calculated RAD (p-value = 0.78), and the yield impact of AvgNT becomes statistically insignificant if RAD is accounted for in the model (p-value = 0.18). This suggests that the impact of night temperature on yield is primarily associated with increased development rates and, consequently, reduced RAD. The model with RAD had the most empirical support, as indicated by the lowest AIC (table 2).

Crop growth stage analyzes
Yields were also compared to temperature stress indices by crop growth stage (EG, Sen, and G). Coefficients for KDD, Canopy T KDD and TT30 were not significant for any of the three crop growth stages (α = 0.05), including during Sen when many physiological processes contributing to yield formation in maize are active. We then examined RAD and TT34 for each of the three crop growth stages, as well as interactions among EG, Sen, and G TT34 (SD Results figure S6, table S2). We found that TT34 does significantly interact with yield, but that yield response to TT34 is complex. High temperatures during G slightly but significantly reduce yield gains that would be expected from increased  RAD during G (p-value = 0.008), but yield response to G TT34 was modified by EG temperature regime, as indicated by a significant crossed effect between EG and G TT34 (SD Results figure S6, table S2). When high values of TT34 occurred as isolated events during either EG or during G, yield is slightly but significantly reduced. However, yield response was both significant (p-value = 4.7 × 10 −4 ) and positive when both EG and G TT34 are high (figure 4).

Yield variance components
To better understand the factors controlling yield variability in irrigated maize, we examined the amount of yield variance explained by climate and management (figure 5). All significant climate interactions explained only a small amount of yield variance (less than 6% of residual variance; SD Results table S4), especially relative to yield variance explained by significant crop management parameters (Cultivar GDD to Maturity, Planting Rate and Planting Date; figure 5; SD Results table S4). Selection of longer season cultivars was associated with more positive response to high temperatures both in EG and G. While earlier planting dates were associated with increased temperature sensitivity during G (perhaps because of an association between earlier planting and cooler conditions in EG), yield gains associated with earlier planting were often sufficient to offset any yield losses from high temperatures, especially if coupled with a higher planting rate.

Adaptation to higher growing season temperatures
We were able to estimate the potential benefits of farmer adaptations relative to the estimated climate impact on irrigated maize yields by comparing 'typical' management versus 'adaptive' management values (figures 6(a) and (b)). 'Typical' management includes the mean planting date, the mean planting rate, and an average-performing cultivar of the median maturity class ( figure 6(a)). 'Adaptive' management considers a cultivar of an above-average maturity class (0.5 standard deviations longer than the mean), an earlier  than average planting date (0.5 standard deviations earlier than the mean), and an above-average planting rate (0.5 standard deviations above the mean) ( figure 6(b)). With these slight changes in management parameters that are associated with 'adaptive' management, our model projects yields that are higher at all Figure 5. Yield variance estimates calculated from significant random and fixed effects variance components. The random effects (other genetics, inter-annual variance, farm-level variance, residual variance) are shown in orange and significant fixed effects (previous crop, cultivar GDD to maturity, location-level, RAD and TT34, and planting rate and planting date,) are shown in blue (SD Results tables S1 and S3). n = 1929. temperatures compared to 'typical' management (figures 6(a) and (b)). The highest yields for 'adaptive' management are seen during years with high EG and G TT34. This suggests that, even in the highest yielding systems, small, agronomically plausible changes in current cultural practices in the US Corn Belt may be sufficient to offset yield declines in irrigated maize projected from increased air temperatures alone.

Discussion
While recent research has indicated that projected declines in maize yields due to high temperatures (Schlenker and Roberts 2009, Lobell et al 2013, Deryng et al 2014 are also likely linked with moisture stress (Roberts et al 2012, Lobell et al 2013, few studies have separated the effects of temperature and moisture stress. In this study, we analyzed irrigated maize in order to control for moisture stress and isolate the effects of temperature on yields. Our analyses indicate that irrigated maize production appears to be more limited by management and RAD than by high temperature. No significant irrigated yield response was observed to high temperature stress indices calculated over the entire growing season for a broad range of seasonal temperatures experienced in the US Corn Belt. This was not due to unusually low temperatures-record or near record seasonal average temperatures were reported in 2012 in many locations in the study area (State of the Climate (NOAA) 2012).
There were statistically significant yield responses to temperature stress indices when examined by crop growth stage. However, these were very weak and, if temperatures were consistently high throughout the growing season, positive. Grain-fill (G) appears to be the most sensitive crop growth stage for high temperature impacts on irrigated maize yields, where high temperatures appear to only slightly offset yield increases from increased RAD. But a direct negative impact of high temperatures during G is only evident if temperatures were low during EG. In contrast, there was a positive yield response to temperature when the crop experienced above-average temperatures in both EG and during G This finding is consistent with reported acclimation to high temperature in maize (Crafts-Brandner and Salvucci 2002, Sinsawat et al 2004, Butler and Huybers 2013, and appears to be a trait that is conserved across maize cultivars (results not shown). A model including the high-temperature acclimation response explains about 30% of the spatial variability in the yield contest data (SD Results table S4), suggesting that it may be a component of previously identified spatial acclimation to high temperature stress (Butler and Huybers 2013). Acclimation to heat stress, or acquired thermo-tolerance, in maize cropping systems may be an important consideration for predicting yield response to shifting temperature regimes.
Although many previous studies have focused on temperature impacts during reproduction, high temperatures occurring during the reproduction growth stage (Sen) were not significantly related to yields. This lack of response in the current population of maize cultivars we examined is consistent with findings of similar studies on irrigated maize (Rattalino Edreira et al 2011, Lobell et al 2013, and may be due to breeding advances that have resulted in increased stress tolerance in maize (Duvick 2005).
High night temperatures were associated with decreased yields. Since we modeled the duration of the crop growth stages using thermal-time, increased night temperatures led to decreased RAD per unit thermal-time. Decreased RAD as a result of accelerated development rates has been implicated as the cause of yield declines associated with elevated night temperature (Cantarero et al 1999, Hatfield et al 2014 and is the basis for shifting to longer season maturity classes both on-farm (Kaiser et al 1993) and in studies projecting the impact of farmer adaptations to climate change on maize yields (Chang 1981, Kaiser et al 1993, Tao and Zhang 2010, Moore and Lobell 2014. Other mechanisms for yield decline associated with high temperatures in the literature include increased rates of night respiration, but our data did not show evidence for a mechanism of night temperature-related yield decline other than accelerated phenological development, as the impact of night temperature became insignificant when radiation was accounted for in the model.
We also see evidence of adaptation to observed temperature impacts-both in maize management and in maize genetics-that are associated with yield responses of sufficient magnitude to offset any climate-related yield losses observed in our dataset (figures 6(a) and (b); SD Results table S3). Of particular interest are adjustments of planting dates and planting rates (Kucharik 2008, continued breeding for stress tolerance during the reproduction growth stage (Sen) (Duvick 2005) and the possibility of shifting to longer season cultivars in response to projected climate change (Kucharik 2008, Grassini et al 2011.
The lack of strong, negative response in irrigated maize yields to temperature despite record or nearrecord high air temperatures experienced in the study region in 2012 suggests that recently proposed heat stress thresholds (calculated from air temperature Cicchino et al 2010, Lobell et al 2011, Roberts et al 2012 are too low. Lack of temperature and VPD sensitivity in irrigated maize suggests that previously identified temperature/yield relationships in non-irrigated maize may have been the result of interactions between these variables and low soil moisture in a given climate system. As maize yield response to temperature is dynamic and highly leveraged by soil moisture status, additional experimental studies combined with carefully parameterized crop process models for separate and combined impacts of temperature and moisture stress may be necessary to accurately project yield trends and evaluate mitigation strategies under shifting temperature regimes.

Conclusion
Cumulative incident solar radiation had the greatest impact on irrigated maize yields in the western US Corn Belt for the study period (2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012). High daytime temperatures had a minimal negative impact on irrigated maize yields despite record or near-record high growing season temperatures in 2012. High night-time temperatures did negatively impact irrigated maize yields. However, our analysis suggests that this was due to increased crop development rates and, therefore, reduced RAD.
We found evidence of acquired thermo-tolerance in irrigated maize-high daytime temperatures early in the growing season followed by high daytime temperatures during grain-fill positively impact irrigated maize yields. This finding is consistent with reports of acquired thermo-tolerance in maize and takes on particular significance for projections of maize yield responses to climate change. These findings also suggest that current thresholds for heat stress in maize may be too low. We assessed high temperature stress indices that have been reported to negatively impact maize yields. However, these were not associated with decreased irrigated maize yields in the study region from 2005 to 2012.
Our results suggest that, for irrigated maize cropping systems, the combination of maize crop management (planting rate, planting rate) and appropriate maize genetics (e.g., selection of longer-season maize cultivars) can offset any climate-related yield losses observed in our dataset.