Responses of crop yield growth to global temperature and socioeconomic changes

Although biophysical yield responses to local warming have been studied, we know little about how crop yield growth—a function of climate and technology—responds to global temperature and socioeconomic changes. Here, we present the yield growth of major crops under warming conditions from preindustrial levels as simulated by a global gridded crop model. The results revealed that global mean yields of maize and soybean will stagnate with warming even when agronomic adjustments are considered. This trend is consistent across socioeconomic assumptions. Low-income countries located at low latitudes will benefit from intensive mitigation and from associated limited warming trends (1.8 °C), thus preventing maize, soybean and wheat yield stagnation. Rice yields in these countries can improve under more aggressive warming trends. The yield growth of maize and soybean crops in high-income countries located at mid and high latitudes will stagnate, whereas that of rice and wheat will not. Our findings underpin the importance of ambitious climate mitigation targets for sustaining yield growth worldwide.

where the suffixes r and t denote country and year, respectively; t r NFC , is the N fertilizer consumption per unit agricultural area (kg N ha -1 yr -1 ); The coefficients were determined country by country using the data for 1961-2012. These were significant in many crop-producing countries and indicated that the N fertilizer consumption increases with income growth and/or with the increase in the scarcity of per capita agricultural land in most cases (Supplementary Table S2).
In the second step, the country annual N fertilizer consumption per unit agricultural area ( t r NFC , ) was converted to the N application rate for a crop (kg N ha -1 yr -1 ): where the suffix p denotes the crop. In the third step, the country annual crop-specific N application rates ( where the suffix i indicates the region of the world (Africa, South America, Central America, North America, Europe, Asia and Oceania); ' , , p t r Napp is the N application rate for a crop, adjusted to be comparable with other data sources (kg N ha -1 yr -1 ); and ) were truncated to not exceed a certain level (50 kg N ha -1 yr -1 for soybean and 200 kg N ha -1 yr -1 for the remaining crops), as the literature has reported the plateau of the N application rate in developed regions 55,56 ; the N application rate for maize in the US has leveled off at approximately 150 kg N ha -1 yr -1 since 1980 (ref. 57). The use of these adjusting factors ensures that the N application rates used in this study are consistent with other data ( Supplementary Fig. S12).
In the final step, the adjusted country annual N application rate for a crop ( ) was translated into a grid-cell N application rate, where the suffix g indicates the grid cell; and The grid-cell N application rates were truncated at the aforementioned levels.
For the future period (2010-2100), we used the SSP country annual per capita GDP data. The agricultural area was kept constant at the 2010 level, whereas time-varying SSP population data were used to derive the country annual per capita agricultural area.
Knowledge stock of agricultural technologies. The country annual knowledge stock of the historical period was modeled using data on country annual GDP, agriculture value added and total R&D expenditure for the period 1960-2015 obtained from the World Development Indicators 50 . The future assumptions were derived using the SSP country annual GDP data 53 . The knowledge stock estimated here is not crop-specific, but in the latter procedure, it was translated into crop-specific assumptions regarding the use of improved technologies and management (see The use of improved technologies and management systems).
The country annual agricultural R&D expenditure, where t r GDP , is the gross domestic production (constant 2005 USD yr -1 ); t r AGV , is the agriculture value added (% of GDP); and t r RDE , is the total R&D expenditure (% of GDP). The comparison between the calculated agricultural R&D expenditure and reported data for the selected countries during the period 1980-2005 (ref. 58) showed good correspondence with relatively large errors in some developed countries ( Supplementary Fig. S13), ensuring that our assumption regarding the historical agricultural R&D expenditure is reasonable.
The literature 59 employed a knowledge stock accumulation equation based on the conventional capital accumulation equation. Another study 60 also used this equation to take the duration of research into account. Based on these studies, the country knowledge stock for a year, t r R , (constant 2005 USD yr -1 ), was calculated: where δ is the obsolescence rate of technological knowledge; and the lag time between research and technology adoption by farmers was set at six years. The value of δ was set at 0.1, as in ref. 60. Although lag time may vary by country and over time, other works use a similar order of lag time from seven to nine years 61,62 . Caution is necessary in the use of a constant lag time and obsolescence rate for long-term analysis, as it becomes increasingly difficult to find new ways to increase yields in regions where yields are already high 61 ; in the US, there are concerns that returns of agricultural R&D have been declining 63 ; in our modeling, the contributions from international agricultural research organizations to the knowledge stock in developing countries are found to be lacking.
In the future period, the SSP country annual GDP data were used. Agriculture value added and total R&D expenditure were kept constant at their mean 2006-2010 levels, although this assumption may be somewhat unrealistic because the GDP share of agriculture has decreased and that of total R&D expenditure has increased with economic growth 50 . Additionally, agricultural R&D expenditure has continued to increase at different rates in high-and middle-income countries 64 .
The use of improved technologies and management systems. The knowledge stock was translated into the use of improved technologies and management systems on farm fields. We assumed that improved varieties are central to new technologies but did not exclude other technologies, as the literature reports that improved varieties explain one-fourth to half of the past yield growth, and improved management, including the increased use of synthetic fertilizers, irrigation, chemicals and machinery, and improved input-use efficiency are the reasons for the remaining portion of the yield growth 65,66 .
In the crop model, we associated the knowledge stock with the curvature factor value that represents a crop's tolerance to stresses (heat, cold, water deficit and water excess) to represent the increased use of improved technologies and management systems with economic growth. Nitrogen deficit stress was also associated with the knowledge stock to represent the reported increase in N-use efficiency 55,56,67 . The form of the stress function (equations (27)-(31)) and the method of calibration were the same across the stress types and crops (see also Tolerance to stresses in Supplementary Note); hence, we used the heat stress as an example for explanatory purposes.
The heat stress factor for the i-th day after sowing, i γ heat, , is computed (equation (28)): where heat α is the curvature factor for heat stress; Ti is the daily mean temperature (°C); and To and Tu are the crop-specific optimal and maximum temperatures for growth (°C), respectively. The values of the stress factor ( i γ heat, ) range from 0 under the severely stressed condition to 1 under the optimal condition.
Here, we assumed that the mean heat stress factor of the crop duration ( t γ heat, ) can be approximated based on the mean temperature during the same period ( t T , °C): We also assumed that if the yield gap in a region is dominantly caused by a single stress type, then the mean stress factor due to the most dominant stress type should explain the yield gap. The yield gap used here was defined as the ratio between actual and attainable yields (Ya and Yp, t ha -1 , respectively). Using these assumptions, equation (12) is rewritten for a region where heat stress is dominant: Then, the value of the curvature factor for heat stress ( heat α ) is algebraically obtained: We associated grid-cell values of the curvature factor for heat stress with the knowledge  Table S4) and no adjustment for the knowledge stock; ii) Grid cells in which the modeled N deficit stress was minor were specified. The threshold for this was the 80th percentile value of the mean modeled N deficit stress factor for the period 1996-2005 over the global harvested area of a crop; iii) Grid cells in which the modeled cold stress was minor were specified using a certain threshold for cold stress (80th percentile), as was done for the N deficit stress. This step was repeated for the remaining stress types (water deficit and water excess); iv) The grid cells specified in step ii) and iii) were combined to select the grid cells in which all stresses other than the heat stress were minor; v) Values of the curvature factor for the heat stress in the selected grid cells were calculated according to equation (12") and associated with the mean knowledge stock values for the period 1996-2005 ( Supplementary Fig. S14); and vi) Steps ii) to v) were repeated for the remaining stress types to determine their coefficient values. For all crops and stress types examined here, negative slope values consistently emerged ( Supplementary Fig. S14), indicating that the curvature factor value decreases (i.e., the crop's tolerance increases) as the knowledge stock increases. These results are consistent with the fact that improved varieties often accompany higher tolerance to suboptimal conditions 66,69,70 . Although the values of the coefficients were not always significant because of the inconsistent spatial resolution between the two variables (i.e., grid-cell curvature factor versus country knowledge stock), we used the relationships addressed here in the crop model ( Supplementary Fig. S15).

Irrigation intensity.
The extent of crop-specific irrigated and rainfed areas was obtained from the Monthly Irrigated and Rainfed Crop Area around the year 2000 dataset referred to as MIRCA2000 (ref. 47). In addition, the annual growth rates of the area equipped for irrigation for 1900-2005 were obtained from the global historical irrigation dataset referred to as HID 71 . Although HID offers grid-cell estimates of the annual growth rate of an irrigated area, no crop-specific information is available. We therefore assumed that the annual growth rate of the irrigation intensity (the ratio of irrigated area to harvested area) was the same across the crops, although this assumption may be unrealistic for a certain location-crop combination.
The growth rates of irrigation intensity for 2006-2012 were extrapolated using a linear regression based on the data for the period 1998-2005. We assumed that the harvested area derived from MIRCA2000 did not change with time, but the irrigation intensity in the historical period did change with time according to HID. Thus, this study partially accounted for the contribution to historical yield growth from the increased irrigation intensity. In the future period, the irrigation intensity was kept constant at the 2010 level.

Supplementary Note
Crop model description A schematic overview of the CYGMA model is shown in Supplementary Fig. S11. The following sections describe the modeled processes of development, growth, yield formation, stresses and soil moisture in detail.

Development.
The crop total GDD requirement from sowing to harvest, GDDc (°C days), is calculated using the mean annual GDD with a base temperature of 0 °C, GDDa (°C days) 72 : For each year t, the mean annual GDD is calculated using the daily mean temperature over the last 10 years (t-10:t-1). a and b are the empirical coefficients (Supplementary Table S4).
The development of a crop for the i-th day after sowing, frGDD (0=sowing and 1=harvest), is expressed as a fraction of the crop growing season: where T is the daily mean temperature (°C); and Tu and Tb are the maximum and base temperatures for growth (°C), respectively (Supplementary Table S4). In this model, we only consider spring wheat (i.e., the vernalization in winter wheat is not considered).
Growth. The model simulates the daily growth of leaves and canopy height. However, the growth of roots is not considered, and the time-constant root depth is used in the model. The leaf area index is calculated based on ref. 73. Before leaf senescence becomes the dominant growth process: where ΔLAI is the increment of the leaf area index for day i (m 2 m -2 ); LAIi and LAIi-1 are the leaf area indices for days i and i-1, respectively (m 2 m -2 ); and frGDD, sen is the fraction of the growing season at which senescence becomes the dominant growth process (Supplementary Table S4). The daily increment of the leaf area index is given as follows: where frLAImax, i and frLAImax, i-1 are the fractions of the maximum leaf area index under the optimal condition for days i and i-1, respectively; LAImax is the maximum leaf area index (m 2 m -2 , Supplementary Table S4); and γ is the most dominant stress for day i (see Stresses in Supplementary Note). The fraction of the maximum leaf area index, frLAImax, is related to the fraction of the crop growing season: where l1 and l2 are the empirical coefficients (Supplementary Table S4). Once leaf senescence becomes the dominant growth process, the leaf area index decreases as the fraction of the growing season increases: where max I LA ′ is the maximum leaf area index achieved in a growing season (m 2 m -2 ). The canopy height of a crop is calculated 73 : where h is the canopy height for day i (m); hmax is the maximum canopy height of a crop (m, Supplementary Table S4); and frLAImax is the maximum leaf area index corresponding to a given fraction of the growing season. Once the maximum canopy height is reached, the canopy height remains constant until harvest.
Yield formation. Once the crop reaches anthesis, the model starts to calculate the amount of the harvestable part of a crop (yield): where yldi and yldi-1 are the yields for days i and i-1 (t ha -1 ), respectively; and Δyldi is the increment of yield for day i (t ha -1 ). The daily increment of yield is calculated: where Δbio/1000 is the increase in crop total biomass for day i (t ha -1 ); frΔbio is the fraction of the increment of the total biomass allocated to the harvestable part of a crop; and γ is the most dominant stress for day i. The partition fraction for the harvested part is calculated 74 : where frGDD, ant is the fraction of the growing season at which anthesis occurs; and p1 and p2 are empirical coefficients (Supplementary Table S4). The increment of the potential total crop biomass is calculated using the radiation-use efficiency for year t, RUE ((kg ha -1 ) (MJ m -2 d -1 ) -1 ), and the intercepted photosynthetically active radiation for day i, PAR (MJ m -2 d -1 ): . (24) The intercepted photosynthetically active radiation is calculated 73 : , (25) where SR is the downward shortwave radiation (MJ m -2 d -1 ); and k is the light extinction coefficient (Supplementary Table S4). The relationship between the annual mean CO2 concentration for year t, CO2 (ppm), and radiation-use efficiency is modeled 73 : where r1 and r2 are empirical coefficients (Supplementary Table S4).
Stresses. The model accounts for five different stress types: N deficit, heat, cold, water deficit and water excess. These stress factors vary between 0 under the severely stressed condition and 1 under the optimal condition. The heat, cold, water deficit and water excess stresses are calculated day by day, whereas the N deficit stress is computed on an annual basis due to the lack of N application rate data at a finer temporal resolution. The N deficit stress factor for year t, γNdef, is calculated: where Napp is the N application rate (kg N ha -1 year -1 ); Nappo is the optimal N application rate; Nappmin is the shape parameter used to determine the minimum N deficit stress; and Ndef α is the curvature factor for the N deficit stress that represents the crop's tolerance to the N deficit stress (Supplementary Table S4). The N deficit stress changes with the annual update of the N application rate, but it is constant across days within a growing season. A similar simplification is found in ref. 75.
The heat stress factor, γThigh, is a function of the daily mean temperature (T, °C), the optimal temperature (To, °C), the maximum temperature for growth (Tu, °C) and the curvature factor for heat stress ( heat The relationship used to calculate the cold stress factor, γTlow, is as follows: where Tb is the base temperature (°C) and cold α is the curvature factor for cold stress. The form of the function used here is the same as that used in ref. 73; however, the heat and cold stresses are separately modeled in this study.
The water deficit stress factor, γWdef, is modeled using an actual-potential evapotranspiration ratio, Ea/Ep, and the curvature factor for the water deficit stress, Under the irrigated condition, water is applied on day i+1 to increase the root-zone soil moisture to 80% of the plant-extractable soil water capacity when the water deficit stress factor for day i is below 0.9.
The relationship used to calculate the water excess stress factor, γWexs, is as follows: where SW is the root-zone soil moisture for day i (mm); SWmax is the plant-extractable water capacity of the soil (mm); and Wexs α is the curvature factor for the water excess stress. The water excess stress factor is calculated after three consecutive days of continued soil saturation (soil moisture exceeds 90% of the plant-extractable water capacity of the soil), as in ref. 76. The form of the function for the water deficit and excess stresses is based on ref. 77.
The calculated stress factors due to the N deficit, heat, cold, water deficit and water excess are adjusted for the knowledge stock to account for the use of improved technologies and management systems in developed countries. Taking the water deficit stress as an example, the adjustment is conducted through the modification of the curvature factor value: α is the curvature factor for the water deficit stress for a grid cell and  Table S4). The calculated curvature factor value is truncated to not exceed the lower and upper bounds (Supplementary Table S4), ensuring the biological limits of crop growth ( Supplementary Fig. S15). Given the form of parameterization, the crop's tolerance to the water deficit increases as the knowledge stock increases. The stress factors associated with the N deficit, heat, cold and water excess are adjusted in a similar manner. The most dominant stress for day i, γ, is then selected to decrease the daily potential increment of leaf area and yield: Soil moisture. The root-zone soil moisture is calculated based on the water balance equation used in ref. 13, which takes rainfall, snow melt, evapotranspiration, runoff, ground water loss through deep percolation and plant-extractable water capacity of the soil into account. The rainfall-snowfall separation and snow accumulation, melt and sublimation are calculated using the daily maximum and minimum temperatures and precipitation as the inputs of the snow cover submodel 78 . Surface runoff, subsurface runoff and ground water loss are calculated based on ref. 79 after revising this submodel to operate with a daily step instead of the original monthly step. Potential evapotranspiration is calculated by using a variant of the Penman-Monteith equation 73 . Actual evapotranspiration is determined by comparing the root-zone soil moisture and potential evapotranspiration. The plant-extractable water capacity of the soil data estimated using the soil texture, soil organic content, and plant-root (or soil profile) depth were obtained from ref. 46.

Calibration
Crop total GDD requirement. We addressed the relationships between the mean annual GDD for the period 1996-2005 and the crop total GDD requirement in 2000, as in ref. 72. Sowing and harvesting dates were obtained from MIRCA2000 (ref. 47) with the assumption that the sowing (harvesting) date was the middle day of the reported sowing (harvesting) month. The daily mean temperature data obtained from S14FD 41 were used to calculate the GDD. The crop calendar information provided by MIRCA2000 (ref. 47) is crude (a single value was often assigned for all of the cropland in a country). Therefore, for the irrigated cropping system, we selected a single grid cell with the most extensive irrigated area of a crop for each country. The same calculation was conducted for the rainfed cropping system. Using the grid-cell sowing and harvesting dates and temperature data, the relationship was specified for each crop and cropping system ( Supplementary Fig. S16).
Sowing date. The modeled sowing date changes on an annual basis with any update of the climate bin because refs. 72 and 80 suggested the importance of the temperature and moisture regimes in determining the sowing date. For the calibration, first, we calculated the climate bin using the 10-year mean annual GDD and aridity index (annual precipitation divided by annual potential evapotranspiration of reference crop, P/PET 81 ) for the period 1996-2005 ( Supplementary Fig. S17). Second, the mean sowing date was calculated for each climate bin and cropping system (irrigated or rainfed) using the reported sowing date in 2000. Third, for the irrigated cropping system, we selected a single grid cell with the most extensive irrigated area of a crop for each country (as presented in Crop total GDD requirement in Supplementary Note) and recorded the corresponding climate bin. Fourth, the mean sowing date across countries was calculated for each climate bin after separating the data from around the world into data from the northern and the southern hemisphere. Finally, if the sowing date for a climate bin was lacking, it was extrapolated using the data in the four neighboring climate bins ( Supplementary Fig. S18). The same calculation was performed for the rainfed cropping system. This calibration led to reasonable sowing dates similar to but more spatially detailed than those of MIRCA2000 (ref. 47) (Supplementary Fig. S19).
Tolerance to stresses. Each crop's tolerance to stresses is modeled to change as the knowledge stock increases. The relationship between the tolerance (represented by the curvature factor values in equations (27)-(31)) and knowledge stock was specified for each crop and stress type using the spatial data in 2000 ( Supplementary Fig. S14). This indicates that we applied the relationship addressed using the difference across the developed and developing countries in 2000 to our outlook. We assumed that the developing countries in 2000 grew their economies according to the socioeconomic assumption and improved the technologies and management used in their crop production systems to be more tolerant, as seen in developed countries in 2000. More details on the calibration are presented in The use of improved technologies and management systems in Supplementary Methods.

Limitations of the crop model
The model used in this study does not necessarily represent all characteristics of plant breeding and improvements in agronomic management during the last few decades. The increased N fertilizer input is largely due to the characteristics of modern semi-dwarf varieties, which can respond better to high N input than traditional tall varieties to become liable to load and to contribute to yield growth 66 . We considered the increased N fertilizer input but not other important characteristics of modern semi-dwarf varieties, such as increased harvest index and shortened plant height 66 . Additionally, the contributions of farm machinery to weed and pest control, timeliness of planting and harvest efficiency to yield growth 82 are lacking in the model.
We considered the increased crop tolerance. This is reasonable because new varieties often accompany higher tolerances to suboptimal conditions 66,69,70 . However, the consideration of some changes in agronomic management that accompany new varieties is lacking. For instance, during the last 50 years, plant density in the American Midwest has increased because of biotechnology-driven increase in crop tolerance to pest, disease, drought, and herbicide and the improvement in soil management 70,83 . Furthermore, the model does not fully consider changes in some traits that, in reality, accompany plant breeding for higher yield. If maize is taken as the example, such traits include ear height (there is a reducing trend in ear height), leaf angle (a more upright leaf orientation), staygreen (delayed leaf senescence), grain-filling period (newer hybrids have a longer period of grain fill with faster final dry-down rate), and kernel weight (increased weight per kernel) 70 . Additionally, plant breeding offers efficient light capture and harvest indices close to their theoretical maxima, leaving the efficiency of conversion of the intercepted light into biomass as the only remaining major prospect for improving yield potential in the coming decades 84 ; however, this view is not incorporated into our model and outlook.

Post-processing of crop model output
The model outputs from the historical run were compared with the reported yields 51 . We calculated the modeled global and country mean yields by combining the simulated grid-cell yields of irrigated and rainfed cropping systems. In the model, the harvested year could vary by cropping system because a different cropping system obtained from MIRCA2000 (ref. 47) sometimes indicated a different cropping season. This is the case, for example, in the tropics, where monsoons play a critical role in determining water availability. For a consistent comparison, the following procedure was used: i) We averaged over the simulated grid-cell yields produced using the irrigated and rainfed cropping systems, with both being harvested in the same year in the model, using the irrigated and rainfed areas in 2000 (ref. 47) as the weights; ii) The calculated grid-cell cropping-system mean yields were averaged using the grid-cell harvested area in 2000 (ref. 68) as the weight to obtain the country mean yield; iii) Global mean modeled and reported yields were calculated from the country mean yields using the FAO country harvested area 51 as the weight; and iv) We averaged the modeled and reported global and country mean yields over two consecutive years (t-1 and t) to obtain the 2-year moving averaged data. As noted in ref. 85, the FAO's definitions state that "yield data are calculated by dividing production by harvested area" and that "when the production data available refers to a production period falling into two successive calendar years and it is not possible to allocate the relative production to each of them, it is usual to refer production data to that year into which the bulk of the production falls". Therefore, a year-to-year comparison between the reported and modeled country mean yields is difficult to justify, except for countries where a single annual harvest is dominant. The 2-year moving averaging ameliorates the errors to some extent. For the future and noCC runs, however, the 2-year moving average was not used because decadal mean modeled yields were analyzed. The harvested area used as the weights was kept constant at the 2010 level for the future and noCC runs.

Global temperature change and cumulative CO 2 emission
We reconstructed the relationship between the cumulative total global CO2 emission from 1870 and the global decadal mean combined land and ocean surface temperature anomaly relative to preindustrial levels (1850-1900) ( Supplementary Fig. S3). Cumulative CO2 emission was computed using the reported emission estimates from fossil-fuel burning, cement manufacture and gas flaring for the period 1870-2009 (ref. 86) and RCPs for the period 2010-2100 (ref. 19). The global temperature change that occurred for the period 1961-2009 was calculated using S14FD 41 with the assumption that the global temperature change for the period 1986-2005, relative to 1850-1900, was 0.61 °C 87 . For the future period (2010-2100), the bias-corrected data of five GCMs 41 (see Climate and other inputs in the main text) were used. The overall relationship represented by the regression line in Supplementary Fig. S3 shows 1 °C and 5 °C warming from preindustrial levels for cumulative CO2 emissions of 1,238 GtCO2 and 7,652 GtCO2, respectively.

Assumptions regarding future agricultural technologies and management
For all four crops, it was assumed that the future N application rates would gradually increase with income growth and scarcity of per capita agricultural area, with variations by country and SSP, and the increase in the N application rate would slow down or approach a plateau, as exemplified by China and India ( Supplementary Fig. S20). These patterns are basically consistent with the assumption used in other works 55 . Interestingly, the difference in N application rates between SSP1 (rapid technological change) and SSP2 (intermediate technological change) was negligible, whereas the N application rate under SSP3 (slow technological change) was always lower than those under SSP1 and SSP2. However, caution is necessary because, in reality, the N application rate begins to gradually decrease after the plateau, reducing the environmental stress or producing high-quality grains at the expense of yield.
The knowledge stock of agricultural technologies assumed in this study rapidly increased with economic growth, with the exception of China for the period 2070-2100 ( Supplementary Fig. S21). This tendency was consistent across SSPs. The differences in the knowledge stock across SSPs were more prominent than those in the N application rate. Although the knowledge stock in developing countries will rapidly increase in the future, the knowledge stock in these countries (except China and India) in 2100 will still be lower than that of the US in 2010.

World
Bank. World Development Indicators.
Supplementary  6 For maize, rice and wheat, the default curvature factor values for N deficit were determined to follow the shape of the stress function presented in the literature. The default values for soybean were determined with the assumption that N stress is roughly half that of other crops. 7 The default curvature factor values for water deficit and water excess were determined to follow the shape of the stress functions presented in the literature.
Supplementary Figure S1. Global and country mean yields of four crops for selected major crop-producing countries for the period 1961-2012, calculated using the modeled and FAO-reported data. The correlation (r), p-value (p) and root-mean-square error relative to the mean reported yield for the period 1961-2012 (RMSE) are also presented.
Supplementary Figure S2. Summary of the comparisons between the modeled and reported yields for four crops. A box plot in the left panel summarizes the correlations calculated between the modeled and reported data for the period 1961-2012 (the sample size varies by country but is greater than 20) across countries; the vertical line of a box plot indicates the 90% probability interval; a box indicates the 50% probability interval; a horizontal line in a box indicates the median; a triangle on the right-hand side of a box plot indicates the correlation for the global mean yield; and the numbers below a box plot indicate the number of countries examined and the number of countries with a significant correlation. The box plot in the right panel summarizes the root-mean-square error relative to the mean reported yield for the period 1961-2012 (RMSE), as in the left panel.
Supplementary Figure S3. Global decadal mean surface temperature anomaly relative to preindustrial time period (1850-1900) as a function of cumulative total global CO2 emissions from 1870. Some decadal means are labeled for clarity (e.g., 2100 indicates the decade 2091-2100). Colored solid lines with dots indicate the future projection represented by the ensemble mean for each RCP (calculated using 5 GCMs). The colored shaded area indicates the ensemble spread (from minimum to maximum) for each RCP. The black solid line indicates the historical past calculated by using the S14FD meteorological forcing dataset. The gray dashed line indicates a regression line fitted to the data from the historical past and four RCPs. See Global temperature change and cumulative CO 2 emission in Supplementary Note for more details. Supplementary Figure S11. Schematic illustrating the CYGMA global gridded crop model. A blue ellipse indicates a physical variable: Tave, Tmax and Tmin, daily mean, maximum and minimum temperatures, respectively; P, precipitation; SR, downward shortwave radiation; RH, relative humidity, WS, wind speed; CO2, atmospheric carbon dioxide concentration; GDD, growing-degree days and P/PET, aridity index. A green ellipse indicates a socioeconomic variable: R&D, research and development; T&M, technologies and management; and Airri and Arain, irrigated and rainfed areas, respectively.
Supplementary Figure S14. Relationships between the grid-cell curvature factor value and mean knowledge stock in 1996-2005 by stress type and crop. Orange dots indicate grid-cell data; the gray shaded area indicates the density of data; and the red line indicates the regression line fitted to the data. The asterisks indicate the significance: ***, 0.1%; **, 1%; and *, 5%.
Supplementary Figure S15. Models of five different stress factors for four crops associated with N deficit (Ndef), cold, heat, water deficit (Wdef) and water excess (Wexs). Tave indicates the daily mean temperature; Ea/Ep indicates the ratio of actual-potential evapotranspiration rates; and SW/SWmax indicates the ratio of root-zone soil moisture relative to the plant-extractable water capacity of soil. The stress factor ranges from 0 under the severely stressed condition to 1 under the optimal condition.  Table S4 for the values of a and b.
Supplementary Figure S17. An example of the climate bins calculated using the S14FD meteorological forcing dataset for the period 1996-2005. The maps presented here were created with the Generic Mapping Tools (GMT) 49 version 4.5.12 (https://www.soest.hawaii.edu/gmt/) using data sources described in Supplementary Note. Figure S18. Sowing dates of four crops specific to the climate bins, cropping systems (irrigated or rainfed) and hemispheres (NH, northern hemisphere; and SH, southern hemisphere) used in the model. A large tile indicates that a data value was calculated using MIRCA2000. A small tile indicates that a data value was extrapolated from the four neighboring climate bins.

Supplementary
Supplementary Figure S21. Assumptions regarding the knowledge stock of agricultural technologies for selected major crop-producing countries in the historical and future periods.