# Assessing the capacity of agricultural research and development to increase the stability of global crop yields under climate change

Toshichika Iizumi, Toru Sakai, Yoshimitsu Masaki, Kei Oyoshi, Takahiro Takimoto, Hideo Shiogama, Yukiko Imada, David Makowski

► To cite this version:

Toshichika Iizumi, Toru Sakai, Yoshimitsu Masaki, Kei Oyoshi, Takahiro Takimoto, et al.. Assessing the capacity of agricultural research and development to increase the stability of global crop yields under climate change. PNAS Nexus, 2025, 4 (4), 10.1093/pnasexus/pgaf099 . hal-05203558

HAL Id: hal-05203558

https://hal.inrae.fr/hal-05203558v1

Submitted on 19 Dec 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Assessing the capacity of agricultural research and development to increase the stability of global crop yields under climate change

Toshichika Iizumi IDa,*, Toru Sakai IDb, Yoshimitsu Masaki IDa, Kei Oyoshic, Takahiro Takimoto IDa, Hideo Shiogama Id, Yukiko Imadae and David Makowski IDf

Edited By Joann Whalen

## Abstract

Agricultural research and development (R&D) has increased crop yields, but little is known about its ability to increase yield stability in the context of increasingly frequent extreme weather events. Using a grid yield dataset, we show that from 2000 to 2019, the SD of yield anomalies for maize, rice, wheat, and soybean increased in 20% of the global harvested area. Based on random forest models relating yield anomaly to climate, soil, management, and public R&D expenditure, we show that cumulative agricultural R&D expenditure, proportion of growing season exposed to optimal hourly temperatures, and dry and very wet days are key factors explaining crop yield variability. An attribution analysis based on large ensemble climate simulations with and without human influence on the global climate shows that unfavorable agroclimatic conditions due to climate change has increased SD, while higher R&D expenditure has led to more contrasting trends in SD over 2000–2019. Although R&D has continued steadily in most countries, this study indicates that the progress made in R&D since 2000 may have lagged behind the unfavorable effect of climate change on yield variability.

## Significance Statement

Improved technologies and management practices—fostered by agricultural research and development (R&D)—have helped to increase yield stability, but it is not known whether R&D investments have been sufficient to fully offset the increase in yield variability induced by climate change. Based on a global dataset, our analysis quantifies the relationships between R&D expenditure, climate change, and yield stability for major agricultural commodities worldwide. Our analysis reveals that R&D investments have probably not been high enough since 2000 to fully mitigate the increased variability of yields resulting from climate change. Future research is needed to establish a causal relationship between R&D and yield stability underlying the empirical results found in this study.

## Introduction

Among the four pillars of food security, namely availability, access, utilization, and stability (1), ensuring the production stability is crucial for preventing any deleterious consequences cascading to agrifood systems downstream. The disruptions to food supply chains caused by recent conflicts and pandemic reaffirm the importance of stability in the interconnected world (2–4). In the past, production stability has been improved by optimizing growing areas, harvestable area fraction, annual number of harvests, and cultivar tolerance to abiotic and biotic stresses (5, 6). Increasing yield stability remains an important avenue to make

crop production more climate-resilient, especially because yield growth has been the predominant contributor to production increase in many parts of the world in recent decades, while the increase in cultivated areas has been a much smaller contributor since 2000 (5, 7).

From 1981 to 2010, yield stability for major crops, namely maize, rice, wheat, and soybean, decreased in 9–22% of the global harvested area in response to increased crop exposure to heat, yet 19–33% of the area saw stabilizing yields (8). After the study by Iizumi and Ramankutty (8) published in 2016, many production shocks associated with extreme climate events have been reported (9–12). Research and development (R&D) on climate-resilient

<figure><img src="image_2.png" /></figure>

<footnote>aInstitute for Agro-Environmental Sciences, National Agriculture and Food Research Organization, 3-1-3 Kannondai, Tsukuba, Ibaraki 305-8604, Japan</footnote>

<footnote>bJapan International Research Center for Agricultural Sciences, 1-1 Ohwashi, Tsukuba, Ibaraki 3058-686, Japan</footnote>

<footnote>cEarth Observation Research Center, Japan Aerospace Exploration Agency, 2-1-1 Sengen, Tsukuba, Ibaraki 305-8505, Japan</footnote>

<footnote>dNational Institute for Environmental Studies, 19-2 Onogawa, Tsukuba, Ibaraki 305-8506, Japan</footnote>

<footnote>eAtmosphere and Ocean Research Institute, The University of Tokyo, 5-1-5 Kashiwanoha, Kashiwa, Chiba 277-8564, Japan</footnote>

<footnote>fUniversité Paris-Saclay, INRAE, AgroParisTech, UMR 518 MIA, Palaiseau, France</footnote>

<footnote>*To whom correspondence should be addressed: Email: iizumit@affrc.go.jp</footnote>
agrifood systems is expanding, including breeding targeting emerging agroclimatic conditions (13) and climate-smart agriculture. However, although R&D innovations have been extensively tested at experimental stations, there is limited evidence showing that advances in R&D and producers' adoption of improved seeds, technologies, and management practices have reduced climate risks at large spatial scales (14). Despite improvements in breeding and weather forecasting in the last decades, some studies suggest that the yield sensitivity to extreme climate events has not been substantially reduced in high-yielding regions (15, 16). It is therefore crucial to update our understanding on yield stability and its drivers. In this paper, we contribute to this effort by studying how climate change and agricultural R&D have affected yield stability since 2000.

Recent advances in climate attribution science and resulting climate model simulations (17–19) allow researchers to estimate the contributions of human-induced climate change to average yields (7, 20) and crop failures (21, 22). However, contrasting to the accumulated evidence showing that climate change has slowed down the growth of average yields (7), little is known about climate-induced changes in yield stability. In addition, it remains difficult to assess the capacity of R&D on reducing climate risks, although this information is often requested by research project funders (23). To fill these gaps, we present here a global analysis of the changes in yield stability due to climate change and R&D from 2000 to 2019, using a large ensemble climate model simulation (17) and machine learning models that relate interannual variation in yields of maize, rice, wheat, and soybean at the 0.5° grid scale (~55 km at the equator) to climate, soil, management, and agricultural R&D expenditure. The four crops considered here are of high importance as they provide nearly two-thirds of the world agricultural calories (24).

# Results

## R&D capacity to increase and stabilize crop yields

Yield anomaly—the departure of annual yield from trend of average yield—is a key metric measuring year-to-year yield fluctuation. The analysis of country-level yield anomalies over 2000–2019 showed that increases in public R&D expenditure are proportional to increases in average yields. Considering the four crops taken together, countries with large areas under cultivation (≥1 Mha) had invested considerably more in R&D than countries

with smaller cropping areas (Fig. 1A). Average yield and yield anomaly SD are positively correlated (Fig. 1B), but, for a given average yield level, countries with higher levels of R&D expenditure often achieved higher stability (lower SD and/or lower coefficient of variation [CV] = SD/average) than countries with lower levels of R&D expenditure (Fig. 1C).

## Recent patterns of agricultural R&D expenditure

The level of cumulative agricultural R&D expenditure varied considerably between countries, although increasing trends are quite common (Fig. 2A). Compared with annual R&D expenditure, their cumulative values, which take into account the time lag between research and adoption of new technologies by producers and technology obsolescence, are relevant to explain the increase in average yield and yield stability. Since 2000, annual R&D expenditure was the highest in China, followed by the United States of America (Fig. 2B). China's annual spending increased rapidly from 3.6 billion USD ($B) in 2000 to 22.9$B in 2019, while the US annual spending increased from 4.2$B to 5.2$B for the period. The level of annual R&D expenditure in other countries was much less than these two countries. Although cumulative R&D expenditure increased in most countries, for some countries, cumulative R&D expenditure stagnated (e.g. in France) (Fig. 2A and C). When the level of annual R&D expenditure was not sufficient to maintain the level of cumulative R&D expenditure, which is subject to annual rates of technology obsolescence, cumulative R&D expenditure declined with time (e.g. in Russia and Sudan).

## Recent changes in yield stability

To analyze the evolution of yield stability, SD was selected as the main metric because SD was more sensitive than CV in detecting climate-induced changes in yield anomalies. Under climate change, yield anomaly SD for a certain period may increase as a result of both occurrences of poor harvests (due to severe agroclimatic conditions) and good harvests (due to favorable agroclimatic conditions, or to elevated carbon dioxide concentrations, which enhance photosynthesis). However, with CV, such a change in SD can be masked by an increase in average yield because CV decreases when average yield increases faster than yield anomaly SD over time.

Our grid-wise analysis showed that 20% of the global harvested area (129.2 Mha, out of 646.2 Mha in 2010) had increased yield anomaly SD in 2000–2019 (Fig. 3B), while areas with increased

<figure><img src="image_3.png" /><figcaption>**Fig. 1.** Relationship between R&D, yield increase and stability. A) Natural log of the cumulative annual public R&D expenditure for agriculture since 1995 vs. average yield. B) Average yield versus the yield anomaly SD. C) Average yield versus the yield anomaly CV (CV = SD/average). These results were obtained over 2000–2019 for 115 crop-producing countries. The selected countries have an annual production of 0.1 Mt or more for at least one of the crops considered here (maize, rice, wheat, and soybean). The cumulative R&D expenditure is in units of constant 2015 million USD ($M, log-transformed). The size and color of the symbols indicate the total area harvested and accumulated agricultural R&D expenditure, respectively. The solid lines show the best-fitted nonparametric quantile regression for the median that divide data into two equal portions. For A to C), all slope values are significant (*P* < 0.01; *F* test).</figcaption></figure>
<figure><img src="image_4.png" /><figcaption>**Fig. 2.** Patterns of agricultural R&D expenditure. A) Global map of trends in cumulative R&D expenditure in 2000–2019. B) Country annual R&D expenditure in constant 2015 USD. C) Country cumulative R&D expenditure since 1995. Classes of increasing (decreasing) cumulative R&D expenditure are classified by growth rate, i.e. “rapidly” and “slowly,” respectively, if a given annual rate is faster and equal to or slower than the median value. The median annual rate is 4.8% for the increasing trend and –1.5% for the decreasing trend; both are relative to the 2015 level. The colors in the time series B, C) correspond to those in the map A).</figcaption></figure>

CV were only 58% of the areas with increased SD (75.8 Mha; Fig. 3A). The overall increasing trends in average yields over the same period explained this difference in results between SD and CV (Fig. 3C).

## Modeling the effects of climate and R&D on yield anomaly

Machine learning random forest (RF) models were trained with current climate conditions (represented by the retrospective meteorological forcing data; see Materials and methods) to attribute the changes in modeled yield anomaly SD to the R&D and climate-related explanatory variables (Table 1). The models explained 19% (rice) to 43% (soybean) of the variance of the grid yield anomaly between the locations, seasons, and years when evaluated using out-of-bag data. Intermediate percentage of explained variance was obtained for maize and winter wheat (26% for both) and spring wheat (39%). Concerning the model ability to reproduce interannual yield variations, the leave-one-year-out cross-validation led to average correlations ranging from 0.639 (major-season rice) to 0.793 (soybean) and to average root mean squared errors ranging from 0.66 (spring wheat) to 0.81 (major-season rice) in SD units (Figs. S4 and S5).

Importantly, the RF model projections obtained with d4PDF factual climate simulations (representation of the actual climate simulated by the climate model, including human-induced climate change) were able to reproduce the geographic pattern of the changes in yield anomaly SD for the crops, with moderate agreement with the actual pattern derived from the grid yield dataset (Cohen's kappa [$\kappa$] of 0.55; Fig. 4A and C). The agreement with the actual pattern was as high as $\kappa=0.71$ when the RF models were forced by the retrospective meteorological forcing data (Fig. 4A and B). We also checked that d4PDF factual and counterfactual climate simulations fell within the ranges covered by the training dataset as the RF model performance often diminish outside the training data (only 0.18% of cases across all climate variables, locations, seasons, years, and ensemble members were outside the ranges covered by the training dataset).

## Variable importance

The RF models were used to rank the relative importance of R&D, climate, soil, and management variables in explaining the yield anomaly. Two methods were considered to determine the importance of variables (see Materials and methods). Cumulative R&D expenditure was identified as the primary factor, irrespective of

the crops and methods (Fig. 5). N application rate was often detected as the next important factor after R&D expenditure for the crops except soybean (Fig. S6). Temperatures ranging from 15 to 35 °C, which encompass the exposures to temperature fluctuations around optimal thermal conditions for crops, were identified as important, albeit with some crop-specific variation. For maize, rice, and soybean, the relatively higher temperatures from the fifth bin (15–20 °C) to the ninth bin (≥35 °C) were more important. For winter and spring wheat, the relatively lower temperatures from the fourth bin (10–15 °C) to the fifth bin (15–20 °C) showed stronger effects. Among the precipitation variables, the first bin (<1 mm day⁻¹) and the fifth bin (≥30 mm day⁻¹) tended to the most important, with variations between the crops and methods. On average, soil organic carbon (SOC) was one of the most important soil variables. Irrigation intensity was less influential.

## Climate effect

Effects of climate and R&D on the changes in yield stability at SD level were estimated by comparing the yield anomalies simulated by the RF models for different scenarios. Two yield anomaly time series were generated using d4PDF factual and counterfactual climate simulations, representing the historical climate (including human-induced climate change) and the preindustrial climate (without human-induced climate change), respectively (fc and ct runs; Fig. 6). R&D expenditures were fixed to those observed in both scenarios. We classified the change in agroclimatic condition due to climate change as “unfavorable” (“favorable”) when the simulated yield anomaly SD was higher (lower) with the factual climate than with the counterfactual climate.

The global grid-wise attribution revealed that unfavorable agroclimatic conditions were relatively more prevalent in the cropland area where an increase in yield anomaly SD occurred than in the area with decreasing yield anomaly SD. The area with unfavorable agroclimatic conditions was estimated to account for the half (0.5) of the area with an increase in yield anomaly SD; this proportion (0.5) was derived by comparing the area with both unfavorable agroclimatic conditions and increasing SD (10% of the global harvested area; Fig. 7A and D) and the area with increasing yield anomaly SD all agroclimatic conditions taken together (20% of the global harvested area; Fig. 7C and D). This proportion (0.5) is twice as high as the proportion (0.25) of area with unfavorable agroclimatic conditions and decreasing SD (4%) in the area with decreasing yield anomaly SD all agroclimatic conditions taken together (16%) (Fig. 7D). In other words,
<figure><img src="image_5.png" /><figcaption>B Change in yield anomaly SD</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>C Change in average yield</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>D Number of crops</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>**Fig. 3.** Relative counts of crops with specific pattern of yield change. A) Increase in yield anomaly CV. B) Increase in yield anomaly SD. C) Increase in average yield. D) Number of crops considered here that are harvested in 2010 (the midpoint of the study period). Statistical significances of monotonic trends in 2000–2019 are examined using the two-sided Mann–Kendall test. Maize, rice, wheat, and soybean are considered together. The results for each crop are available in Figs. S1–S3.</figcaption></figure>

the frequency of unfavorable agroclimatic conditions is twice as high in areas where the yield anomaly SD is increasing than in areas where the yield anomaly SD is decreasing.

Consistent with the above result, favorable agroclimatic conditions were relatively more common in the area where a decrease

in yield anomaly SD occurred than in the area with increasing yield anomaly SD. Favorable agroclimatic conditions were found in a quarter (0.25) of the area with decreasing yield anomaly SD; this proportion (0.25) compared the area with both favorable agroclimatic conditions and decreasing SD (4%) and the area with decreasing yield anomaly SD all agroclimatic conditions taken together (16%) (Fig. 7D). This proportion (0.25) is 1.7 times higher than the proportion (0.15) between the area with both favorable agroclimatic conditions and increasing yield anomaly SD (3%) and the area with increasing yield anomaly SD all agroclimatic conditions taken together (20%).

## R&D effect

The R&D effect was estimated by comparing two series of yield anomaly simulations, i.e. a series simulated with the actual R&D expenditures in 2000–2019 and another series simulated without R&D investment (fc and fc.em runs; Fig. 6). The results showed that a significant yield-stabilizing effect of R&D was found in a large part (0.44) of the cropland area where a decrease in yield anomaly SD occurred. This proportion (0.44) was derived by comparing the area with a significant R&D effect of 7% and the area with decreasing yield anomaly SD of 16%, both relative to the global harvested area (Fig. 7B, C, and E). This proportion (0.44) is 4.4 times greater than the proportion (0.10) calculated between the area with a significant R&D effect (2%) and the area with increasing yield anomaly SD (20%).

On the contrary, the yield-stabilizing effect of R&D was not found to be significant in a dominant proportion (0.85) of the area where an increase in yield anomaly SD occurred. This value (0.85) was derived from the area with a nonsignificant R&D effect of 17% and the area with increasing yield anomaly SD of 20% (Fig. 7E), which is 1.7 times greater than the value (0.50) calculated from the area with a nonsignificant R&D effect (8%) and the area with decreasing yield anomaly SD (16%). Importantly, the area with a significant R&D effect became smaller when estimated over a more recent, shorter period than over the entire study period (Fig. 7F).

## Sensitivity to different R&D assumptions

The estimated area shares for the climate and R&D effect categories were found to be robust to the different assumptions used in the R&D variable calculation. We took into account whether the cumulative R&D expenditure was for agriculture or for the four crops, whether the time lag between research and adoption of technologies by producers was relatively shorter (6 years) or longer (12 years), and whether the technology obsolescence was relatively slower (10% per year) or faster (20% per year) (see Materials and methods). The main results described in the former subsections (specifically, the fact that R&D is the primary explanatory variable explaining yield anomaly, and unfavorable agroclimatic conditions and/or nonsignificant R&D effect are proportionately more prevalent in the cropland area experiencing a decrease in yield stability [higher SD]) were robust to the different R&D assumptions (Figs. S9 and S10).

## Region- and crop-specific characteristics

We identified the regions and crop species for which R&D were able to compensate for the negative impacts of climate change on yield stability after 2000. For the four crops taken together (479.8 Mha), R&D had no significant effect in 75% of this area. However, R&D was found to have a significant stabilizing effect on large shares of croplands in several specific crop species and
<table><caption>**Table 1.** Variables used in the RF models.</caption><thead><tr><th>Category</th><th>Symbol</th><th>Description</th><th>Unit</th></tr></thead><tbody><tr><td>Crop yield</td><td>Y</td><td>Z-scored yield anomaly</td><td>SD</td></tr><tr><td rowspan="10">Climate</td><td>t1</td><td>First hourly temperature bin (hourly temperature [T] &lt; 0 °C)</td><td>Fraction<sup>a</sup></td></tr><tr><td>t2</td><td>Second hourly temperature bin (0 ≤ T &lt; 5)</td><td>Fraction</td></tr><tr><td>t3</td><td>Third hourly temperature bin (5 ≤ T &lt; 10)</td><td>Fraction</td></tr><tr><td>t4</td><td>Fourth hourly temperature bin (10 ≤ T &lt; 15)</td><td>Fraction</td></tr><tr><td>t5</td><td>Fifth hourly temperature bin (15 ≤ T &lt; 20)</td><td>Fraction</td></tr><tr><td>t6</td><td>Sixth hourly temperature bin (20 ≤ T &lt; 25)</td><td>Fraction</td></tr><tr><td>t7</td><td>Seventh hourly temperature bin (25 ≤ T &lt; 30)</td><td>Fraction</td></tr><tr><td>t8</td><td>Eighth hourly temperature bin (30 ≤ T &lt; 35)</td><td>Fraction</td></tr><tr><td>t9</td><td>Ninth hourly temperature bin (35 ≤ T)</td><td>Fraction</td></tr><tr><td>p1</td><td>First daily precipitation bin (daily precipitation [P] &lt; 1 mm)</td><td>Fraction<sup>b</sup></td></tr><tr><td>p2</td><td>Second daily precipitation bin (1 ≤ P &lt; 5)</td><td>Fraction</td></tr><tr><td>p3</td><td>Third daily precipitation bin (5 ≤ P &lt; 20)</td><td>Fraction</td></tr><tr><td>p4</td><td>Fourth daily precipitation bin (20 ≤ P &lt; 30)</td><td>Fraction</td></tr><tr><td>p5</td><td>Fifth daily precipitation bin (30 ≤ P)</td><td>Fraction</td></tr><tr><td rowspan="10">Soil</td><td>soc</td><td>Topsoil organic carbon content</td><td>kg C m<sup>-2</sup></td></tr><tr><td>wsc</td><td>Available water storage capacity of the topsoil</td><td>mm H<sub>2</sub>O per m of the soil unit</td></tr><tr><td>clay</td><td>Topsoil clay fraction</td><td>% weight</td></tr><tr><td>silt</td><td>Topsoil silt fraction</td><td>% weight</td></tr><tr><td>sand</td><td>Topsoil sand fraction</td><td>% weight</td></tr><tr><td>bd</td><td>Topsoil bulk density</td><td>kg dm<sup>-3</sup></td></tr><tr><td>cec</td><td>Cation exchange capacity of the clay fraction in the topsoil</td><td>cmol per kg</td></tr><tr><td>ph</td><td>Topsoil pH in water</td><td>-log(H<sup>+</sup>)</td></tr><tr><td>irr</td><td>Irrigation intensity (the irrigated area divided by harvested area)</td><td>Fraction</td></tr><tr><td>napp</td><td>Nitrogen application rate</td><td>kg N ha<sup>-1</sup> year<sup>-1</sup></td></tr><tr><td>R&D</td><td>agks</td><td>Knowledge stock (log natural of cumulative public R&D expenditure on agriculture since 1995). See Materials and methods for variants calculated using different assumptions</td><td>log($)</td></tr></tbody></table>

regions, especially for soybean (42% [significant R&D effect] vs. 58% [nonsignificant]), in sub-Saharan Africa (41 vs. 59%), East Asia (35 vs. 65%), and North America (29 vs. 71%) (Fig. 8A).

Climate change had no significant effect on yield stability in almost half of this area (43%), but led to unfavorable agroclimatic conditions (higher SD) and favorable agroclimatic conditions (lower SD) in 36 and 21% of the area, respectively. The negative effect of climate change on yield stability was particularly pronounced for maize (45% [unfavorable] vs. 22% [favorable]) and soybean (57 vs. 23%), and in North America (81 vs. 3%), followed by East Asia (33 vs. 18%). On the contrary, the shift to favorable agroclimatic conditions was relatively prominent in South Asia (10 vs. 36%) (Fig. 8B).

## Discussion

Yield stability has increased in 16% of the global cropland. Although not mutually exclusive, nearly half of this area is associated with a significant stabilizing effect of R&D, and a quarter of this area is associated with a shift to favorable agroclimatic conditions. In contrast, yield stability has decreased in 20% of the global cropland. Half of this area is associated with a shift to unfavorable agroclimatic conditions, while about nine-tenths of this area is associated with a lack of significant R&D effect, although these areas may overlap. These results provide insights into the association between investment in agricultural R&D, climate change, and resulting changes in yield stability. In addition, based on these results, we could formulate the hypothesis that, in the global cropland area where yield stability has decreased, the progress made in R&D since 2000 has lagged behind the negative impacts of climate change. This question needs to be addressed in future research.

Although R&D effect is relatively more prevalent in the area with increasing yield stability, significant R&D effect is also

detected in some part of the global cropland with unchanged and decreasing yield stability (8 and 3%, respectively). In this study, the R&D effect is considered significant when R&D expenditures have contributed to a decrease of yield anomaly SD over the study period. Even if the R&D effect is found to be significant, decreasing or unchanged yield stability may occur if the R&D effect is not strong enough to reverse the trend in yield anomaly SD. However, without any R&D effect, the yield stability status in these areas would probably deteriorate further.

It should be noted that the R&D effect on yield stability identified here is empirical rather than causal, as some confounding factors potentially correlated with R&D expenditure were not considered in this study. Process-resolving analyses would be useful to determine the causal effects of individual technologies and practices and confirm our empirical evidence of the yield-stabilizing effect of R&D, which could ultimately support the implementation of effective interventions. Although applied to study average yield and not yield stability, a relevant example is presented in Ref. (26), which decomposes the historical yield increases into individual drivers, including increases in the use of land, water, and other inputs, and increases in the efficiency of input use, with an analysis of policy implications.

## Important factors of yield anomaly

The yield-stabilizing effect of R&D is not totally surprising. For instance, new varieties bred through R&D generally have greater tolerance to suboptimal growing conditions and higher yield potential than earlier varieties (27–29). In addition to breeding, R&D informs producers about practices to reduce climate risks, such as sowing date, soil, water, and fertilizer management, and pests, diseases, and weed control.

The temperatures in the growing season, as well as occurrences of dry and heavy rainfall days, are identified to be

<footnote><sup>a</sup>Fraction of exposure hours to total hours from planting to harvesting. The same applies to other temperature bins. <sup>b</sup>Fraction of exposure days to total days from planting to harvesting. The same applies to other precipitation bins.</footnote>
<figure><img src="image_9.png" /></figure>

<figure><img src="image_10.png" /></figure>

<figure><img src="image_11.png" /><figcaption>**Fig. 4.** Actual and reproduced changes of yield anomaly SD. A) Actual pattern (same as in Fig. 3B). B) Reproduction obtained with the trained RF models using the retrospective meteorological forcing data. C) Reproduction using d4PDF factual climate simulation. In C), a cell is classified as "increasing" if at least one-third or more of the 100 climate simulation members have a significant increasing trend in the factual yield anomaly SD series. Cohen's kappa statistics (κ) calculated for each of B) and C) against A) are shown as an indication of the similarity of the geographic pattern between a given pair.</figcaption></figure>

important climatic factors explaining yield anomaly. This result is consistent with previous findings (30–32). Under warming, a decrease in the frequency of near-optimal temperatures is accompanied by an increase in the frequency of supraoptimal temperatures, which decreases yield stability (33). The effects on yield stability from heavily wet conditions identified here are more uncertain than that the effects of dry conditions, although there is a growing body of literature reporting the negative impacts of waterlogging and resulting excessive soil water (34).

The importance of N application in explaining yield anomaly can be explained by previous findings that abiotic stresses due to growing season climate become increasingly important stressors as N inputs increase (35, 36). The importance of SOC found in this study confirms previous findings that drought damage is mitigated by SOC build-up (37). Despite the initial expectation that irrigation would increase yield stability, the contribution of irrigation appears to be limited in this study. The facts that the existence of irrigation equipment does not necessarily ensure water availability during droughts (38) and that yield sensitivity at

<figure><img src="image_12.png" /><figcaption>**Fig. 5.** Ranking of variable importance according to the RF models. The bars show the average rank calculated over the crops (the closer the rank is to 1, the more important the variable is in explaining yield anomaly). The thin lines indicate the SD of the rank between the crops. Two variable importance methods (permutation and node purity) are used to account for uncertainty. The results of individual crops are available in Fig. S6.</figcaption></figure>

regional levels to heat no longer decreases as irrigated area saturates (39) partly explain the counterintuitive outcome.

## Limitations and uncertainties

There are a few limitations to this study. First, we only consider public expenditures and do not include expenditures by private sector and international organizations. Therefore, the R&D effect estimated in this study may be underestimated to some extent, given that the private sector's share of total agricultural R&D expenditure can be substantial, particularly in high-income countries (30%) compared with low- and middle-income countries (5%) (40). Moreover, the role of international organizations for agricultural R&D is crucial in low-income countries, and investments in agricultural R&D for one country may generate spillover benefits for neighboring countries and regions with similar agroclimatic conditions. Data on private R&D expenditure are difficult to obtain. Even for public expenditure, the information on the technology on which R&D expenditure was spent is rarely available (41). Future research should fill these data gaps to distinguish the relative contributions of individual technologies and management practices.

Second, a direct comparison with previous estimates of climate-induced changes in yield stability is not possible due to a lack of relevant literature. There is literature showing changes in average yield (20, 23). However, none of these studies assessed changes in yield stability and the R&D effect.
## Experiment 1. fc run vs. cf run

<figure><img src="image_13.png" /></figure>

## Experiment 2. fc run vs. fc.em run

<figure><img src="image_14.png" /><figcaption>**Fig. 6.** Attribution of changes in yield stability. Model simulations were conducted to estimate the effects of climate and R&D on the changes in yield stability in 2000–2019.</figcaption></figure>

The remaining limitations and uncertainties are mostly related to the underlying data. Third, this study is based on a single grid yield dataset. Different yield datasets often produce different results. Subnational census statistics are expected to be more reliable than satellite- or model-based yield datasets, including one adopted here. However, the spatial resolution of census statistics is coarser than the resolution of d4PDF climate simulation, which compromises the ability to capture the yield impacts from extreme agroclimatic conditions. In addition, some recent grid yield datasets (42) use temperature and precipitation as inputs in their yield estimation procedure, positing the circularity issues because we relate yield anomaly to temperature and precipitation. Fourth, attribution outcomes may be influenced by the choice of meteorological forcing datasets and climate models (43, 44). Fifth, our analyses are based on the harvested area in 2010, although the area for the four crops in 2022 increased 12% since 2010. This increase in area would include not only a response to increasing demand but also a shift of areas suitable for cultivation in response to the observed warming (45, 46). However, we used the grid harvested area map for 2010, the midpoint of the study period, as no annual crop-specific grid harvested area map is available at the global scale. Last, future research should consider more management variables than those considered here, such as pesticide application rates, which are thought to increase yield stability (47).

# Materials and methods

## Experimental design

This section describes the data on crop yields, climate (actual, factual, and counterfactual), and agricultural R&D expenditure, as well as the method used to detect changes in yield stability. Also included are the machine learning models and the yield impact attribution analyses.

## Crop yields

We obtained country annual yield data from the FAO statistical database for maize, rice, wheat, and soybean (48). In addition, we used the 0.5° grid global dataset of historical yields (GDHY) (49). In the GDHY dataset, yield data are available for two seasons for maize and rice (major and secondary) and for wheat (winter and spring). For soybean, data are available for a single main

season. The GDHY dataset is a hybrid of satellite vegetation index and national census statistics reported by FAO. The information on accumulated net primary production over the crop-specific growing period, calculated using remotely sensed leaf area index and fraction of photosynthetically active radiation as well as re-analysis downward shortwave radiation and crop-specific radiation-use efficiency, is used together to spatially disaggregate FAO country average yields to the grid level. Grid yield values are available in 76–92% of the global harvested area with variations by crop and lacking when crop calendar information (50) is missing.

## Climate

We used 0.5° grid daily data on maximum and minimum air temperatures and precipitation, covering three distinct climate conditions in 2000–2019, namely *actual*, *factual*, and *counterfactual*. The *actual* climate refers to observed conditions. The *factual* climate represents a modeled approximation of the actual climate, influenced by both human activities and natural forces (e.g. volcanic eruptions and solar activity). In contrast, the *counterfactual* climate is a modeled preindustrial climate without any appreciable human impacts on the global climate.

The actual climate data were obtained from the global retrospective meteorological forcing data. We used the JRA55-CDFDM-S14FD forcing data (called JCS for simplicity), which is a bias-corrected Japanese 55-year Reanalysis (JRA55) using the cumulative distribution function-based downscaling and bias-corrected method (CDFDM) and the reference climatology from S14FD forcing dataset (51).

For the factual and counterfactual climates, we used a large ensemble of long-term climate simulations with and without historical trends in external forcing to the global climate, called d4PDF (the Database for Policy Decision Making for Future Climate Change) (17, 52). The climate simulations were performed using the Meteorological Research Institute Atmospheric General Circulation Model, version 3.2 (MRI-AGCM3.2), with a grid interval of 60 km (53). Each of the factual and counterfactual climate simulations has 100 ensemble members associated with slightly different initial conditions and small perturbations in sea surface temperatures that represent observational uncertainties. Comparisons between the factual and counterfactual climate simulations with large ensemble members enable researchers to
<figure><img src="image_15.png" /><figcaption>**Fig. 7.** Climate and R&D effects on yield anomaly SD. A) Relative counts of crops showing that climate change in 2000–2019 increased the yield anomaly SD, compared with what would occur under counterfactual climate. B) Relative counts of crops showing that actual R&D had a significant effect in decreasing yield anomaly SD, compared with what would occur in the absence of investments to R&D after 2000. The differences in annual rate of changing in yield anomaly SD in A) and B) are, respectively, tested using the two-sided Wilcoxon rank-sum test (*P* < 0.05). C) Relative counts of crops showing an increase in yield anomaly SD derived using the grid yield dataset (same as in Fig. 3B). D, E) The pie charts show the percentage of the global harvested area in 2010 (25) with color-coded categories for changes in yield anomaly SD as well as the climate and R&D effects. F) Area share for the R&D effect categories estimated from year to year. In F), levels of R&D expenditure, N application, and irrigation were kept unchanged after the year indicated in the x-axis to predict yields with the RF models. The results obtained for individual crops are available in Figs. S7 and S8.</figcaption></figure>

conduct a robust attribution analysis of the effects of human activities on average climate conditions, extreme climate events (54–57), and crop yields (58, 59). The daily outputs from the climate model were interpolated to the 0.5° resolution and bias-corrected using the CDFDM along with the reference climatology from S14FD forcing dataset.

## Agricultural R&D expenditure

Public agricultural R&D includes government, higher education, and nonprofit entities, but excludes the private for-profit sector (40). It includes salaries, operating and program costs, and capital

investments for all entities, excluding the private for-profit sector, involved in agricultural R&D expenditure in a country.

For the Organization for Economic Co-operation and Development (OECD) member economies, we sourced annual data on gross domestic expenditure on R&D expressed as a percentage of gross domestic production (GDP) from the Main Science and Technology Indicators (MSTI) database (60). Agricultural R&D expenditure for OECD member economies was calculated by multiplying the GDP share of domestic R&D expenditure with the agricultural share of GDP and the total GDP, both obtained from the World Bank (61). For 88 low- and middle-income countries,
<figure><img src="image_16.png" /><figcaption>**Fig. 8.** Areas shares of R&D and climate effect categories by crop and region. A) The R&D effect is classified into two categories, i.e. R&D after 2000 has a significant or nonsignificant effect, in increasing yield stability in 2000–2019 under climate change. B) The climate effect is classified into three categories, i.e. climate change, relative to preindustrial levels, leads to a favorable, unfavorable, or no significant change in agroclimatic condition in terms of yield stability. The region codes are as follows: EAs, East Asia; EEuCAs, Eastern Europe and Central Asia; LAmCa, Latin America and Caribbean; MENAf, Middle East and North Africa; NAm, North America; SAs, South Asia; SEAsO, Southeast Asia and Oceania; SSAf, sub-Saharan Africa; and WEu, Western Europe (see also Fig. S11).</figcaption></figure>

data were obtained from the International Food Policy Research Institute (IFPRI)'s Agricultural Science and Technology Indicators (ASTI) database (40). The country annual R&D expenditure data collected were adjusted for inflation to be expressed in 2015 USD.

## Detection of changes in yield stability

Changes in yield stability were detected based on the method of Iizumi and Ramankutty (8). For each grid cell, an annual yield time series (1981–2019) was detrended with a cubic smoothing spline representing time-evolving average yield (Fig. S12A). Although the present study analyzed changes in yield stability between 2000 and 2019, a longer yield time series was used for a more robust identification of yield trends. Yield anomalies (in tonnes per hectare) were computed as the differences between the original yield time series and the fitted trends (Fig. S12B). By applying a 9-year centered moving window, an annual time series of the SD of yield anomalies was derived. A regression line was fitted to the recent portion (2004–2015) of the yield anomaly SD series, and the estimated slope was tested for a monotonic trend using the two-sided Mann–Kendall trend test (62) (Fig. S12C). Note that the yield stability trends analyzed here were based on the yield anomalies over the 2000–2019 period (Fig. S12C). A significant positive (negative) slope indicated a decreasing (increasing) yield stability. Monotonic trends in average yield and the CV time series of yield anomaly in the same period were tested in the same manner as the SD time series.

To aid visualizing multiseason and multicrop results as a single map, the detected changes in yield stability for individual seasons and crops were combined according to the following rules. The major season result was selected as the representative when both major (winter) and secondary (spring) seasons were operated. When only one season was operated, the result for that season was considered even when it was the secondary (spring)

season. Then, the number of crops with specific yield change (either a significant increase in yield anomaly CV, a significant increase in yield anomaly SD, or a significant increase in average yield), out of four crops, namely maize, wheat, rice, and soybean, was counted and divided by the number of crops grown for a given location (Fig. 3D).

## RF yield anomaly models

Machine learning models were developed to relate yield anomaly to climate, soil, management, and R&D expenditure. We used RF regression models (63) and did not explore alternative machine learning approaches as RFs are straightforward to implement and comparable to other approaches in many applications (64, 65). RFs are generally insensitive to correlations between explanatory variables in terms of predictive accuracy and overfitting (64), which further supports their use in this study. A RF algorithm deals with collinearity by randomly sampling a small subset of input variables to build each splitting node of each tree. As a result, the learning algorithm can only choose among a small number of inputs at each splitting node, thus reducing the risk of collinearity.

We calculated the z-score values of yield anomaly by dividing annual yield anomaly by the SD of yield anomaly for the same period as that used to estimate yield trends (Fig. S12D). The z-scored yield anomalies were related to climate, soil, management, and R&D variables (Table 1). Details of the explanatory variables are provided in the subsequent sections.

A specific model was constructed for each crop. Data from major and secondary maize (rice) seasons were combined and used to develop a single maize (rice) model. However, the winter and spring wheat models were developed separately due to potential differences in the yield responses to environmental conditions because winter and spring wheat has different growing seasons and physiological characteristics. Model fitting was conducted using
the statistical software R (66) using the randomForest package (67) with the following settings: number of trees (ntree = 500), number of predictors sampled at each split (mtry = 3), and minimum size of terminal nodes (nodesize = 5). These hyperparameter settings were selected based on the package's default values.

For model validation purposes, we used the leave-one-year-out cross-validation technique and assessed the Pearson's correlation coefficients and the root mean squared errors calculated between the actual and simulated yield anomalies. Furthermore, the model evaluation was done to see whether the RF models with d4PDF factual climate simulation could reproduce the actual global pattern of the changes in yield anomaly SD derived from the grid yield dataset. We also checked whether values of d4PDF factual and counterfactual climate simulations did not exceed the ranges covered by the training data (i.e. the JCS retrospective meteorological forcing) to ensure the reliability of the RF models. Last, we assessed the importance of the explanatory variables. We used two methods to rank model input variables according to their importance: the permutation method and the node purity method that rely on the residual sum of squares and the Gini index as importance metrics, respectively.

## Climate variables

Temperature and precipitation were considered to be the primary climatic drivers of yield anomalies. To account the effects of short-term extreme temperatures, relative counts of hourly temperature exposure for a growing season were calculated using nine bins with constant 5 °C intervals (T < 0, 0 ≤ T < 5, 5 ≤ T < 10, ..., 25 ≤ T < 30, 30 ≤ T < 35, and 35 ≤ T). Hourly temperature values were obtained by fitting a sine curve to daily maximum and minimum temperatures. Relative counts of daily precipitation exposure during the growing season were also calculated using 5 bins with unequal intervals (P < 1, 1 ≤ P < 5, 5 ≤ P < 20, 20 ≤ P < 30, and 30 ≤ P). The growth periods were fixed based on the reported planting and harvesting dates (50). For winter wheat only, the period before the completion of vernalization was excluded from the calculation by identifying the first date at which the fraction of growing-degree days to the crop total thermal requirements exceeds 0.1; this was estimated using the crop phenology model with the crop-specific base (0 °C) and maximum temperature (26 °C) for winter wheat. These climate variables were calculated for each of the actual, factual, and counterfactual climate conditions.

## Soil variable

We focused on SOC content as a main soil variable, given the fact that in arid and semiarid regions of the world, high levels of SOC improve soil water holding capacity and moderate drought damage (37). The topsoil (0–30 cm) SOC data were obtained from the Reggridded Harmonized World Soil Database v1.2 (68) and aggregated to the 0.5° resolution. In addition to SOC, we also used fractions of clay, silt and sand, bulk density, cation exchange capacity of the clay fraction, pH, and available water storage capacity, all at the topsoil, to consider possible effects of soil characteristics on yield anomaly.

## Management variables

We considered two variables related to producer's crop management, i.e. irrigation and N fertilizer. Irrigation mitigates drought damages and heat stress through evaporative cooling of the crop canopy (69). The irrigation intensity used in this study considered changes due to the expansion of the irrigation-equipped

area based on the global historical irrigation dataset for the period of 1900–2005 (HID) (70) and crop-specific variations based on the global monthly irrigated and rainfed crop areas around the year 2000 (MIRCA2000) (71). The data were aggregated to the 0.5° resolution. Due to the lack of data, irrigation intensity values for the period 2006–2019 were extrapolated using a linear regression using the data from 1998 to 2005, although this assumption may not be totally accurate for some regions experiencing decreases in irrigation area, such as Russia (38). In addition to irrigation, we used the grid annual N application rate for each crop (72) as it is known that abiotic stresses become important as N inputs increase (35, 36).

## R&D variable

We incorporated the effect of agriculture R&D and resulting technologies and management practices that help stabilize yields into the RF models. To do so, an economic indicator was used. The indicator, Rᵢ,t, is the accumulated agricultural R&D expenditure for country i from 1995 to year t and expressed in constant 2015 USD. The accumulation took into account the duration of research and the obsolescence of technology:

$$R_{i,t} = E_{i,t-6} + (1 - \delta)R_{i,t-1} \quad (1)$$

where E is the country annual agricultural R&D expenditure (constant 2015 USD) and δ is the annual rate of technology obsolescence (= 0.1 or 10% per year) (72). The lag time between research and technology adoption by producers was set at 6 years, although lag time may vary by country and over time (73). As the R&D variable is of key importance in this study, we performed the sensitivity analysis on this variable, as described in the later subsection.

## Attribution analyses

Simulation experiments were conducted using the RF models and d4PDF factual and counterfactual climate simulations. As shown in Fig. 6, the fc and ct runs compared the factual and counterfactual climates and determine the climate effect to yield stability under actual R&D conditions. The fc and fc.em runs compared yield stability with and without R&D after 2000 under climate change.

## Climate effect

The changes in yield stability associated with human-induced climate change were detected through a comparison of fc and ct runs (Fig. 6). Using 100 samples of the slope value of yield anomaly SD derived from the RF models forced by d4PDF climate simulation, a histogram was derived for each of the factual and counterfactual yields (Fig. S13). When the average slope value of the factual yield anomaly SD is higher (lower) than that of the counterfactual yield anomaly SD, it indicates that climate change contributed to decreasing (increasing) yield stability in 2000–2019. The two-sided Wilcoxon rank-sum test (74) was used to test the null hypothesis that the two types of slopes were equal.

## R&D effect

The R&D effect after 2000 to yield stability under climate change was identified through a comparison of the fc and fc.em runs (Fig. 6). We assumed that the fc run, which included the time-dependent R&D and management variables, reflects advances in R&D and associated practices. Then, we performed a counterfactual run in which advances in R&D and management practices are lacking to varying degrees. To do this, we fixed the R&D and
management variables to constant values. In the first model specification (M1), values of the R&D and management variables were set at the 2000 levels to calculate yield anomalies from 2000 to 2019, while in the fifth model specification (M5), they were set equal to their yearly values until 2004 and were held constant thereafter (Fig. S14). This model setting resulted in 20 model specifications per crop (M1 to M20). For each year in which advances in R&D and practices were eliminated in the RF models, a histogram of the slope value for the corresponding model specification was derived and compared with a histogram of the slope value for actual R&D and practices, each consisting of 100 members from the RF models forced by d4PDF factual climate simulation (Fig. S15). When the average slope value for actual R&D and practice was significantly smaller than that without advances in R&D and practice, it indicates that changes in R&D and practices after 2000 increased yield stability. The absence of a significant difference indicates that changes in R&D and practices since 2000 did not have a significant effect on increasing yield stability under climate change.

## Sensitivity analysis

We examined the sensitivity of variable importance and attribution results to the different assumptions used in the R&D variable calculation. The country annual R&D expenditure data collected were for agriculture and not limited to the four crops considered here. Therefore, we first examined whether the use of R&D expenditure for these crops, not for agriculture, affected results. To that end, we assumed that annual R&D expenditure designated for the four crops was proportional to the share of value derived from production of the four crops to total value of agricultural production. We collected country annual data on value of production for the four crops and entire agriculture from FAO (48) and derived them for each country and year. Second, we considered shorter (6-year) and longer (12-year) time lags between research and adoption of technologies by producers. Finally, we considered relatively slower (10% per year) and faster (20% per year) technology obsolescence by setting $δ$ value in Eq. 1 to 0.1 and 0.2, respectively.

We performed five experiments in the sensitivity analysis: Ex. 1—R&D expenditure for agriculture (c0), shorter time lag between research and adoption (r6), and slower technology obsolescence (o1) (denoted “c0r6o1”); Ex. 2—R&D expenditure for agriculture, longer time lag (r12), and slower obsolescence (“c0r12o1”); Ex. 3—R&D expenditure for the four crops (c1), shorter time lag, and slower obsolescence (“c1r6o1”); Ex. 4—R&D expenditure for the four crops, longer time lag, and slower obsolescence (“c1r12o1”); and Ex. 5—R&D expenditure for the four crops, shorter time lag, and faster obsolescence (“c1r6o2”). For each experiment, the RF models were built and then the variable importance and the area shares for the climate and R&D effect categories were estimated.

In the main text, we presented the results from Ex. 1 (c0r6o1). We found that the assumptions used in the R&D variable calculation affected the level of but not the pattern of cumulative R&D expenditure (Fig. S16). The difference in the level of cumulative R&D expenditure was absorbed in the RF models, as ordinary regression models would do, and did not affect the variable importance and attribution results.

## Acknowledgments

T.I. thanks Jui-Min Wang for her analysis in the early stages of this study.

## Supplementary Material

Supplementary material is available at PNAS Nexus online.

## Funding

T.I. was supported by the Environment Research and Technology Development Fund (grant number: JPMEERF23S21120) of the Environmental Restoration and Conservation Agency provided by the Ministry of the Environment of Japan, and Grants-in-Aid for Scientific Research were received from the Japan Society for the Promotion of Science (grant numbers: 22H00577 and 23H00351). Y.I. and H.S. were supported by the Advanced Studies of Climate Change Projection (SENTAN, grant number: JPMXD0722680395) of the Ministry of Education, Culture, Sports, Science and Technology of Japan. D.M. and T.I. were supported by INRAE-NARO Linkage Call 2023.

## Author Contributions

T.I. conceptualized the study, developed the methodology, and conducted formal analysis. T.S., Y.M., K.O., and T.T. curated crop data and assisted with the investigation. H.S. and Y.I. curated climate data and assisted with the investigation. D.M. assisted with methodology development. All authors contributed to writing the original draft.

## Preprints

This manuscript was posted on a preprint: https://doi.org/10.21203/rs.3.rs-3895355/v1.

## Data Availability

The country yield and area harvested data from the FAO statistical database can be obtained from https://www.fao.org/faostat/en/#home. Gross domestic R&D spending data for some developed countries can be obtained from the OECD's MSTI database: https://data-explorer.oecd.org/. Agricultural R&D expenditure data for some developing countries are accessible from the IFPRI's ASTI database: https://doi.org/10.7910/DVN/9OXBIB. The grid yield dataset is accessible at the Data Integration and Analysis System (DIAS): https://doi.org/10.20783/DIAS.564. The actual climate data and the bias-corrected d4PDF factual and counterfactual climate simulations can be, respectively, obtained from https://doi.org/10.20783/DIAS.523 and https://doi.org/10.20783/DIAS.544.

## References

1. FAO, IFAD, UNICEF, WFP, WHO. *The state of food security and nutrition in the world* 2020. Rome: FAO, 2020.

2. Nóia Júnior RS, et al. 2022. Needed global wheat stock and crop management in response to the war in Ukraine. *Glob Food Secur.* 35:100662.

3. Aday S, Aday MS. 2020. Impact of COVID-19 on the food supply chain. *Food Qual Saf.* 4(4):167-180.

4. Gaupp F. 2020. Extreme events in a globalized food system. *One Earth*. 2(6):518-521.

5. Zhu P, et al. 2022. Warming reduces global agricultural production by decreasing cropping frequency and yields. *Nat Clim Change*. 12(11):1016-1023.
6 Wei D, Gephart JA, Iizumi T, Ramankutty N, Davis KF. 2023. Key role of planted and harvested area fluctuations in US crop production shocks. *Nat Sustain.* 6(10):1177–1185.

7 Bezner Kerr R, et al. Food, fibre, and other ecosystem products. 2022. In: Pörtner H-O, et al., editors. *Climate change 2022: impacts, adaptation and vulnerability*. Cambridge and New York: Cambridge University Press. p. 713–906.

8 Iizumi T, Ramankutty N. 2016. Changes in yield variability of major crops for 1981–2010 explained by climate change. *Environ Res Lett.* 11(3):034003.

9 Ben-Ari T, et al. 2018. Causes and implications of the unforeseen 2016 extreme yield loss in the breadbasket of France. *Nat Commun.* 9(1):1627.

10 Thirumalai K, DiNezio PN, Okumura Y, Deser C. 2017. Extreme temperatures in Southeast Asia caused by El Niño and worsened by global warming. *Nat Commun.* 8:15531.

11 Beillouin D, Schauberger B, Bastos A, Ciais P, Makowski D. 2020. Impact of extreme weather conditions on European crop production in 2018. *Philos Trans R Soc Lond B Biol Sci.* 375(1810):20190510.

12 White RH, et al. 2023. The unprecedented Pacific Northwest heat-wave of June 2021. *Nat Commun.* 14(1):727.

13 Cooper M, Messina CD. 2023. Breeding crops for drought-affected environments and improved climate resilience. *Plant Cell.* 35(1):162–186.

14 Berrang-Ford L, et al. 2021. A systematic global stocktake of evidence on human adaptation to climate change. *Nat Clim Change.* 11(11):989–1000.

15 Lobell DB, et al. 2014. Greater sensitivity to drought accompanies maize yield increase in the U.S. Midwest. *Science.* 344(6183):516–519.

16 Miller SJ, Clarke E, Logan Mathews S. 2023. Little evidence of avoided yield loss in US corn when short-term forecasts correctly predict extreme heat. *Environ Res Lett.* 18(12):124041.

17 Mizuta R, et al. 2017. Over 5,000 years of ensemble future climate simulations by 60-km global and 20-km regional atmospheric models. *Bull Am Meteorol Soc.* 98(7):1383–1398.

18 Otto FEL. 2017. Attribution of weather and climate events. *Annu Rev Environ Resour.* 42(1):627–646.

19 Gillett NP, et al. 2021. Constraining human contributions to observed warming since the pre-industrial period. *Nat Clim Change.* 11(3):207–212.

20 Romanovska P, Undorf S, Schauberger B, Duisenbekova A, Gornott C. 2024. Human-induced climate change has decreased wheat production in northern Kazakhstan. *Environ Res Clim.* 3:031005.

21 Funk C, et al. 2018. Anthropogenic enhancement of moderate-to-strong El Niño events likely contributed to drought and poor harvests in Southern Africa during 2016. *Bull Am Meteorol Soc.* 99(1):S91–S96.

22 Verschuur J, Li S, Wolski P, Otto FEL. 2021. Climate change as a driver of food insecurity in the 2007 Lesotho-South Africa drought. *Sci Rep.* 11(1):3852.

23 FAO. *The impact of disasters on agriculture and food security 2023—avoiding and reducing losses through investment in resilience*. Rome: FAO, 2023.

24 Ray DK, Mueller ND, West PC, Foley JA. 2013. Yield trends are insufficient to double global crop production by 2050. *PLoS One.* 8(6):e66428.

25 Yu Q, et al. 2020. A cultivated planet in 2010—part 2: the global gridded agricultural-production maps. *Earth Syst Sci Data.* 12(4):3545–3572.

26 Fuglie KO, Morgan S, Jelliffe J. *World agricultural production, resource use, and productivity, 1961–2020*. U.S. Department of Agriculture, Economic Research Service, 2024. Report No.: EIB-268.

27 DeBruin JL, Schussler JR, Mo H, Cooper M. 2017. Grain yield and nitrogen accumulation in maize hybrids released during 1934 to 2013 in the US Midwest. *Crop Sci.* 57(3):1431–1446.

28 Ahrends HE, et al. 2018. Genetic yield gains of winter wheat in Germany over more than 100 years (1895–2007) under contrasting fertilizer applications. *Environ Res Lett.* 13(10):104003.

29 Shew AM, Tack JB, Nalley LL, Chaminuka P. 2020. Yield reduction under climate warming varies among wheat cultivars in South Africa. *Nat Commun.* 11(1):4408.

30 Schmidt M, Felsche E. 2024. The effect of climate change on crop yield anomaly in Europe. *Clim Resil Sustain.* 3(1):e61.

31 Kim W, Iizumi T, Nishimori M. 2019. Global patterns of crop production losses associated with droughts from 1983 to 2009. *J Appl Meteorol Climatol.* 58(6):1233–1244.

32 Vogel E, et al. 2019. The effects of climate extremes on global agricultural yields. *Environ Res Lett.* 14(5):054010.

33 Tigchelaar M, Battisti DS, Naylor RL, Ray DK. 2018. Future warming increases probability of globally synchronized maize production shocks. *Proc Natl Acad Sci U S A.* 115(26):6644–6649.

34 Kim Y-U, et al. 2024. Mechanisms and modelling approaches for excessive rainfall stress on cereals: waterlogging, submergence, lodging, pests and diseases. *Agric For Meteorol.* 344:109819.

35 Müller C, et al. 2018. Global patterns of crop yield stability under additional nutrient and water inputs. *PLoS One.* 13(6):e0198748.

36 Ahrends HE, et al. 2021. Nutrient supply affects the yield stability of major European crops—a 50 year study. *Environ Res Lett.* 16(1):014003.

37 Qiao L, et al. 2022. Soil quality both increases crop production and improves resilience to climate change. *Nat Clim Change.* 12(6):574–580.

38 Mehta P, et al. 2024. Half of twenty-first century global irrigation expansion has been in water-stressed regions. *Nat Water.* 2(3):254–261.

39 Zaveri E, Lobell DB. 2019. The role of irrigation in changing wheat yields and heat sensitivity in India. *Nat Commun.* 10(1):4144.

40 IFPRI. *Agricultural science and technology indicators: 2019 global food policy report annex table 1*. IFPRI, 2019.

41 Nalley L, Dixon B, Chaminuka P, Naledzani Z, Coale MJ. 2018. The role of public wheat breeding in reducing food insecurity in South Africa. *PLoS One.* 13(12):e0209598.

42 Wu H, et al. 2023. AsiaRiceYield4km: seasonal rice yield in Asia from 1995 to 2015. *Earth Syst Sci Data.* 15(2):791–808.

43 Ruane AC, et al. 2021. Strong regional influence of climatic forcing datasets on global crop model ensembles. *Agric For Meteorol.* 300:108313.

44 Mengel M, Treu S, Lange S, Frieler K. 2021. ATTRICI v1.1—counterfactual climate for impact attribution. *Geosci Model Dev.* 14(8):5269–5284.

45 Sloat LL, et al. 2020. Climate adaptation by crop migration. *Nat Commun.* 11(1):1243.

46 Zhang G, et al. 2017. Spatiotemporal patterns of paddy rice croplands in China and India from 2000 to 2015. *Sci Total Environ.* 579:82–92.

47 Maggi F, Tang FHM, la Cecilia D, McBratney A. 2019. PEST-CHEMGRIDS, global gridded maps of the top 20 crop-specific pesticide application rates from 2015 to 2025. *Sci Data.* 6(1):170.

48 FAO. FAOSTAT. FAO, 2023. [accessed 2025 Mar 28]. https://www.fao.org/fastat/en/#home
49 Iizumi T, Sakai T. 2020. The global dataset of historical yields for major crops 1981–2016. *Sci Data*. 7(1):97.

50 Sacks WJ, Deryng D, Foley JA, Ramankutty N. 2010. Crop planting dates: an analysis of global patterns. *Glob Ecol Biogeogr*. 19(5): 607–620.

51 Iizumi T, et al. 2024. Increasing heavy rainfall events and associated excessive soil water threaten a protein-source legume in dry environments of West Africa. *Agric For Meteorol*. 344:109783.

52 Ishii M, Mori N. 2020. d4PDF: large-ensemble and high-resolution climate simulations for global warming risk assessment. *Prog Earth Planet Sci*. 7(1):58.

53 Mizuta R, et al. 2012. Climate simulations using MRI-AGCM3.2 with 20-km grid. *J Meteorol Soc Jpn*. 90A:233–258.

54 Shiogama H, et al. 2016. Attributing historical changes in probabilities of record-breaking daily temperature and precipitation extreme events. *Sci Online Lett Atmos*. 12:225–231.

55 Imada Y, Watanabe M, Kawase H, Shiogama H, Arai M. 2019. The July 2018 high temperature event in Japan could not have happened without human-induced global warming. *Sci Online Lett Atmos*. 15A:8–12.

56 Kawase H, et al. 2019. Contribution of historical global warming to local-scale heavy precipitation in western Japan estimated by large ensemble high-resolution simulations. *J Geophys Res Atmos*. 124(12):6093–6103.

57 Imada Y, et al. 2020. Advanced risk-based event attribution for heavy regional rainfall events. *NPJ Clim Atmos Sci*. 3(1):37.

58 Iizumi T, et al. 2018. Crop production losses associated with anthropogenic climate change for 1981–2010 compared with pre-industrial levels. *Int J Climatol*. 38(14):5405–5417.

59 Sultan B, Defrance D, Iizumi T. 2019. Evidence of crop production losses in West Africa due to historical global warming in two crop models. *Sci Rep*. 9(1):12834.

60 OECD. MSTI database. OECD, 2023. [accessed 2025 Mar 28]. https://stats.oecd.org/Index.aspx?DataSetCode=MSTI_PUB

61 World Bank. World development indicators. World Bank, 2023. [accessed 2025 Mar 28]. https://databank.worldbank.org/reports.aspx?source=world-development-indicators#

62 Hipel KW, McLeod AI. Time series modelling of water resources and environmental systems. Amsterdam: Elsevier, 1994.

63 Breiman L. 2001. Random forests. *Mach Learn*. 45(1):5–32.

64 Folberth C, et al. 2019. Spatio-temporal downscaling of gridded crop model yield estimates based on machine learning. *Agric For Meteorol*. 264:1–15.

65 Chen M, Guilpart N, Makowski D. 2024. Comparison of methods to aggregate climate data to predict crop yield: an application to soybean. *Environ Res Lett*. 19(5):054049.

66 R Core Team. R: a language and environment for statistical computing. Vienna: R Foundation for Statistical Computing, 2023.

67 Liaw A, Wiener M. 2002. Classification and regression by randomForest. *R News*. 2(3):18–22.

68 Wieder WR, Boehnert J, Bonan GB. 2014. Evaluating soil biogeochemistry parameterizations in Earth system models with observations. *Glob Biogeochem Cycles*. 28(3):211–222.

69 McDermid S, et al. 2023. Irrigation in the Earth system. *Nat Rev Earth Environ*. 4(7):435–453.

70 Siebert S, et al. 2015. A global data set of the extent of irrigated land from 1900 to 2005. *Hydrol Earth Syst Sci*. 19(3):1521–1545.

71 Portmann FT, Siebert S, Döll P. 2010. MIRCA2000—global monthly irrigated and rainfed crop areas around the year 2000: a new high-resolution data set for agricultural and hydrological modeling. *Glob Biogeochem Cycles*. 24(1):GB1011.

72 Iizumi T, et al. 2017. Responses of crop yield growth to global temperature and socioeconomic changes. *Sci Rep*. 7(1):7800.

73 Baldos ULC, Fuglie KO, Hertel TW. 2020. The research cost of adapting agriculture to climate change: a global analysis to 2050. *Agric Econ*. 51(2):207–220.

74 Hollander M, Wolfe DA, Chicken E. Nonparametric statistical methods. 3rd ed. Hoboken: John Wiley & Sons, 2015.