Abstract Crop yields are projected to decrease under future climate conditions, and recent research suggests that yields have already been impacted. However, current impacts on a diversity of crops subnationally and implications for food security remains unclear. Here, we constructed linear regression relationships using weather and reported crop data to assess the potential impact of observed climate change on the yields of the top ten global crops–barley, cassava, maize, oil palm, rapeseed, rice, sorghum, soybean, sugarcane and wheat at ~20,000 political units. We find that the impact of global climate change on yields of different crops from climate trends ranged from -13.4% (oil palm) to 3.5% (soybean). Our results show that impacts are mostly negative in Europe, Southern Africa and Australia but generally positive in Latin America. Impacts in Asia and Northern and Central America are mixed. This has likely led to ~1% average reduction (-3.5 X 1013 kcal/year) in consumable food calories in these ten crops. In nearly half of food insecure countries, estimated caloric availability decreased. Our results suggest that climate change has already affected global food production.

Citation: Ray DK, West PC, Clark M, Gerber JS, Prishchepov AV, Chatterjee S (2019) Climate change has likely already affected global food production. PLoS ONE 14(5): e0217148. https://doi.org/10.1371/journal.pone.0217148 Editor: Young Hoon Jung, Kyungpook National University, REPUBLIC OF KOREA Received: October 31, 2018; Accepted: May 6, 2019; Published: May 31, 2019 Copyright: © 2019 Ray et al. This is an open access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited. Data Availability: Data underlying the study are government reported data and are publicly available from the links in S1 Text. The authors confirm they did not have any special access to this data that others would not have. Yield data for the studied regions for the top four global crops are in S8 Table. Funding: DKR, PCW and JSG were partially supported by funding from the Gordon and Betty Moore Foundation and the Belmont Forum/FACCE-JPI funded DEVIL project (NE/M021327/1), and the Institute on the Environment. MC was supported by the Wellcome Trust, Our Planet Our Health (Livestock, Environment and People - LEAP), award number 205212/Z/16/Z. SC was partially supported by the National Science Foundation (NSF) grants DMS-1622483 and DMS-1737918. The funders had no role in study design, data collection or analysis, decision to publish, or preparation of the manuscript. Competing interests: The authors have declared that no competing interests exist.

Introduction Previous assessments of climate change impact on crop yields commonly combine future climate scenarios and process-based crop models to project future yields for a limited number of crops for 2050 or later [1–4]. At higher levels of warming, strong yield losses are predicted in lower latitudes especially for maize and wheat crops [2]. Although these results provide insights into long-term future changes, there are large uncertainties in both the modeled climate projections [5] and in the crop model parameters [6–8]. Hence, the distant time horizon, small number of crops, and coarse resolution limit the results’ utility for stakeholders and policy makers to develop goals and strategies. Assessing the impacts of recent climate change complements long-term forecasts and identifies which crops and places are already at greater risk. Since the 1970s, global surface temperature warmed at an average of 0.16°C to 0.18°C per decade [9], a rate higher than any period since the industrial revolution. During that same period, we find that the growing season temperature, over all harvested areas for the top ten global crops–barley, cassava, maize, rice, oil palm, rapeseed, sorghum, sugarcane, soybean and wheat–increased 0.5°C to 1.2°C (S1 Table; S1–S4 Figs). Growing season precipitation changes were more variable; from a decrease of 3.4 mm averaged over all sugarcane harvested croplands to an increase of ~19 mm averaged over all oil palm harvested croplands (S1 Table; S1–S4 Figs). Given the complexities in modeling crop growth response to regional variability in climate and management, previous process based crop modeling efforts have had difficulty reproducing historical crop yields for the few major crops generally studied such as maize, rice, and wheat [7, 10]. Empirical (statistical) global estimates of recent climate change impacts, which some studies have shown to perform equally well as process-based modeling in assessing future impacts on crop yields [11, 12], is available, but only at the coarse national scale [13, 14] or at selected locations, and only for the top few crops [4, 15]. To enable subnational analysis, we first constructed a database of harvested area, yield, weather, and climate statistics building on methods developed previously and accessing publicly available data from individual countries [16]. Crop statistics were compiled from 1974–2013 for ten crops across ~20,000 political units globally (S1 Text). These ten crops account for ~83% of global kilocalorie production from all croplands [17] and represent several major crop group types. (Crops were not distinguished between varieties and managements—see section 1 of S1 Text for more details). Weather and climate statistics (seasonal and annual, normal and extreme, precipitation and temperature) were calculated for each political unit using the CRU TS4.01 gridded global dataset ([18], see S1 Text for further information on how political unit level weather and climate information was constructed). Historical climate is defined as the 30-year average weather prior to 1974. Current climate is defined as the historical climate plus the addition of the linear trend of the weather for the 35 years ending in 2008, from the year 1974. We constructed statistical models relating the observed yields to observed weather at each political unit from 1974 to 2008. (The period 2009 to 2013 was set aside for out-of-sample cross validation of model predicted crop yields against observations, at each political unit). We used a time-series analysis (Methods below) in which the influence of technology and management changes were accounted for in linear and quadratic time terms, whereas the normal and extreme temperature and precipitation variations and their interactions were represented with linear and quadratic terms consistent with previous studies [16, 19–21]. The analysis included only those political units that had crop data for each year from 1974 to 2008. We base quantitative conclusions on models that are statistically significant (ANOVA F-statistic at p < 0.05 level) and are applied within the domain of observed variables used in model construction (Methods below; S5 and S6 Figs). We applied models in the political units in which they were constructed. We confirmed that yield predictions were made only with historical and current climate conditions that were within the range of the observed weather used for model construction (S6 Fig). Further, testing the model residuals showed Gaussian nature of residuals everywhere and their general white noise nature (absence of autocorrelation) (S7 Fig) indicating that these models were robust relative to data quality issues. Cross validation of the model predicted yields, against observations, showed low average (2009 to 2013) harvested-area-weighted errors globally (S2 Table), and generally also low errors at the political units for individual crops (Methods below and S5 Fig). The model coefficients of determination (R2) at the global level ranged from 0.76 (sorghum) to 0.87 (rice) (S3 Table; S8 Fig). The potential impact of climate change at each political unit is the difference in crop yield under current, and under historical climate conditions. This yield change information at the political unit was then converted to production changes using the average harvested area information over the 2003 to 2008 period at each political unit. Then all the political unit level climate-driven production changes in a country were summarized to give the country-level climate-driven production changes. Finally, we translated the country-level production changes to country-level consumable food calorie changes. In this process we accounted for harvested calories that returned as food calories after getting processed through animals as feed, through food processing, and calories directly consumed with little to no processing (S1 Text).

Discussion and conclusions The Intergovernmental Panel for Climate Change—Assessment Report (AR) 5 [2] (IPCC AR5) on climate change impact on crop yield/production notes that between the AR4 and AR5 the required connection of climate change impact to food security impact was missing. Our study directly addresses this by translating the potential impact of recent climate change on crop yields (Fig 1) to consumable food calories change in each country (Table 1, S4 Table). While translating crop production change to consumable food caloric change, we accounted for the current dietary consumption pattern of individual crops per country, including bookkeeping for directly and indirectly consumed calories in each country [34]. We found that out of the studied 53 countries with a hunger index of serious, alarming, or seriously alarming [32] in 2008, 27 countries (or ~51%) had decreased consumable calories due to mean climate changes (~ -0.4% in these ten crops or ~ -0.3% across all food calories consumed across all 53 countries studied–S4 Table). Though we detected subnational level crop yield and production changes (Fig 1), determining consumable food caloric changes at that level would require data not available globally: subnational dietary patterns, evaluation of the climate change impact in the entire food supply chain [35] and socio-economic conditions [36]. These individual issues should be explored in future studies to understand mean climate change impacts on local scale food security. Linear and extreme trends in weather were captured using linear and quadratic terms whereas correlations in warmth and moisture were captured using interaction terms. We did not explicitly model increased atmospheric carbon dioxide (CO 2 ) impact as CO 2 is (1) highly correlated with time [11] and further (2) several lines of investigations [37] show that the science of CO 2 effect on crop yields is not settled. Impacts from changes in other variables were beyond the scope of this study and should be considered in future studies. However, high model R2 values, and the Gaussian nature and white noise of model residuals (indicating that residual errors were normally distributed and not auto-correlated), and cross validation against the yields of 2009 to 2013 showing low errors show that models were robust enough to answer the questions on mean precipitation and temperature change on crop yields (Fig 2; S5–S8 and S11 Figs and S2, S3 and S6 Tables). PPT PowerPoint slide

PowerPoint slide PNG larger image

larger image TIFF original image Download: Fig 2. Six examples showing construction of the regression model relating observed yield (top panels) to the independent variables (middle and bottom panels–only the seasonal average observation is plotted) for example crops and political units (filled black circles) over 35 years 1974 to 2008. The models were then used for predicting yields for historical (filled blue circles) and current (filled red circles) climate conditions with time terms switched off as we are only interested in the difference to yield from difference in climate. Out-of-sample predictions do not occur as the historical and current conditions are bounded within the training weather conditions (middle and bottom panels for seasonal conditions). Yield predictions for individual years 2009 to 2013 (out of sample) are shown in open green circles and observed yield in open black circles with 5 year average error reported as follows for the specific political unit (noted in the figure) in a country: (a) barley (France) -9.3% error, (b) cassava (Brazil) 2.4% error, (c) maize (USA) -13.3% error, (d) soybean (USA) -5.3% error, (e) rice (Bangladesh) 18.4% error, and (f) wheat (China) -6.2% error. https://doi.org/10.1371/journal.pone.0217148.g002 The global summary of the mean climate changes impact on crops contained in the AR5 [2] report uses country-, regional-, and farm-level studies. These studies varied across methods, scales, and time periods and were for the top four global crops–maize, rice, wheat and soybean. Consequently a range of impacts was found; the rice and soybean 10th to 90th percentiles estimates covered both positive (benefit) and negative (reduction) yield impacts from climate trends. Similar multi-method analyses provide conflicting results across regions and crops [38] though recent analysis in wheat that compared point and grid based crop model simulations, with statistical regression approach, were consistent [11]. We tracked ~20,000 political units globally for 10 crops, providing more detail on the spatial resolution and a larger number of crops than previously studied [3, 11, 14, 19]. Further, this is the first observational global study reporting the impact of current climate change on the yields of the top ten global crops (Fig 1, Table 1), and six of them, barley, cassava, oil palm, rapeseed, sorghum and sugarcane are the next important global crops after maize, rice, wheat, and soybean for daily dietary calories. Some such as cassava and sorghum are a major source of calories in food insecure regions (~3.5% and ~3.4% of total food calories provided respectively). Sugarcane, oil palm and rapeseed are important commodity crops with the latter two having greatly expanded in recent decades. Another implicit advancement in our report is from computing the harvested area weighted weather statistics from reported harvesting information each year at the unit of the study ~20,000 political units. We find a range of impacts of mean climate change on crop yields (Fig 1, Table 1 and S4 Table) and production in different regions. We found that crop yields across Europe, Sub-Saharan Africa and Australia had in general decreased because of climate change, though exceptions are present. Similar variations are seen in other crops and regions all over the world. They are indicative of the underlying variations of agronomic growing conditions that range from the agro-meteorological to crop management. Precipitation variations for example are much heterogeneous in their trends (S3 and S4 Figs) and thus captured well using time series analysis per political unit compared to the more homogeneous patterns of temperature trends (S1 and S2 Figs). Our analysis focuses on historical precipitation and temperature change impacts on crop production and food security. Future studies should explore impacts from extreme temperature changes (for example determining thresholds [39, 40] and exposure to killing degree days [41]), extreme precipitation impacts [42] both historically, as well as for the expected future warming [43] and intensification of the hydrological cycle for larger number of crops as well as for livestock. Various agronomic changes and advancements can also mask climate change adaptation measures: irrigation expansion can occur for stabilizing crop production, expansion into new areas, or for double cropping in the dry season; but they can also serve as a measure to counter the effects of extreme heat [44]. Crop planting dates can advance in colder regions of the world as agronomic techniques advance but also because climate is getting warmer on average [45]. These nuances should be explicitly modeled in future studies to tease out such individual contributions. Other related variables that could be considered in future global studies are radiation use efficiency changes [46], weed, pest and pathogen infestation changes [47], soil moisture variability on crop productivity [48–49] and other important biophysical changes to determine their relative contributions to crop yields at the local to the global scale. Changes in both current and future variability in climate [43] and changes in future extreme events [39–42] need further analysis as well as the changing strategy of farmers especially in globally traded commodity crops. Lastly, crop yields and production are not only impacted from climate change, but also drive climate change [33]. Our findings here illustrate that climate change has potentially already affected global production of the ten largest crops and the production of consumable food calories in specific countries and globally. The approach used here complements the long-range projections, avoids many of the challenges faced by process-based models, and analyzes a broader set of crops. Although recent climate change has likely reduced overall consumable food calories in these ten crops by ~1% (or ~0.5% across all food calories), there is much variability among crops and regions. These findings can be used to target interventions for adaptation to climate change through better management, crop breeding, and switching crops as climate continues to evolve.

Methods Data We studied the top ten global crops from around the world where they are commonly harvested. Not all crops are harvested everywhere and each year. To conduct the study we used two datasets 1) climate and weather, and 2) crop yields and harvested areas. Climate and weather variables were all derived from the Climate Research Unit (CRU) TS4.01 [18] dataset and its previous data versions have been commonly used in previous studies [13, 16, 20]. The CRU TS4.01 data is a major upgrade and supersedes all previous versions of the CRU data. We mapped the half-degree CRU derived climate and weather information to the crop harvested areas of each political unit (S1 Text section 1.1). We carried out two major revisions to our crop data that were accessed from public sources [16] and links to the sources are provided in S1 Text section 1. We increased the spatial resolution tracked to ~20,000 political units, and from four major crops being tracked–maize, rice, wheat, and soybean–to tracking the top ten global crops–barley, cassava, maize, oil palm, rapeseed, rice, sorghum, soybean, sugarcane, and wheat. Further details of the crop data and the countries now being tracked at the sub-national levels (one or two levels below the national scale) are provided in S1 Text section 1. The new political units were mapped using information from the public source GADM (https://gadm.org/). Inclusion of sub-national information enables capturing geographical differences, which is missed when studies are restricted to country level political units, especially applicable in cases dealing with geographically large countries. We could not entirely remove this challenge i.e. countries such as Kazakhstan were studied as a single political unit, but in 86 large countries (S1 Text) we conducted the study at one to two administrative divisions down from the country level. To show the effects of studying large countries as a single unit versus sub-national analysis consider the impact of climate change on rice in Cambodia and Malaysia. When analyzed at the country level, rice yield change of -0.03, and +0.15 tons/ha/year was found for Cambodia and Malaysia respectively. Analysis conducted at the subnational scale when summarized to the country level however shows -0.12, and -0.06 tons/ha/year rice yield change for Cambodia and Malaysia respectively. The difference is because when studied at subnational level we are able to isolate the impacts among regions. In the case of Cambodia in the rice bowl central plains region (like Kampong Thom district) and surrounding districts large rice yield losses occur. In the southern districts along the Gulf of Thailand (districts like Koh Kong) climate change benefitted rice yields, but the gains were not enough to overcome the losses elsewhere, leading to overall national level rice yield losses. Similar effect of geography is seen in Malaysia between the eastern (gains in rice yield) and western parts (losses in rice yields) of peninsular Malaysia on the two parts of the natural divide of the Titiwangsa Range; this subnational signal can only be captured through sub-national high-resolution analysis. Thus where possible higher resolution analysis will lead to more precise results. Analysis Our statistical model is a 15-parameter equation, relating crop yields to weather variables at each political unit and of the form: (1) where, t = time (year). The time term is included to account for technological/management changes. We have linear (t) and squared time (t2) terms to account for slow and rapid changes respectively as similarly used in previous studies [16, 19]. α 1 and α 2 are the coefficients associated with the linear and squared time terms respectively. Time terms account for omitted variable bias. sP and sT are terms to account for the contribution of gradual/linear fluctuations in the main crop growing season average precipitation (temperature) conditions on crop yield; α 3 and α 4 are respectively the coefficients associated with the seasonal precipitation and temperature. P and T respectively represent precipitation and temperature. sP2 and sT2 are terms that account for the contribution of extreme/quadratic seasonal precipitation and temperature conditions respectively to crop yields. α 5 and α 6 are the coefficients associated with these two terms. The contribution of the interaction between linear and extreme seasonal precipitation and temperature to crop yield at the political unit is from the sP*sT and sP2*sT2 terms respectively; α 7 and α 8 are the coefficients associated with these two terms. aP and aT terms are similar to the seasonal sP and sT terms respectively, but for capturing the contribution of annual (one year prior to harvest) weather conditions to yields; α 9 and α 10 are the coefficients associated with these two terms. Similarly we have terms aP2 and aT2 to account for contribution of annual extreme weather to yields, and annual interaction terms aP*aT and aP2*aT2. Coefficients associated with them are α 11 , α 12 , α 13 and α 14 . Annual terms account for antecedent weather conditions, secondary, and third season crops, and staggered crop production, and is similar to previous modeling setup [16]. The constant of the regression is k and ε is the error term. Each political unit had this form of the model. Setting up our regression as time series analysis allows us to more accurately capture the effect of precipitation variations. This is because precipitation variations are highly heterogeneous and previous attempts that used panel regressions are less sensitive to precipitation variations [50]. The setting up of this form of the model is uniquely to answer only questions on climate trends/mean climate change impact on crop yields for historical and current conditions. To answer other questions such as identifying temperature thresholds beyond which crops have severe yield losses other forms of model specifications are required, such as an eighth-order polynomial function of temperature [39]. A linear fit using observed weather variables (Fig 2) with observed yield (Eq 1 above) was thus constructed for each crop and political unit (the fitting method is QR decomposition). The model was next tested for significance at the p < 0.05 levels (ANOVA F-statistics; NULL model is the average yield). If the linearly fitted model was significant it was then used to conduct predictions of yield for historical and current climate condition (which are in-sample predictions, as historical and current climate condition are within the range of observed weather–Fig 2). The difference in the predicted yield is the likely observed impact on crops yield due to mean climate change/climate trends at the level of the political unit (see more elaboration in S1 Text sections 2). Residuals were found to be normally distributed and not autocorrelated (S7 Fig), with overall R2 of the models generally greater than 0.8 (S3 Table, S8 Fig), indicating that the models could be used to answer questions on climate change. The lowest (minimum) R2 value in the linear regressions retained for this analysis from anywhere globally was in wheat (0.42 for the political unit of Ryazan oblast in the Central Federal region of the Russian Federation; the root-mean-squared error here was 0.4 tons/ha whereas the observed average yields are ~2.7 tons/ha/year). In all these cases, the goodness of fit test rejects the hypothesis of a lack of fit to the data (at the p = 0.05 level). Cross validation of model yield predictions against observations (2009 to 2013) showed low prediction errors (S2 Table). We do not report the Akaike or Bayesian Information criteria (AIC/BIC) here since these merely reinforce the goodness-of-fit results. Coefficient confidence intervals may be obtained by routine methods, and are not reported here for brevity, but are available. Low model R2 that exists in some regions for some crops such as in Ryazan oblast for wheat should be borne in mind when interpreting results for those areas. Consistent with this approach, we restrict use of our models to climate conditions within the bounds of the observed weather conditions under which the models were constructed i.e. we conducted within-sample predictions and avoided out-of-sample predictions (see S1 Text section 2 and S5 and S6 Figs). Predictions being from statistical analysis, these are not deterministic, and the globally harvested area averaged model Mean Squared Errors (MSE) are provided in S6 Table and mapped in S11 Fig for each crop. Sensitivity to individual model coefficients (sensitivity = coefficient/change in variable) connected to seasonal temperature, squared seasonal temperature, precipitation and squared seasonal precipitation are provided in S12 Fig for maize, rice and wheat that shows how differently each term behaves; the temperature only and precipitation only change impact is provided in S9 and S10 Figs and the full impact in Fig 1. The model formulation per crop and political unit that we determined are provided in S7 Table and the corresponding maize, rice, soybean and wheat yield data for the study period is provided in S8 Table. We have used a relatively simple statistical framework to tease out the main effects of changed climatic conditions on crop yields. A potential limitation of the analysis presented is that regression models were separately estimated within each geographical unit, thus we ignore spatial auto-correlations. Given our spatial scales–even though finer scaled compared to the country scale–each political unit nevertheless are hundreds to thousands of squared kilometers spatially–diffusion of agronomic knowledge and other socio-economic variables i.e. spatial autocorrelation effects would be relatively limited compared to analyses at much more finer spatial scales such as at village levels. Future studies should explore spatial and temporal multi-scale effects [50]; we do not consider potential physical interaction between the explanatory variables at different scales in different geographical units. We would also like to caveat against interpreting our regression-based analysis as causality. Needless to say, detailed studies using technically more complex statistical models, including causal models, and more extensive model diagnostics, is needed. Ignoring spatial autocorrelation does not necessarily bias the final results, however inclusion of such effects may improve the precision of the estimates. On the other hand, spatio-temporal statistical models can be computationally and mathematically very challenging. Future publications should try to address spatio-temporal dependency, diagnostics associated with statistical modeling and other issues mentioned above in global-scale studies. For computing production change we multiplied the climate change driven yield change with the 5-year average harvested areas (2004–2008); for computing percentage yield change with respect to recent yields we similarly computed a 5-year average yield (2004–2008) at each political unit. Country- and global-scale numbers are similarly all harvested area weighted. Finally we determined the change in food calories in a country as the sum of all of individual crop impacts out of the maximum of ten crops studied; most countries did not harvest all ten crops. Details are provided in S1 Text. Mapping After determining the yield change at the political unit due to climate change we mapped the result back at the political unit level. Each political unit contains grid cells that harvested a crop and we assign all these grid cells the same yield change impact. The results are thus valid only at the scale of the political unit for which data were available. In cases where sub-national data were not available there is a marked difference at the state or country borders. For example, Kazakhstan was studied at the country level only and therefore results are valid only at that scale and the contrast with for example the Russian Federation boundary are from mapping. In regions of the world where higher resolution information was available boundary effects are not present e.g. US-Mexico boundary for maize, India-Nepal-Bangladesh boundary for wheat or rice.

Acknowledgments The Global Landscapes Initiative team and others at the Institute on the Environment provided valuable feedback throughout the research and writing process. DKR, PCW and JSG were partially supported by funding from the Gordon and Betty Moore Foundation and the Belmont Forum/FACCE-JPI funded DEVIL project (NE/M021327/1), and the Institute on the Environment. MC was supported by the Wellcome Trust, Our Planet Our Health (Livestock, Environment and People—LEAP), award number 205212/Z/16/Z. SC was partially supported by the National Science Foundation (NSF) grants DMS-1622483 and DMS-1737918. The funders had no role in study design, data collection or analysis, decision to publish, or preparation of the manuscript. We thank the three anonymous reviewers whose constructive comments and suggestions led to further improvement in the manuscript.