Quantifying the indirect effects of different air pollutants on crop yields in North China Plain

High concentrations of air pollutants could affect crop yield directly via influencing crop physiology and indirectly via regulating climate conditions. How multiple air pollutants affect crop yields directly and indirectly remains elusive in the North China Plain (NCP), which is the hotspot of air pollution and crop production. Here, we construct data-driven models to quantify the indirect effects of four major air pollutants on the yields of winter wheat and summer maize through modulating climate variations in the NCP during 2005–2020. Our results show that wheat has a greater negative sensitivity to absorbing aerosol optical thickness (AAOD), ozone concentration (O3), sulfur dioxide concentration (SO2) and nitrogen dioxide concentration (NO2) than maize. The indirect effects of AAOD, O3 and SO2 in November and O3 in April control wheat yield variation, while the indirect effects of AAOD and O3 in June, O3 and SO2 in July, O3 and NO2 in August, and O3 in September dominate maize yield variation. Our results emphasize the indirect effects of air pollutants on crop yield via regulating climate variations, which has great implications for improving our understanding of air pollution-climate-crop interactions and guiding targeted adaptation and mitigation efforts.


Introduction
High concentrations of air pollutants are not only detrimental to human health but also have large impacts on crop production.Based on scale site experiments, the damages of elevating surface ozone on crops were first documented by research projects conducted in the United States (US) in the 1980s (Heck 1989) and Europe in the 1990s (Krupa et al 1998).Subsequently, more and more studies have quantified the negative impacts of air pollution on crop yields worldwide (Avnery et al 2011a, 2011b, Feng et al 2015, Mcgrath et al 2015, Pleijel et al 2019, Stevens et al 2019, Zhao et al 2020).Recent studies suggest that air pollutants may have contributed more to crop yield variations than climate change (Burney and Ramanathan 2014, Tai et al 2014, He et al 2022).For example, Burney and Ramanathan (2014) illustrated that 90% of the relative yield change loss in wheat between 1980-2010 can be attributed to short-lived climate pollutants in India.Tai et al (2014) showed that the effect of ozone is much larger than the temperature effect under RCP 4.5 and RCP 8.5 scenarios for the major wheat-producing regions around the world.He et al (2022) demonstrated that the effects of particulate pollution mitigation surpassed climate change impacts and controlled the crop yield variation during 2006-2018 in most provinces of China.
The negative effects of air pollutants can be exerted by entering crop leaves through stomata or adhering to crop leaves, which could directly accelerate crop leaves' senescence and reduce crop photosynthesis and yield (Agrawal et al 2003, Bell et al 2011, Singh et al 2014, Schauberger et al 2019).Besides the direct impacts, changes in air pollutants could indirectly influence crop yield via regulating climate variations (e.g.radiation, temperature, precipitation) during critical crop growth stages.For example, ozone, which is an important greenhouse gas in the atmosphere, not only affects temperature but also modulates the formation of clouds with subsequent influences on radiation and precipitation (Hauglustaine et al 1994, Shindell et al 2005, Fry et al 2012, Bai et al 2016).Aerosols which are categorized into absorbing and scattering groups have important feedback on regional climate (Li et al 2022b).Specifically, absorbing aerosols (e.g.black carbon, dust, organic carbon, etc.) can heat the atmosphere and alter atmospheric stability by absorbing incoming solar radiation, while scattering aerosols composed of S-oxide and N-oxide can reflect radiation to cool the atmosphere as well as the surface (Wang 2013, Sanap and Pandithurai 2015, Kuniyal and Guleria 2019).Besides, both absorbing and scattering aerosols can affect the optical properties and lifetime of clouds by serving as cloud condensation nuclei (Lohmann andFeichter 2005, Huang et al 2006).In addition, higher diffuse radiation induced by scattering aerosols can improve the light use efficiency of crops, but the resultant lower direct radiation could cause a yield loss (Zhang et al 2017(Zhang et al , 2018b)).Using satellite observations, He et al (2023) showed that aerosols have large impacts on plant photosynthesis in Europe by perturbing radiation, temperature, vapor pressure deficit and soil moisture.Therefore, understanding the direct and indirect effects of multiple air pollutants on crop yield is critical for targeted air pollution control measures, which have large benefits for global and regional crop production (Tai et al 2014, Zhang et al 2018b, Lobell and Burney 2021).Van Dingenen et al (2009) estimated the global impact of surface ozone and found that 40% of global crop yield loss occurs in China and India.Mills et al (2018a) estimated that ozone during 2010-2012 caused a global yield reduction annually by 12.4%, 7.1%, 4.4% and 6.1% for soybean, wheat, rice and maize, respectively.Lobell et al (2022) evaluated the global effects of nitrogen dioxide (NO 2 ) on crop yield and estimated that reducing NOx levels to the current fifth percentile could enhance wheat and maize yields in China by ∼25% and ∼15%, respectively.Several studies have also investigated the sensitivity of crop yield to air pollution at the regional scale, and most are focusing on the impacts of ozone  (Ahmed S, 2015, Lobell andBurney 2021).These global and regional-scale studies consistently highlighted the negative effects of air pollutants on crop yield in the North China Plain (NCP), which is home to 300 million people and is not only an important crop production region in China but also a hotspot of severe air pollution (Wang et al 2006, Yin et al 2019).To date, however, the direct effects of multiple air pollutants via influencing crop physiological processes and indirect effects via regulating regional climate are not well separated for different crop growth stages in the NCP.
To fill the above gaps, we present the first quantitative study on the indirect effects of multiple air pollutants on crop yield for different crop growth stages in the NCP.Besides O 3 , the index of absorbing aerosol optical thickness (AAOD) is used to represent absorbing aerosols, while SO 2 and NO 2 are selected to account for the effects of scattering aerosols, as their oxides account for more than 70% of scattering aerosols in NCP (Yang et al 2015, Zou et al 2018).Specifically, three scientific questions are addressed in this study: (1) what are the differences in the effects of temporal variations of four air pollutants (O 3 , AAOD, SO 2 and NO 2 ) on crop yield during different growth stages across the NCP? (2) how do the four air pollutants exert impacts on crop yield via various influencing paths?(3) what is the relative contribution of the indirect effects of air pollutants to crop yields?

Study area
The NCP (31 Beijing, Tianjin, Hebei Province, Shandong Province and Henan Province.It is a major grain production area in China, accounting for about 27% and 58% of the country's total maize and wheat production by 2020 (www.stats.gov.cn/),respectively.Besides, the NCP is home to about 300 million people and accounts for about 25% of China's total gross domestic product.With rapid growth in the steel, oil refining and petrochemical industries, the NCP has become a hotspot region of severe air pollution where the total exhaust emissions in 2020 account for about 17% of the entire country (figure 1(a)).More importantly, the relatively dry climate in winter makes it difficult to remove air pollutants (Guo et al 2016), and most cities in the region are heated in winter via burning coals which further intensifies air pollution (Ma et al 2011).The overlap of crop production and air pollution makes it a good case study region for quantifying the effects of air pollution on crop yield.

Crop yields and major growing periods
Crop yield census data are obtained from the statistical office website of the province and city (table S1).To be consistent with the time range of available high-resolution air pollutants data, wheat and maize yields data during 2005-2020 are collected for analysis.Following the literature (Li and Wang 2018, Wan et al 2018), the major growth period of winter wheat includes November (seedling and tillering), February and March (greening), April (jointing and tasseling) and May (filling), while the main growth stages of summer maize are June (seedling), July (jointing), August (pregnancy and flowering) and September (filling).

Air pollutants data
Due to data limitation of long-term ground observations of air pollutants during the NCP, we obtain satellite-based datasets of AAOD, total column ozone concentration (O 3 ), total column sulfur dioxide concentration (SO 2 ), and total column NO 2 concentration during 2005-2020 from the ozone monitoring instrument (OMI) sensor product (https://giovanni.gsfc.nasa.gov/),which is jointly developed by the Netherlands, Finland, and the National Aeronautics and Space Administration in the US (Levelt et al 2018).Specifically, we select AAOD data in the 'OMAEROe' dataset based on the near-ultraviolet band algorithm, O 3 data in the 'OMDOAO3e' dataset based on the differential optical absorption spectroscopy algorithm, NO 2 data in the 'OMNO2d' dataset with 30% cloud removal; and SO 2 data in the 'OMSO2e' dataset based on improved principal component analysis (PCA) algorithm.All the above data are in dobson units (Du) except for AAOD, which is dimensionless.The temporal resolution of these satellite data is daily average and spatial resolution is 0.25 • × 0.25 • .
To validate the satellite data, we obtained the ground observation-based near-surface datasets of O 3 , SO 2 , NO 2 from Tang et al (2020), which are only available for the period from 2013 to 2019.By comparing the satellite-based and ground-observed datasets during the overlapping period, we found that the temporal change patterns of total column concentrations of the air pollutants were similar to those of the near-surface conditions, except for O 3 .Therefore, we replaced the total column O 3 dataset with the OMI-MLS tropospheric ozone column concentration dataset (Ziemke et al 2006) for the subsequent analysis, which exhibits a similar change pattern with the nearsurface observation data (figure S1) and has been used to explore the negative effects on crop productivity and terrestrial gross primary productivity (Mahmood et al 2020, Zhu et al 2022).Statistically significant correlations between the temporal variations of OMI and surface observations are found for the O 3 , SO 2 and NO 2 with a high correlation coefficient of 0.91, 0.88 and 0.87, respectively.This suggests that the satellitebased OMI dataset can be well used for our assessment, since we focus on the sensitivity of crop yield variation to the temporal fluctuations in air pollution concentrations.

Climate variables data
For meteorological data, the monthly mean nearsurface daily maximum air temperature (T max , • C) is obtained from the Chinese meteorological station data (http://data.cma.cn/).The monthly mean daily precipitation data (Pre, mm) is from the high-precision (1 km) gridded product (Qu et al 2022), which is developed through interpolation of the Chinese meteorological station data.The direct photosynthetically active radiation (DPar) and diffuse photosynthetically active radiation (DFPar) are derived from Ryu et al (2018), which are available at 5 km spatial resolution at a 4 day interval.These climate variables are selected for our analysis, because they are not only relevant to crop growth, but also respond to the variation of air pollutant concentrations.

Methods
In our study, we synthesize various datasets of crop yields, air pollutants and climate factors and build several novel data-driven models to quantify the indirect effects of air pollution on crop yield via regulating the climate.The overall framework of the methodology is shown in figure 2, and details on the key steps are introduced in the following sub-sections.

Data preprocessing
All the above grid data are reprojected and resampled to 1 km with the software of ArcGIS, while the in-situ gauge temperature data are interpolated into 1 km resolution grid data by kriging interpolation.The gridded data of air pollutants and climate variables are masked firstly by the harvest area map (Luo et al 2020), and then spatially aggregated to the county scale at which census yield data is available.The daily air pollutant data is averaged over each crop growth stage at the monthly scale to derive the annual time series of monthly mean air pollution values.The counties where crop yields data are missing for 10 yr are not included in the research.
Since our study focuses on the effects of temporal variations of air pollutions (rather than trend effects) on crop yield variability, the time series of crop yield, climate and air pollutant variables are all detrended to minimize the trend effects.Instead of using a specific predetermined function such as linear or polynomial regression against time, all data is detrended using the first-difference method (i.e.year-to-year changes), which reflects the change in a variable over two consecutive years.The detrending method can minimize the confounding influences of long-term changes in technology and management (Prabnakorn et al 2018, Wang et al 2022a).Further, the detrended data is standardized to allow comparison of the influences of variables with different units on crop yields.

Sensitivity analysis
We develop a statistical crop model to explore yield sensitivity to the changes in four air pollutants as follows: where ∆A, ∆O, ∆S, ∆N and ∆Y represent the firstorder differential and standardized values of AAOD, O 3 , SO 2 , NO 2 , and crop yields for the county j in the year t, respectively.The first-order differential and standardization methods are adopted to assess the response of yield variability to one standard deviation change in air pollutants (see section 2.3.1).Following previous studies (Liu and Desai 2021, Lobell and Burney 2021, He et al 2022), the dummy variable β j is adopted to account for the fixed effects across counties (wheat: 284 counties; maize: 287 counties).S A , S O , S S and S N are the regression coefficients, which indicate the sensitivity of crop yield to AAOD (A), O 3 (O), SO 2 (S) and NO 2 (N), respectively.For example, the value of S O represents S O standard deviation change in crop yield in response to one standard deviation change in Ozone.ϵ jt is the error term.For each crop growth stage (wheat: November, February, March, April and May; maize: June, July, August and September), a fixed effects model is constructed to understand the intra-seasonal differences in the effects of air pollution.
To address the effect of multicollinearity between air pollutants (e.g.Ozone-NO 2 interactions), we employ the partial least squares regression (PLSR) to fit the fixed effects model, which combines the advantages of PCA, typical correlation analysis, and linear regression without requiring normal assumptions on data distribution (Geladi and Kowalski 1986).In this study, the PLSR is performed with the 'pls' package in R to obtain the sensitivity of crop yields to four air pollutants for each growth stage (i.e. S A , S O , S S and S N ).

Path analysis
To investigate the indirect effects of air pollutions, path analysis which is based on a path structure consisting of variables with potential connectivity (Li 1975) is adopted.It has been widely used in previous studies to explore the relationship between multiple variables (Henseler and Sarstedt 2013, Kimm et al 2020, Yu and Leng 2022).In our path analysis, it is assumed that air pollutant affects crop yields indirectly by influencing regional climate (i.e.air pollutants → climate factors → crop yields) and the path coefficients are derived from the regression coefficients obtained by PLSR.The indirect impacts of AAOD, O 3 , SO 2 and NO 2 refer to the sensitivity of crop yield to air pollutants via regulating climate conditions and can be displayed as P A , P O , P S and P N , respectively.Besides, we simplify the impacts of different paths as W(M) air pollutant ± climate factor , where 'W' and 'M' represent wheat yield and Maize yield, respectively; 'air pollutant' represents AAOD (A), O 3 (O), SO 2 (S), or NO 2 (N); 'climate factor' represents T max (T), DPar (DP), Pre (P) or DFPar (DFP); and the symbol of '±' means that air pollutant affects crop yields positively (+) or negatively (-) via influencing climate factor.

Random forest model (RFM)
The RFM algorithm is based on the algorithm of bagging by generating several regression trees, based on which the mean value of all trees is adopted as the prediction result (Breiman 2001).In this study, the RFM is performed with the 'randomForest' package in R and the relative contribution of each explanatory variable is calculated by 'IncNodePurity' .The construction of RFM in this study is based on the abovementioned influencing paths.First, the interannual variations in each climate factor are simulated by the four air pollutant variables [Model 1: air pollutant indicators → climate factor].Second, the interannual variations in crop yield are simulated by the four climate variables [Model 2: climate factors → crop yield].The relative contribution indicates the portion of yield variability that can be explained by air pollution via regulating climate conditions.We display the relative contribution of each explanatory variable in different paths as RC A-B which means that variable A contributes to the variation of variable B.

Temporal and spatial changes in crop yields and air pollutant indicators
Figures 3(a)-(f) show the temporal changes in wheat and maize yields and the seasonal variation of air pollutants in the NCP during 2005-2020.Besides, the temporal changes of air pollutants and climate factors are exhibited in figure S2.Overall, the median yields of wheat and maize vary within a slight range of 5−7 t ha −1 with an upward trend of 0.05 t ha −1 a −1 (p < 0.01) and 0.02 t ha −1 a −1 (p > 0.1), respectively (figures 3(a) and (b)).High levels of AAOD (0.06−0.08),SO 2 (0.6−0.8 Du), O 3 (42-52 Du) and NO 2 (0.6-1.0 Du) are found for the major growing seasons of wheat and maize (figures 3(c)-(f)), and evident correlations between yields and air pollutants are observed in much of NCP (figures S3-S4).Spatially, the southern Hebei province, northern Henan province and western Shandong province are hotspot areas with high crop yields ranging of 6.00−8.00t ha −1 and serious air pollution (figures 3(g)-(l)).

Indirect effects of air pollutants on crop yields during different growth stages
The indirect effects of air pollutants on wheat and maize yields via various paths are shown in figures 5 and 6.Overall, the dominant indirect influences of air pollutants on yields are dependent on the specific growth stage.Specifically, in November (figure 5

Relative contributions of each explanatory variable in different paths
The dominant paths that regulate crop yield variation are obtained via path analysis, and the relative contributions of each explanatory in dominant paths gained by RFM are exhibited in (figures 7(a)-(h)).Specifically, AAOD exerts an indirect effect which contributes to 53% of wheat yield variation by regulating 58% of T max variation (RC T-W = 16%), 65% of DPar variation (RC DP-W = 34%), 35% of Pre variation (RC P-W = 33%) and 58% of DFPar variation (RC DFP-W = 17%).O 3 generates an indirect effect which contributes to 24% of wheat yield variation by regulating 18% of T max variation (RC T-W = 16%), 14% of DPar variation (RC DP-W = 34%), 36% of Pre variation (RC P-W = 33%) and 25% of DFPar variation (RC DFP-W = 17%) in November (figure 7(a)).In February and March, 32% of wheat yield variation is indirectly explained by AAOD which modulates 26% of T max variation (RC T-W = 22%),

Dominant roles of indirect effects of air pollutants on crop yields during different growth stages
Taking the NCP as a whole, the positive indirect effect of AAOD and the negative indirect effects of O 3 and SO 2 dominate wheat yield variations in November with contributions even exceeding the direct effects (figure 8(a)).Spatially, the indirect effect of

Indirect effects of air pollutants on crop yields
In our study, we use statistical models to assess the indirect effects of air pollutants on crop yield, since a fully coupled process model with suitable representations of air pollution-climate-crop interactions is not readily available.We build the statistical model based upon our physical knowledges of air pollutionclimate-crop interactions from the literature.
Overall, the effects of air pollutants on climate show large spatiotemporal variations (tables S4-S5 and figures S7-S10), which are broadly consistent with previous studies (Tao et al 2012, Bellouin et al 2020).Specifically, in most crop-growing areas of NCP, an increase in AAOD could give rise to a decline  8(a)).Our results further show that such indirect effects of air pollutants on crop yield could exceed 50% of its net effects in the NCP, though the relative contributions depend on the crop growth stage.The revealed patterns of the indirect effects of air pollutants have great implications for more targeted adaptation and mitigation efforts.Indeed, for counties where the indirect effects of air pollutants on crop yields are dominant, it implies that measures addressing the heat and water stresses induced by air pollution should be deployed.In contrast, in counties where the indirect effects of air pollutants are smaller than the direct effects, it suggests that new breeding or molecular biology techniques are required to improve the tolerance of crops to air pollutants to resist the disruption of their physiological processes (Teixeira et al 2011).
Besides the main paths examined in this study, air pollutants could exert impacts on crop yields via other mechanisms.For instance, although the particles of AAOD attached to crop leaves could negatively affect the function of leaves (Burkhardt and Pariyar 2014), Wang et al (2002) reported that dust aerosol particle deposition was beneficial to crop growth by contributing to the alkaline neutralization of acidic soils.In addition, the S-oxide and N-oxide in aerosols as the form of acid rain could further reduce crop yield by influencing soil properties (Wang et al 2018, Du et al 2020).Thus, the mechanisms underlying the impacts of air pollutants on crop yields are more complex and need to be further investigated in the future.

Differences in the sensitivity of wheat and maize yields to air pollutants
In this study, it is shown that wheat is more sensitive to O 3 than maize in the NCP, which is consistent with previous studies indicating larger wheat yield loss than maize under high O 3 concentrations (Avnery et al 2011b, Feng et al 2022).The different sensitivity of the two crops to air pollutants is related to their physiological characteristics and differences in light, temperature, and water requirements.Indeed, the C 4 crop (e.g.maize) is more photosynthetic and more tolerant to heat and drought than the C 3 crop (e.g.wheat) (Lychuk et al 2017, Kukal and Irmak 2020, Ozeki et al 2022).Meanwhile, the concentrations of AAOD, SO 2 and NO 2 are higher during the growing season of wheat, which may lead to a stronger negative sensitivity of wheat yield to the three air pollutants.
Previous studies showed that ozone often causes large crop yield losses, but such conclusions depend on the data and methods.For example, Lobell and Burney (2021) found that PM and NO 2 appear more damaging to crops than O 3 and SO 2 .He et al (2022) showed larger wheat yield loss in response to 0.1 AOD increase than 5ppbv ozone increase.Liu and Desai (2021) showed Ozone and aerosol reduced maize yield by 15.90% (90% CI: 14.44%-17.53%)and 18.88% (90% CI: 10.97%-27.04%).Indeed, there exists uncertainty in comparing the sensitivity of crops to different air pollutants with different units.In our study, we standardize all data to eliminate the effect of different units before our assessment.It is shown that obtained NO 2 exerts the largest negative effects on wheat yield in November, April, May, and AAOD generates the largest adverse effect on wheat yield in February and March, while NO 2 exerts the largest positive effects on maize yield in June and September, O 3 generates the largest adverse effect on maize yield in July, and AAOD produces the largest favorable effect on maize yield in August.
The different effects of air pollutants on climate factors combined with the varying requirements of crops for light, heat and water in each crop growing stage lead to distinct spatiotemporal patterns of the indirect effects of air pollutants on crop yields (figures S7-S13).Generally, crops are more vulnerable to drought and light deficiency during the reproductive stage (from tasseling to filling) with stronger transpiration and photosynthesis (Butler and Huybers 2015, Feng et al 2019, Hu et al 2019).Our results confirm this by showing that air pollutants exert a larger indirect influence on crop yields during this stage (figures S12-S13).In addition, the direct effects of different air pollutants on crop yields also vary among crop growth stages.We can see that O 3 has a negative direct effect on maize yield in July but a positive direct influence in June, August and September (figure 8(b)), which might be due to the differences in stomatal conductance and other agronomic managements at different growth stages which subsequently affect the uptake of gaseous air pollutants (Mills et al 2018b).Overall, the combination of indirect impacts of air pollutants by affecting climate factors and direct impacts via influencing crop physiology would exacerbate the spatiotemporal heterogeneity of net effects on crop yields in each growth stage (figures S5-S6).
Many studies have assessed the effects of air pollutants on crop yields averaged over the whole growing season, which may overestimate or underestimate the effects of air pollutants in a given growth stage.For instance, Tie et al (2016) demonstrated that aerosols have caused a significant yield loss of wheat and rice by negatively affecting growing-season mean radiation.Our study also found that AAOD exerts negative effects on wheat yield when considering the growing season mean, but we observed that the effects of AAOD depend on the growth stage (figure S14(a)).Specifically, AAOD generates a negative impact in February and March (greening), but positive influences in November (seeding and tillering), April (jointing and tasseling) and May (filling) (figure S14(a)).In addition, although O 3 has a negative effect on wheat yield when looking at the growing season mean, a positive effect is found in April as its positive indirect effect offsets its adverse direct effect (figure S14(a)).

Limitations of our study
Several limitations of this study should be acknowledged for better interpreting the results obtained from our assessment.First, we used satellite remote sensing products to obtain the long-term time series and high-resolution monthly mean air pollutants data, due to the fact that ground stations covering major cities in China were established at the end of 2012 and long-term time series of ground observations of hourly air pollutants are not available.Unlike the cumulative indicators such as AOT40 of ozone used in previous studies, our results may underestimate the effect of air pollution on crops as the monthly mean air pollutants indicators smooth out the extremely high air pollution concentration.
Second, three types of gaseous pollutants used in this study (O 3 , SO 2 , NO 2 ) mainly reflect the column concentrations, which may not be accurate to characterize pollutant concentrations at the height of crop canopy.We acknowledge the potential differences between total-column and near-surface air pollutant concentrations, but here we mainly focus on the longterm effects of air pollutants variations on crop yield and similar temporal change patterns between total column and near-surface concentrations are found for O 3 , SO 2 and NO 2 for the period 2013-2019.
Third, it is well known that oxides of nitrogen (NO x ) and volatile organic compounds are important precursors for ozone production and an increase in NO x would induce an increase in O 3 since O 3 formation is generally NO x limited (Duncan et al 2010).Meanwhile, the strong oxidation of O 3 will influence the formation process of aerosol particulate matter (Shang et al 2021), which could in turn counteract O 3 production by decreasing the amount of UV radiation reaching the earth's surface (Liu et al 2019).These complex interactions between the four air pollutants lead to multicollinearity of variables in interannual variation, which could be well addressed by the PLSR (Kemalbay and Korkmazoglu 2012).
Fourth, the impacts of air pollutants are assessed without explicitly considering the effects of adaptations such as irrigation, choices of cultivar, sowing data and extreme weather effects at daily and subdaily scales.We examine the main paths through which air pollutants influence crops via regulating T max , DPar, Pre and DFPar by path analysis, which also has limitations by assuming linear relationships among variables (Jeon 2015).Other important mechanisms are not examined explicitly, and process crop models can be used in this regard to complement with the separation of direct and indirect effects of air pollutants on crop yields.
Therefore, we have to acknowledge that this study is not intended to give an accurate estimate of the effects of air pollution on crop yield.Instead, we provide a sensitivity analysis to measure the potential indirect effects of air pollution on crop yields via regulating regional climate in the NCP, which is a hotspot region in terms of air pollution and crop production in China.Future efforts should be made by collecting more ground observations and revisiting the assessment using process crop models together with datadriven models.

Conclusions
High concentrations of air pollution can not only affect crop yield directly but also indirectly via modulating regional climate variations.How crop yields respond to multiple air pollutants during different crop growth stages is not well examined in the NCP, which is an important area for crop production in China but facing severe air pollution.In this study, we present the first quantification of the indirect effects of air pollutants on crop yields in the NCP, and the main findings of this work are as follows: (2) Based on path analysis and RFM, we found that AAOD and O 3 exert the larger indirect impacts on wheat yield in November (P A : 0. (3) Further analysis shows that the indirect effects of AAOD, O 3 and SO 2 on wheat yield variation can exceed the direct effects in November.Similar findings are obtained for O 3 which exerts a larger positive indirect effect than the negative direct effects in April.In June, maize yield variation is mainly dominated by the indirect effects of AAOD and O 3 , both of which could surpass their direct effects.The negative indirect impact of O 3 and positive indirect influence of SO 2 are significant in July, which dominate maize yield variations.In August, maize yield variation is controlled by the negative indirect effect of O 3 and the positive indirect effect of NO 2 .The negative indirect impact of O 3 is also significant in September with contributions surpassing its positive direct effects.
This study highlights the important but leastunderstood roles of the indirect effects of air pollutants on crop yield during different crop growth stages, which has great implications for targeted adaptation and mitigation efforts toward better addressing the challenges of air pollution and food security.While acknowledging the uncertainties from satellite datasets and statistical models, we call for more efforts to be made in the future to collect long-term hourly ground observations of air pollutions and revisit the assessment by developing and using a fully coupled process model with detailed representation of air pollution-climate-crop interactions.

Figure 1 .
Figure 1.Map of the study area.The spatial patterns of total exhaust emissions over China in 2020 (a) (https://data.stats.gov.cn/); the elevation of the study area (b) (www.gscloud.cn/);the percentages of production and planting area of wheat and maize in the study region compared to the whole China (c) (https://data.stats.gov.cn/); the spatial maps of long-term mean climate factors (Tmax: daily maximum air temperature; DPar: direct photosynthetically active radiation; Pre: daily precipitation; DFPar: diffuse photosynthetically active radiation) during 2005-2020 (d) (data information can be seen in section 2.2.3).

Figure 2 .
Figure 2. Schematic overview of quantitative analysis conducted in this study.

Figure 3 .
Figure 3. Spatial and temporal distribution patterns of crop yields and air pollutants in the NCP.(a) and (b) are the temporal changes in wheat and maize yields, respectively; (c)−(f) are the seasonal variation of air pollutants; (g)−(l) are the spatial patterns of long-term mean crop yields and air pollutants during 2005-2020, respectively.The red line in (a and b) denotes the change trend of yields by fitting the median of crop yields data across all counties against time.The horizontal line, hollow square, edges and whiskers of the boxplots indicate the mean, median, interquartile range and 5th-95th percentiles of yields across all counties, respectively.In (g) and (h), the diagonals represent areas where data are not available.

4.
The sensitivity of wheat (a) and maize yields (b) to air pollutants at the regional level during different crop growing stages obtained by PLSR.The ' * ' indicates the regression coefficient passes the 90% significance test.

Figure 5 .
Figure 5.The paths through which air pollutants exert indirect effects on wheat yield during different wheat growing stages (a)-(d).The path coefficients (ρ: −1-1) are shown along with the path arrows, where red (blue) indicates positive (negative) effects, solid (dashed) indicates statistically significant (insignificant) effects, and the thickness of the line indicates the magnitude of coefficients.

Figure 6 .
Figure 6.The same as figure 5 but for maize yield.

Figure 7 .
Figure 7.The relative contributions of air pollutants to climate variation and its feedback effects on wheat yield (a)-(d) and maize yield (e)-(h) during different crop growing stages.The outer circle indicates the relative contributions of four air pollutants to the variations of each climate variable, while the inner circle donates the relative contribution of each climate variable to the variations of crop yield.

Figure 8 .
Figure 8.The net effects (colored bar), indirect effects (masked bar) and direct effects (red box) of four air pollutants on wheat (a) and maize yields (b) in the NCP during different growing stages.The net effects include the direct effects of air pollutants on crop yields and their indirect effects via regulating climate variations.
(O 3 ) (Feng et al 2015, Tian et al 2016, Dong et al 2021, Li et al 2022a, Wang et al 2022b) followed by aerosols (Zhang et al 2017, Zhang et al 2018b, He et al 2022), sulfur dioxide (SO 2 ) and NO 2 041, P O : −0.051) with the RC A-W of 53% and RC O-W of 24%.AAOD and SO 2 caused the larger indirect effects on wheat yield in February and March (P A : −0.008, P S : 0.018) with the RC A-W of 32% and RC S-W of 25%, respectively.In April, the indirect effects of O 3 and NO 2 on wheat yield are larger (P O : 0.034, P N : −0.012) with the RC O-W of 26% and RC N-W of 25%.In May, the indirect effects of AAOD and NO 2 (P A : 0.022, P N : −0.029) on wheat yield contribute 26% and 18% of wheat yield variation, respectively.For maize in June, positive indirect effects of AAOD (P A : 0.013) and negative indirect effects of O 3 (P O : −0.017), were observed with the RC A-M of 22% and RC O-M of 28%.The indirect effects of O 3 and SO 2 are larger both in July (P O : −0.045, P S : 0.018) with the RC O-M of 33% and RC S-M of 16%, and in August (P O : −0.017, P S : 0.009) with the RC O-M of 21% and RC S-M of 24%.In September, AAOD and O 3 exert the larger indirect impacts on maize yield (P A : 0.013, P O : −0.024) with the RC A-M of 21% and RC O-M of 42%.