Significant Reductions in Crop Yields From Air Pollution and Heat Stress in the United States

The joint exposure of plants to surface ozone, atmospheric aerosols, and heat stress can lead to considerable decreases in crop yields. Surface ozone negatively impacts plant photosynthesis while aerosols can have positive or negative effects from its dual impact on light and temperature. Here, using a statistical model, we show that in the United States, as a result of improvements in air quality, the damages caused by ozone and aerosols have decreased since 1980. Historically, relative yield losses due to ozone were 8.7% and 4.8%, and due to aerosols were 11.3% and 23.2% for maize and soybean, respectively. Maize yields are more sensitive to ozone pollution while soybean yields are more sensitive to aerosol pollution. In future RCP 8.5 scenario, absent significant reductions in emissions or improvements in air quality, maize, and soybean would have on average, 58.5% and 36.9% additional yield reductions, respectively, mainly caused by warming. Future climate warming and fossil fuel combustion driven changes to air pollution may have differing impacts on crop yield and should be jointly considered in any assessment of U.S. food security.

Interestingly, they also pointed out that aerosols could reduce soil evapotranspiration and ease plant water stress. This reduction in water loss from evapotranspiration may moderate the lessening of photosynthesis and cause aerosols to have a beneficial effect on final yields in rainfed farmland.
Air pollution levels also coincide with global climate warming, as both are outcomes of fossil fuel related emissions. Lobell et al. (2011) indicated that maize and wheat production declined by 3.8% and 5.5% across the world due to recent climate change, respectively. As for certain nations, those tendencies were large enough to repress sharp growth rate in average yields that arose from technological application, carbon dioxide fertilization (CDF), and other factors. However, anthropogenic aerosols through local cooling effects can overcome heat stress effects and lead to increased plant productivity (Zhang, Goll, et al., 2019). Such trade-offs are meaningful to quantify and pave the essential path for policy approaches to mitigate loss of yields caused by heat stress. Moreover, many studies also have presented the impacts of meteorological parameters on ozone (Bloomer et al., 2009;Han et al., 2020). Increasing weather extremes might fuel ozone pollution but could be influenced by aerosols arising from increased fossil fuel combustion.
While prior studies have examined specific pollutants or climate trends over specific crops and regions, no study has looked at the joint and combined impact of all three (ozone, aerosols, and heat stress) including interactions over a large region for multiple crops over time periods of historic changes in all three and in management related improvements in yield. Here, we obtained historical annual variability of two major air pollutants (i.e., ozone and aerosol) and heat stress, and quantified Pearson correlation among these values and crop yields. Based on these results, we then developed and evaluated an empirical regression model to identify their impacts on crop yield, including maize and soybean, in the United States from 1980 to 2019. We hypothesized that crop yields have suffered from heat stress; however, given air pollution-climate interactions, part of that yield reduction might be compromised or exacerbated in regions suffering poor air quality. We also analyzed how much change of crop yields could be influenced by future temperature and precipitation and test the sensitivity of crop yield in response to air pollution changes on the basis of the statistical model, assuming historical trajectories in air quality. Based on well-established literature (Hong et al., 2020;Lobell & Asner, 2003;McGrath et al., 2015) and pollutant-temperature interaction assumptions (see Materials and Methods for details), we chose the regression model found in Equation 4.
There are several approaches to estimate crop yields responses to climate and air pollution at regional scales, for example, dose-response functions (Mills et al., 2018). However, extrapolating over larger regions on the basis of limited observations at a small number of field locations may be unreliable and bring numerous uncertainties. Statistical regression models that include climatic variables as well as air pollutant metrics can account for potential interactions, possibly giving a better assessment of compound damages or tradeoffs (Hong et al., 2020;Lobell et al., 2011).

Data Collection
Annual yields in bushels per acre (bu 1 acre ) of maize and soybean for counties within the United States were provided by the U.S. Department of Agriculture NASS for the years 1980-2019 (NASS USDA, USDA/ NASS QuickStats query tool, 2014) and converted to kilogram per hectare 1 (kg ha ). The relationships of agricultural production, climate, aerosol, and ozone were analyzed using historical yields, weather, aerosol optical depth (AOD), and ground ozone data, respectively. The monthly average of the daily maximum temperature (T MAX ) along with precipitation (PRCP) at 4,377 stations in the U.S. were obtained from the Global Historical Climatology Network (GHCN) Global Summary of the Month (GSOM), Version 1 (Lawrimore, 2016). We omit the difference between rainfed and irrigating crop yield to focus on total yields. Y. Li et al. (2020) showed that irrigation leads to a considerable cooling on daytime land surface temperature, an increase in the enhanced vegetation index, and 81% higher maize yields compared to rainfed maize. To verify that our results are still valid without accounting for this effect, we made an additional analysis using counties that have separated and available irrigated and rain-fed yields data to test the effects of irrigation, see Supporting Information Text S2 and Figure S8.
Gridded three-hourly instantaneous aerosol optical depth (AOD) at 550 nm was obtained from MERRA-2 reanalysis (Global Modeling and Assimilation Office (GMAO), 2015) and then aggregated to local daylight (8:00-19:59) daily averages from 1980 to 2019. The gridded data originally had a spatial resolution of    0.5 0.625 and was averaged to the county level. The mean of the AOD data for all cells within a county was used as an estimate of the mean AOD for that county. AOD indicates the extinction and absorbing capability of atmospheric aerosols, which serves as a proxy of suspended aerosol concentrations. It provides a more consistent way to account for aerosol-induced impacts on incoming solar radiation than using ground fine particulate matter 2.5 (PM ) observations. MERRA-2 is a reanalysis using Goddard Earth Observing System, version 5 (GEOS-5), an earth system model, and assimilates MODIS and MISR AOD observations (Gelaro et al., 2017). Gueymard and Yang (2019) indicated that globally MERRA-2 performed better than CAMS AOD products, suggesting it is appropriate to use the MERRA-2 AOD products in a long study period compared to others. Furthermore, Song et al. (2018) compared MERRA-2 reanalysis products with surface observations of the diurnal and seasonal variability of AOD and indicated MERRA-2 captures changes of daytime AOD variations.
Hourly ozone data at 3,905 stations in the study region from the Environmental Protection Agency's Air Quality System (EPA AQS) were employed to calculate three widely used ozone cumulative indices projected to measure the plants exposure to ozone (that is, W126, AOT40, and SUM06) (US EPA, 2009;Air Quality System [AQS]). In each metric, lower ozone values are given less weight. These metrics explain that ozone exposure in longer time and lower amount poses smaller threat to vegetation than ozone exposure in shorter time and higher amount. We use the three indices (hereafter, OI) to elucidate ozone damaging on crop yield. W126 is a cumulative indicator of hourly ozone concentrations weighted by a sigmoidal scale. AOT40 and SUM06 are cumulative indicators of hourly ozone concentrations exceeding 0.04 and 0.06 ppm thresholds, respectively. The calculations are in the follows, where h O is the hourly ozone concentration in ppm per hour h during the daytime (8:00-19:59), and n is the total number of hours within a period. To adjust for missing hourly observations, the cumulative indices are summed up within a certain time period, this result is then divided by the numbers of hours in its reported period, and in the end, the result will be multiplied with total hours in the daytime (Hong et al., 2020). Some studies suggest that ozone flux is a better indicator of damaging plants than ozone concentration metrics, as the effects of ozone on the plant mainly depend on its stomatal uptake (Kangasjärvi et al., 1994). Pleijel et al. (2004) and Peng et al. (2019) illustrated that the ozone exposure-and flux-yield relationships with maize and wheat were similar. Thus, we believe that using ozone concentration did not bring many errors into our estimation of yield impacts. We mainly reported results of AOT40 metric in Main text, and others could be seen in Supporting Information Figures S5 and S6. This study only refers to those counties with equal and more than 20 years of available data. Included counties are shown in Supporting Information Figure S3.
In order to assess the impacts of future climate change on crop yields, we combine historical sensitivities with future scenarios of climate from the Multivariate Adaptive Constructed Analogs (MACA) data set (Abatzoglou & Brown, 2012), which used a statistical downscaling method to bias-corrected climate models outputs from Coupled Model Inter-Comparison Project 5 (CMIP5). We use the future Representative Concentration Pathways (RCPs) RCP 4.5 and RCP8.5 scenarios (2051-2089) from these high resolution MACA data sets. To incorporate model uncertainty, we chose an ensemble of 18 different global climate models (GCMs) combinations (Table S4). All gridded model simulations are aggregated to the county level by calculating the mean of all the pixels covering one county.

Statistical Model
To evaluate contribution of changes in management, air quality, and climate on county-level maize and soybean yield, we first organized data in a panel format. We then developed a statistical yield model for both crops using the following regression model (Equation 4), where OI is growing season (June-August) average three ozone cumulative indices (in units of ppm h); AOD is growing season average aerosol optical depth (dimensionless); , T MAX syi and , PRCP syi here are seasonal mean of daily maximum temperature (C), and seasonal mean precipitation (mm), respectively; , T MAX yi and PRCP yi are growing-season maximum temperature (C) and mean precipitation (mm), respectively; the subscripts , y i and s are indices of year, county, and season, respectively;  is the error term. It is worth noting that we use different ozone and aerosol interaction with temperature and precipitation terms. This choice is because we would like to emphasize the aerosol may have more sophisticated and nonlinear effects on crop yields according our correlation analyses (see Results and Discussion for details). We use seasonal average weather variables instead of monthly variables in the model to partially avoid over-fitting. The four seasons are defined as winter (December-February), spring (March-May), summer (June-August), and autumn (September-November). It have been widely recognized the nonlinear effects of climate on agricultural yields (Hong et al., 2020;Lobell & Asner, 2003;McGrath et al., 2015;Ortiz-Bobea et al., 2019), thus consistent with prior studies, seasonal weather variables are presented in the quadratic polynomial form. The county-specific and overall time trend ( CY ) accounts for changes in agronomic practices over time that influence yields, assuming these changes are linear. County fixed effects ( C ) represented by dummy variables explain time-invariant factors that vary across counties, such as differences in soil quality and agronomic practices in different regions. We assume all effects are fixed. Similar models have been used in previous studies of air pollution and climate (Hong et al., 2020;Lobell & Asner, 2003;McGrath et al., 2015).
While growing season climate is the dominant driver of yields, shoulder season climate influences planting and harvest dates, while winter season climate can influence moisture availability in snow or soil frost depths. This is why we use seasonal average weather variables. Our model includes linear time trends that control for cross-sectional differences in management and soil quality across the county, though we recognize that such time control might be not able to represent all situations, and food yields also rely on the consumer demand and governmental policies that are varying and difficult to include in our empirical model (Hong et al., 2020;Tilman et al., 2011). We assume independent random variation among the counties and years. Counties close in space might share similar agricultural practices, climatic and air pollution conditions, leading to a spatial correlation that may lead to overconfidence in our model if not properly accounted (Chen et al., 2016;McGrath et al., 2015;Moore & Lobell, 2015). To test this effect, we regressed an additional spatial autocorrelation model (spatial lag model) to account for the spatial correlation (see Text S3), and do not find a significant effect of spatial correlation on our overall findings ( Figure S9) (Chen et al., 2016;Gelfand et al., 2010).

Data Analysis
We obtained historical annual variability of two major air pollutants (i.e., ozone and aerosol) and heat stress, and calculated Pearson correlations between detrended climate/air pollution and crop yields. We then employ the least absolute shrinkage and selection operator (LASSO) regression method (Tibshirani, 1996) to fit the model, in which predictors are chosen in an automatic manner using the "lars" package in R (Hastie & Efron, 2013), that is, thus independent from researcher choice. This method is well suited to high-dimensional studies and can be tuned to deal with collinearity. Such method works by shrinking all coefficients toward zero and adding a penalty on their size, such that only the stronger signals retain non-zero coefficients and is thought to improve prediction accuracy. We perform an out-of-sample cross-validation as a robustness check, using two-thirds of the data for training, and the other one-third for testing. We compute correlation coefficients (r) between actual and predicted yields for both the trained and tested data. Models in which the tested r was significantly lower than the trained r indicate the risk of over-fitting. Confidence intervals of r are derived from 1000 repeated tests with a random one-third of observations omitted. A similar method has been used in Hong et al. (2020). The results of over-fitting test and model performance were shown in Supporting Information Figures S1 and S2. The collinearity (VIF) and ozone-aerosol interaction analysis for variables is provided in Supporting Information Text S2 and Table S3.
Based on climate and air pollution sensitivities given by the statistical yield models, we then determine the crop yield responses as percentage changes in yields (i.e., relative yield change). For each crop, we calculated the whole country mean yield as the harvested area-weighted mean of yields from all counties and years. The averaged annual yield effects for all counties ozone and aerosols are estimated from the percentage differences between predictions using historical ozone and aerosol levels and hypothetical scenarios with zero ozone and aerosol pollution (J. Burney & Ramanathan, 2014;Hong et al., 2020;McGrath et al., 2015). The impact of historical temperature trends on crop yields is estimated from the difference between predictions using historical temperature and a hypothetical scenario using the 1980-2019 median temperature, to evaluate the impact of annual temperature variability during 1980-2019 (Hong et al., 2020). The averaged annual yield effects for all counties air pollution-climate interactions are estimated from the percentage differences between predictions using historical interactive data and a hypothetical scenario with zero interaction. We test the overall yield response to warming as the percentage difference between predictions using historical temperature (1980-2019) and the future climate scenario (2050-2089).
We acknowledge that separation of interactive effects from a regression model is challenging and to date it appears that no adequate estimation method has been established. Here, we do not intend to explicitly estimate the interactive effects of air pollution and climate, but need to consider them to select and evaluate an appropriate model. Our main focus is on direction of such effects, for example, whether aerosol pollution buffers or exacerbates high temperature impacts on crop yields. Therefore, we simply calculate the interactive effect in Figure 3 by the relative yield change against a conceptual model scenario where a certain interactive terms are assigned to be zero and a normal condition. High values for some of our interactive effect could be explained by the fact that when one interactive term is set to zero, the dominant driver is also zero. For example, in the temperature-ozone interaction term, the temperature is zero, which is not and leads to more than 100% RYC. However, the focus should be on the net change with this statistical model. Importantly it should be noted that the numeric value of interactive effect might vary with the choice of units in our variables (e.g., the temperature could use K or C), but the sensitivity, importance, and the direction of the interaction is insensitive to the choice of baseline or unit, since the model calibration of coefficients would account for that.
Finally, to determine the confidence intervals (CIs) related to these estimates, we use a bootstrapping method to produce distributions by using both crop models generated from re-sampling (with replacement) the observations for one county in one year 500 times, and the 5th and 95th or 10th and 90th percentiles among the 500 bootstrap replicates are selected as 90% and 80% CIs, respectively. The bootstrapping accounts for uncertainties in model formulation by choosing new predictors for each replicate. Non-standardized regression coefficients in median estimates from 500 times bootstrapping were shown in Supporting Information Table S1. The regions in U.S. are in compliance with the definitions of the Agricultural Research Service (ARS), the research arm of the USDA. The ARS has sectioned their work into five geographic regions: Midwest Area, Northeast Area, Pacific West Area, Plains Area, and Southeast Area ( Figure S11).

Historical Trends and Correlations
The United States underwent a significant improvement in crop yields over the past decades (Kucharik et al., 2020). Historical median yield increased from 5,804 to 10,046 kg 1 ha for maize, from 2,053 to 3,212 kg 1 ha for soybean from 1980 to 2019 (Figures 1a and 1c) ( 2 r = 0.218, 0.233 and  0.01 p ). Such patterns were substantially contributed from human interventions, for example, improvement in use of water, fertilizer, application of pesticides, and new varieties of species (Tilman et al., 2011).
Over the time period, climate changed (Figures 1a and 1c). While maximum daily temperature did not show significant trends 2 r = 0.001-0.002,  0.01 p ), precipitation demonstrated increasing trends (  0.01 p ), but with small effect sizes (maize: 0.30 mm 1 yr ; soybean: 0.27 mm 1 yr ). Increased precipitation in these agricultural hot spots arise both from intensification of the water cycle expected from climate change in these regions and possibly from increased transpiration occurring due to expansion of irrigation and intensification of cropping (Mueller et al., 2016;Portmann et al., 2009). Similar effects were found for air quality as for precipitation (Figures 1a and 1c) 1988, 1999, 2003, and 2012, ozone level was higher than surrounding years and associated with hot temperature and low precipitation, suggesting the interactions of ozone formation with meteorological conditions. Atmospheric aerosols have declined since 1980 ( 2 r = 0.248, 0.421, and  0.01 p ), but aerosol optical depth (AOD, see Methods) was anonymously high in 1983 and 1992, largely as a result of the eruptions of volcanoes (Proctor et al., 2018). Overall, regions planted by maize and soybean shared similar trends of climate and air pollution variability, but maize compared to soybean received a larger range of temperature, precipitation, ozone, and AOD (i.e., the medians are similar but CIs are larger for maize).
We found ozone negatively affected maize and soybean yields in the U.S. over the study period ( 2 r = 0.69, 0.47, and  0.01 p ) (Figures 1b and 1d). It is worth noting that the coincidence of negative impacts of ozone and temperature on crop yields can jointly reduce crop yields. In contrast to ozone, maize, and soybean yields did not exhibit linear responses to AOD (  0.1 p ) suggesting atmospheric aerosols have more sophisticated modifications to the environment (Yamasoe et al., 2006). For instance, aerosols can redistribute the ratio of direct beam and diffuse photosynthetically active radiation (PAR), which have different responses for photosynthesis in some plants (Cannich et al., 2012). C4 plants, such as maize, are less sensitive to the enhanced diffuse radiation compared to C3 plants, such as soybean, a theory consistent with our sensitivity tests (Text S1 and Figure S4). In general, we suggest regardless of technological improvement and genetic achievement to yields, natural and other human-made modifications (i.e., air pollution) played crucial roles in regulating crop production.
Our results on ozone damages are consistent with McGrath et al. (2015), which reported the estimated average yield losses from ozone were 9.8% for maize and 5.5% for soybean from 1980 to 2011. The lower estimate in our finding is likely because our study period is longer, and in recent years yield continued to increase, while air pollution has been reduced. The finding of ozone correlation to crop yield (Figures 1b and 1d) is also aligned with Betzelberger et al. (2012), which suggested soybean being exposed in ozone environment for three months caused a linear increase in antioxidant capacity while reducing leaf area, light absorption, specific leaf mass, primary metabolites, seed production, and harvest index.
Crop yield sensitivity to ozone is calculated by a comparison between a scenario that has one addition of ozone exposure (+1  (Heagle, 1989). For soybean, we found a sensitivity to AOT40 of 0.43%  (Text S1 and Table S2). Thus, our current estimated sensitivities of soybean yield to ozone aligned well with previously experimental and statistical studies, across multiple metrics. However, we find that soybean is almost half as sensitive to ozone exposure than maize in three dose metrics-based estimates, which could largely explain why maize is damaged by ozone more than soybean. This result also emphasizes the important role aerosols plays in interaction with ozone and climate impacts on crop yield. This contrasts with earlier studies such as McGrath et al. (2015), who did not account for the effects of aerosol to crop, and as a result found more similar sensitivities to ozone for maize and soybean. The higher ozone sensitivity of maize yield estimates in this study than field and laboratory studies indicated that using a statistical method over a large region may cause unrealistic results, but this disagreement also suggests that crop variants and regional conditions have vital impacts on response to air pollution exposure by plants.
Historical crop yield damage from aerosols has not received much attention in threatening food security in the United States. Previous studies focused on aerosol damages were mostly specific to Asia. For example, J. Burney and Ramanathan (2014)   and (c-d) represent maize and soybean, respectively. The damages refer to relative yield change (RYC). For air pollution (a and c), ozone index (OI) indicates ozone index (only shows AOT40), aerosol optical depth (AOD) indicates aerosol optical depth. Error bars span the 90% confidence interval from 500 times bootstrapping. Red and blue lines are linear fitting for OI and AOD, respectively. Determination coefficients 2 r and p values of linear fittings are indicated. For temperature (b and d), we assume the baseline is median historical temperature . The 0% and −5% change in yield lines are highlighted in blue. than they otherwise would have been, absent climate and pollutant emissions trends, with some densely populated states experiencing 50% relative yield losses. Our model also finds that soybean yield is more sensitive to atmospheric aerosols than maize. Although studies have noted that C4 plant productivity is LIU AND DESAI 10.1029/2021EF002000 9 of 14 . Countrywide air pollution-related variables' relative yield change (RYC, %) to maize and soybean yield from 1980 to 2019. Note: noAOD represents the total aerosol RYC, noAODT means aerosol-temperature RYC, noAODP means aerosol-precipitation RYC; noOI represents the total ozone RYC, noOIT means ozone-temperature RYC, noOIP means ozone-precipitation RYC. The horizontal solid black line and the darker/lighter colored bands correspond to the median and 90/80% confidence interval from 500 times bootstrapping, respectively. Zero RYC line is dashed. AOT40 was used here. The calculation is derived from relative change from a conceptual scenario that certain variables are assigned to be zero and normal condition. inherently more sensitive to the reduction in direct light than C3 (Kanniah et al., 2012), maize has higher canopy height than soybean, which is thought to increase light use efficiency under diffuse light found in higher aerosol loading (Emmel et al., 2020).

Heat Stress Damages
Heat stress also influenced crop yield over the study period. Our results in Figures 1b and 1d demonstrate negative correlations between T MAX and crop yield, both for maize and soybean ( 2 r = 0.46, 0.16 and  0.01 p ), representing the significant damage of heat stress. We also show positive impacts of precipitation on maize and soybean yields ( 2 r = 0.64, 0.32 and  0.05 p ). The Southeast suffered the most heat stress effect (Figures 4c and 4f). In most years, maize and soybean had been significantly reduced by heat stress from 1980 to 2019 in the United States (Figures 2b and 2d). The variations in yields due to temperature are mostly within 10% and do not exhibit clear trends. Maize yield has been less reduced by heat stress with positive effects of temperature in 15 years, while soybean has suffered more frequent heat stress with 28 years of damages. These extreme temperature years are also associated with severe ozone pollution (Figures 1a and 1c), which could play an important role in interacting with heat stress effects on yield.

Interactive and Combined Effects
Consistent with this idea, we found ozone and temperature interactions effects were significant on both crop yields. Each interaction contribution to the crop yield was calculated by the comparison of normal predication and conceptual zero interaction, and divided by normal predication (e.g., OIT = yields yields yields normal zero OI TMAX and OI TMAX norma It is noted that the following numeric value should not be considered as exact estimates of interactive effects, but with a main focus on the direction or sensitivity of each effect (see Materials and Methods for details). Maize and soybean yields are reduced by 297.92% and 204.65%, respectively as a result from the ozone-temperature interaction term (Figure 3). Our findings also show significant effects of ozone-precipitation interaction on soybean yields (11.51%, 80% CI: 8.77%-13.88%), but insignificant on maize yields (−1.66%, 80% CI: −3.39% to 0.47%). Atmospheric aerosols and temperature interactions have caused increases of 162.83% (80% CI: 90.42%-243.66%) and −3.16% (80% CI: −80.14% to 45.57%) for maize and soybean yields, respectively. Maize shows an increase of 1.85% (−4.45% to 9.27%), but soybean has been damaged 27.48% (80% CI: 13.07%-47.95%) by aerosol-precipitation interactions. We only show the results of using unit of C (Figure 3), other unit (e.g., K, see in Figure S10) or baseline (not zeroing of interactive terms, but giving other constant values) could result in different numeric values in calculation. However, the sensitivity, importance, and the direction of the interaction is insensitive to the choice of baseline or unit, since the model calibration of coefficients would account for that. Caution must be given in interpreting the absolute parameter values of the interactive effects.
When we put air quality and heat stress together, we find that heat stress induced reduction has been compromised (by aerosols) or exacerbated (by ozone) for soybean yields (Figure 3). However, soybean and maize experienced the unequal magnitude of these interactive effects, although the aerosol cooling effect in maize yields is not as significant as soybean. A net temperature-associated interaction led to −207.81% yield change for soybean. Meanwhile, a net −135.08% maize yield change has been led by temperature-linked interaction. The direction of these interactions indicate the air pollution has jointly damaged crop yields with temperature.

Implications for Future Crop Productivity
To examine the impacts of future global warming to crop yields, we used climate projections from 18 global climate models in future climate scenarios (2050-2089) by combining our estimated historical responses to compare with the baseline (1988-2019) in CMIP5 RCP 4.5 and RCP 8.5 ( Figure 5). Compared with current climate, in future scenarios, both maize and soybean yields would be damaged by 52.60% and 30.58%, respectively, on average, under RCP 4.5. As expected, larger damages (58.45% and 36.87%) caused by future climate change under RCP8.5. These reductions are mainly contributed by future warming (more than 95%), and slightly by changes of precipitation. Warming also results in considerable increases of temperature-associated interactions (aerosol: 2%; ozone: 20%-50%) under RCP8.5 ( Figure 5). This result indicates that absent with air pollution changes, heat stress might compound with ozone pollution to reduce crop yields. Owing to the larger uncertainty of projections of precipitation, the precipitation-associated interactions are relatively unclear.
Future warming would also go along with rising atmospheric  2 CO 2 CO fertilization may partially offset the warming damages. A recent study indicated that in temperate regions, 2 CO emissions would cause a substantial short-term increase in agricultural production due to the instantaneous response of atmospheric 2 CO concentrations to emissions, but the net effect attenuates within a decade or so to near-zero as the 2 CO fertilization is canceled by the impacts of 2 CO -induced climate change (Shindell et al., 2019). Greater atmospheric 2 CO concentration induced by future fossil fuel combustion brings an offset of global warming on crop yields, and should be integrated into future analysis. Another vital factor subjected to be taken into account is that climate warming could significantly advance spring phenological development and facilitate earlier sowing, especially in mid and high latitudes (Buermann et al., 2018). Human managements (e.g., irrigation) could cause uncertain influences on future food security regionally and globally. Irrigation would supply extra water and have a cooling effect (Y. Li et al., 2020), especially in summer. With increased number of counties reporting these data, it would be excellent future work to develop a framework for further separating yield responses by irrigation type.
LIU AND DESAI 10.1029/2021EF002000 11 of 14 Figure 5. Countrywide soybean (a) and maize (b) yield changes (%) between future (2050-2089) and current (1980-2019) climate, using projections from 18 global climate models in RCP45 and 85. Each dot represents a particular GCM in CMIP5, downscaled by MACA. Vertical lines around each dot represents the 95% confidence interval from 500 times bootstrapping. The horizontal solid red line and the colored bands correspond to the median and limits of each ensemble, respectively. T represents temperature contribution changes including interactions. P represents precipitation contribution changes including interactions. OIT means ozone-temperature contribution changes, OIP means ozone-precipitation contribution changes; AODT means aerosol-temperature contribution changes, AODP means aerosol-precipitation contribution changes.

Conclusion
Overall, our analysis demonstrates a clear impact of air pollution and heat stress on reducing maize and soybean yields from 1980 to 2019 in the United States, but with crop specific responses. Maize and soybean are more sensitive to ozone and aerosol than each other, respectively, though it is linked with local air pollution loading as well. We highlight the importance and needs of including the interactions between air pollution and climate. Our results implicate larger ozone-induced crop yield losses in the future than earlier estimated due to the additional account of ozone-temperature interaction, but also a larger role of ozone regulation in mitigating the warming-induced effects on agriculture and securing the future food supply. From Figures 2a and 2c, we see the consequence of clean air policy, which shows reductions but less so in soybean. Emissions of black carbon (BC), sulfate, etc., aerosols may require further reductions in the regions with an abundance of soybean to decrease impact on yields. Crop yield losses from air pollution have been decreasing since 1980 and these trends remain nowadays with the damages ranging from 5% to 20%. It can be noted that controls over air pollution has generated sufficient benefits to agricultural production, and further reductions in air pollutant emissions would lead to further benefits in crop yields (Figures 4a,4b,4d and 4e). Through the practice of air quality regulation, the U.S. has set a positive example for other parts of the world that may still be suffering from substantial losses in crop yield from air pollution (K. Li et al., 2019;Zhang, Zheng, et al., 2019), conclusively showing an additional co-benefit of air quality improvement on improving food security.

Data Availability Statement
All historical data used are publicly available and open access, with the data sources list in the Materials and Methods. The MACA data could be accessed through https://climate.northwestknowledge.net/MACA/. Replication code is available at https://doi.org/10.5281/zenodo.4922064.