Has the sensitivity of soybean cultivars to ozone pollution increased with time? An analysis of published dose–response data

The rising trend in concentrations of ground‐level ozone (O3) – a common air pollutant and phytotoxin – currently being experienced in some world regions represents a threat to agricultural yield. Soybean (Glycine max (L.) Merr.) is an O3‐sensitive crop species and is experiencing increasing global demand as a dietary protein source and constituent of livestock feed. In this study, we collate O3 exposure‐yield data for 49 soybean cultivars, from 28 experimental studies published between 1982 and 2014, to produce an updated dose–response function for soybean. Different cultivars were seen to vary considerably in their sensitivity to O3, with estimated yield loss due to O3 ranging from 13.3% for the least sensitive cultivar to 37.9% for the most sensitive, at a 7‐h mean O3 concentration (M7) of 55 ppb – a level frequently observed in regions of the USA, India and China in recent years. The year of cultivar release, country of data collection and type of O3 exposure used were all important explanatory variables in a multivariate regression model describing soybean yield response to O3. The data show that the O3 sensitivity of soybean cultivars increased by an average of 32.5% between 1960 and 2000, suggesting that selective breeding strategies targeting high yield and high stomatal conductance may have inadvertently selected for greater O3 sensitivity over time. Higher sensitivity was observed in data from India and China compared to the USA, although it is difficult to determine whether this effect is the result of differential cultivar physiology, or related to local environmental factors such as co‐occurring pollutants. Gaining further understanding of the underlying mechanisms that govern the sensitivity of soybean cultivars to O3 will be important in shaping future strategies for breeding O3‐tolerant cultivars.


Introduction
Ensuring that the rising global population has access to a sufficient and stable food supply is a key international priority for the 21st century. At a time when an estimated 795 million people worldwide are undernourished (FAO, 2015), agricultural productivity is being limited by several factors, including inter alia, rising water scarcity (Falkenmark, 2013), the limited land available for cultivation (Zabel et al., 2014), widespread soil erosion and degradation (FAO, 2011), and the impacts of climate change (Parry et al., 2004). A further threat to agricultural yield comes from rising concentrations of ground-level ozone (O 3 ) (Fuhrer, 2009)a common air pollutant and phytotoxin (Krupa et al., 2001). Ozone is a secondary pollutant formed in photochemical reactions from precursor compounds, the most important of which are nitrogen oxides (NO x ), methane (CH 4 ) and carbon monoxide (CO) (Royal Society, 2008). The global surface background concentration of O 3 more than doubled between the early 1900s and the end of the 20th century (Hough & Derwent, 1990;Parrish et al., 2014), most likely as a result of rising anthropogenic emissions of O 3 precursor compounds from fossil fuel combustion, biomass burning and paddy field cultivation (Brasseur et al., 2001). Projected changes in global surface O 3 for the period 2000-2050 range from a decrease in the 24-h mean of 2.5-7.2 ppb under the optimistic emission scenarios (RCP2.6, RCP4.5, RCP6.0, B1) to an increase in 1.5-6.2 ppb under the more pessimistic RCP8.5 and A2 emission scenarios (IPCC, 2013). Trends in surface O 3 are, however, highly variable geographically, and the most rapid increase is currently occurring in South Asia where surface O 3 concentrations are expected to continue to rise until 2050 under all but one of the emission scenarios (Beig & Singh, 2007;IPCC, 2013). Establishing a thorough understanding of crop and cultivar responses to O 3 , and the incorporation of these responses into crop production models, is therefore needed to quantify the potential impact of O 3 on food supply in different world regions.
Soybean (Glycine max. (L.) Merr.) ranks among the most O 3 -sensitive agricultural crops (Mills et al., 2007). It is the fifth most significant crop in terms of global production (FAO, 2012), is a key source of vegetable protein for humans (Mateos-Aparicio et al., 2008), provides approximately 30% of the world's processed vegetable oil (Graham & Vance, 2003), accounts for 77% of global nitrogen fixation by crop legumes (Herridge et al., 2008) and is an important feed constituent for the livestock and aquaculture industries (Hartman et al., 2011). The crop holds significant economic importance for a number of world economies including the USA, Brazil, Argentina, China and India (FAO, 2014), and world soybean demand is increasing by an average of 2.2% annually (Masuda & Goldsmith, 2009). Ozone exposure reduces the photosynthetic rate, stomatal conductance (g s ), leaf chlorophyll content and leaf starch concentration of soybean (Morgan et al., 2003). Groundlevel O 3 pollution over agricultural land has been estimated to cause an annual reduction in soybean yield ranging between 6 and 16%, and financial losses of $2.0-5.8 billion annually, based on analysis of year 2000 data conducted in two separate global crop loss assessments (Van Dingenen et al., 2009;Avnery et al., 2011a). Soybean crop yield reduction for the year 2030 as a result of O 3 is estimated to be 9.5-15% under the optimistic (B1) scenario or 15-19% under the pessimistic (A2) emission scenario (Avnery et al., 2011b).
The magnitude of O 3 damage to soybean is dependent on the timing of exposure, with greater reductions in photosynthesis and yield being observed when exposure occurs during the reproductive stages of growth (Morgan et al., 2003). Co-occurrence of seasonal peaks in O 3 surface concentrations and the flowering and pod-filling stages could therefore be particularly damaging for yield. Ozone damage occurs when the gaseous pollutant enters the leaf via the stomatal pores, and interacts with cell membranes and walls in the apoplast to yield reactive oxygen species (ROS) (Wilkinson et al., 2012); these directly damage plant tissue through protein oxidation, leading to accelerated senescence and cell death . The widely observed reduction in photosynthetic rate in response to O 3 is not fully understood, but is in part the result of a reduction in the leaf concentrations of chlorophyll and Rubisco (Glick et al., 1995;Fiscus et al., 2005). Ozone has also been observed to reduce nodulation in a range of legume species including soybean (Tingey & Blum, 1973;Reinert & Weber, 1980;Zhao et al., 2012), although this effect is largely thought to be a secondary response as a result of reduced total carbon assimilation and the diversion of assimilates away from the roots (Hewitt et al., 2016).
Dose-response studies for a range of crops have revealed that O 3 sensitivity is a heritable trait (Reinert & Eason, 2000) and is highly variable among species and among cultivars (Ariyaphanphitak et al., 2005;Mills et al., 2007;Mills & Harmens, 2011). The maximum stomatal conductance which a species or cultivar can reach (g max ) is thought to play a role in determining O 3 sensitivity, because greater conductance results in greater O 3 uptake. This view is supported by the observation that wheat cultivar sensitivity to O 3 is positively correlated with g max (Biswas et al., 2008). Furthermore, modern wheat varieties are more sensitive to O 3 than older varieties; this may be a result of selective breeding programmes targeting varieties with a higher g s , as these have a higher rate of CO 2 fixation leading to higher yields (Biswas et al., 2008;Roche, 2015). The detoxification and repair capacity of a plant species or variety is also thought to be important in determining sensitivity : for example, O 3 tolerance of a number of plant species has been seen to positively correlate with greater apoplastic concentrations of ascorbic acid, an antioxidant (Frei et al., 2008(Frei et al., , 2010Frei, 2015). A thorough understanding of how O 3 sensitivity varies among cultivars of the same speciesand the factors which drive these differencesis key in improving assessments of current and future O 3 -induced crop losses. Previous studies in soybean investigating intercultivar variation in O 3 response have typically compared a relatively small number of cultivars from the same geographical region: examples include studies of USA cultivars by Betzelberger et al. (2010Betzelberger et al. ( , 2012 and an investigation of Chinese cultivars by Zhang et al. (2014). Knowledge of which cultivars are most resistant to the effects of O 3 could potentially help plant breeders to develop O 3 -tolerant soybean varieties, which, if adopted by farmers, could mitigate O 3 -induced crop losses.
Much of the research relating to soybean-O 3 responses conducted to date has taken place in the USA, as part of the US National Crop Loss Assessment Network (NCLAN) programme in the 1970s and 1980s (Heagle, 1989) and more recently at the University of Illinois, Urbana-Champaign and USDA Agricultural Research Service SoyFACE facility in Illinois (Long et al., 2005;Betzelberger et al., 2010Betzelberger et al., , 2012Gillespie et al., 2012). Groups in India and China have also studied O 3 responses to soybean in recent years, but these data have, to date, not been pooled to produce dose-response relationships. Response functions for soybean used in global crop loss assessments have therefore been based on experimental data collected only in the USA. Two dose-response functions for soybean have been published: one by Lesser et al. (1990), synthesized from the NCLAN data set; and one by Mills et al. (2007), who combined some of the NCLAN data with more recent dose-response data collected in the USA to update the function. These functions have been applied in a number of different studies in order to estimate O 3 -induced soybean yield reduction globally and the associated financial loss to farmers. Producing these estimates involves combining a dose-response function for soybean with crop distribution and yield maps, growing season dates and modelled O 3 concentrations. The Mills et al. (2007) function was used by Avnery et al. (2011a) in their global assessment of O 3 -induced soybean crop losses . The Lesser et al. (1990) function was used by Wang & Mauzerall (2004) in their soybean yield loss assessment for East Asia and by Van Dingenen et al. (2009) in their global assessment. Both functions were used by Tai et al. (2014) in their analysis of combined O 3 and climate change effects on future soybean production. All of these assessments applied a soybean dose-response function based on data from North America to model yield impacts in Asia. However, a comparison by Emberson et al. (2009) of wheat and rice dose-response data from North America and Asia has shown that Asian wheat and rice cultivars appear to be more sensitive to O 3 than their North American counterparts, possibly due to locally occurring physiological traits associated with sensitivity, such as high g s and low antioxidative capacity . The application of North American dose-response functions in global yield loss assessments for wheat and rice may have therefore underestimated O 3 -induced yield losses in Asia.
This study, to our knowledge, synthesizes all existing data in the scientific literature describing soybean yield response to O 3 , in order to produce a comprehensive and up-to-date dose-response function. We also analyse intercultivar differences in O 3 sensitivity, allowing the most O 3 -sensitive and O 3 -tolerant soybean cultivars to be identified. Additional analysis is also conducted on the dose-response data set, to investigate potential correlations between the degree of O 3 sensitivity observed and (i) the year in which the soybean cultivar was released, to identify temporal trends in sensitivity; (ii) the geographical location of the dose-response experiment, to determine whether sensitivity varies geographically; and (iii) the method of O 3 fumigation used in experimentation, to assess whether experimental design influences the sensitivity observed.

Literature search
A search of the published scientific literature was performed between October 2013 and September 2014 to find all O 3 exposure studies conducted on soybean. The search was conducted using the Science Citation Index Expanded â (Thomson-ISI, Philadelphia, PA, USA). The criteria for inclusion were as follows:  ersen & Lauer, 2004). About 60% of this period is equal to 7.7 weeks, which was rounded to a minimum exposure duration of 8 weeks for the purpose of this study. 4. Yield must have been measured directly, as the pod or seed weight. Response parameters such as total aboveground biomass, photosynthetic rate, percentage leaf damage or the 100-seed weight were not considered to represent the yield response.
The literature search found 28 studies that met the inclusion criteria and were included in the analysis. These studies included experiments investigating 48 cultivars, and when combined, produced a data set comprising 379 data points. A list of all the experimental studies included in our analysis can be found in Table 1, alongside information relating to study sites, cultivars tested and experimental design. Experiments which had used pot-grown soybean were included in our analysis; this was justified given that we found no significant difference in the dose-response relationships exhibited by pot-grown and field-grown soybean (See Appendix S1).

Standardization of O 3 and yield parameters
Dose-response data in the literature were presented using a number of different concentration metrics and yield parameters, as listed above. All O 3 concentration data had to be converted into a standard metric to enable the data to be combined for analysis. The M7 was selected to act as the common O 3 metric in our analysis, because this was most frequently reported in the literature. O 3 values presented in the form of the AOT40, M12 and M24 were converted to the M7 metric using conversion functions calculated using The ICP Vegetation database (described in Fig. 1), which contains O 3 observations measured at the same time and location but using a range of different O 3 metrics. The three different conversion functions which we used to standardize our data to the M7 metric, calculated from the ICP Vegetation database, are shown in Fig. 1 represented using the different O 3 metrics were then plotted against each other, and conversion functions were derived using linear regression. During standardization of the reported O 3 concentrations to the M7, concentrations presented as the M8 were considered to be equivalent to the M7, as the small difference between the two was considered unlikely to add significant uncertainty to our analysis. A total of 205 O 3 concentration values were presented in the soybean data set using the M7 or M8 metrics and did not need to be converted. A total of 125 and 49 data points were presented using the AOT40 and M12 metrics, respectively, and were converted to M7. Not all of the studies included in our analysis used a full 3-month O 3 exposure; for studies which had shorter exposure durations, it was assumed that the 3-month mean would not radically differ from the mean covering a shorter duration, as O 3 exposure in all studies was artificial and therefore would not follow natural seasonal patterns in O 3 concentration. No study which had used an exposure duration of less than 8 weeks (60% of the soybean growing season) was included in our analysis. Only one study -Betzelberger et al. (2012) required conversion of the AOT40 to the M7, and this study used an exposure duration of 3 months. The process of conversion to the M7 metric had the potential to introduce some error into our data set, which we tested for in our statistical analysis.
As with O 3 concentration, yield was reported in the literature using a range of different metrics, and the control O 3 concentration varied considerably between the different experiments. Yield data were therefore standardized following the method described by Fuhrer et al. (1997). For each separate O 3 exposure experiment, linear regression was used to determine the theoretical yield at 0 ppb O 3 , expressed as the M7 metric. In a second step, the theoretical yield at 0 ppb O 3 was used as the reference (relative yield of 1) for calculating relative yields. The range of theoretical yields at 0 ppb O 3 for each study is included in Table 1.

Derivation of species and cultivar dose-response functions
All statistical analyses were carried out using R software (R Development Core Team, 2015). To calculate the overall doseresponse function for soybean, relative yield data from all studies which met the inclusion criteria were pooled and plotted against the seasonal M7. The shape of the distribution was determined by fitting linear, quadratic and Weibull functions to the combined dose-response data set. Goodness of fit of the model best-fit lines was compared by eye and using the Akaike information criterion (AIC). The linear model was found to be the best fit to the data (AIC values for linear, quadratic and Weibull models are reported in the results section). Linear modelling was therefore chosen as the method to be used in the derivation of independent dose-response functions for individual soybean cultivars which had three or more supporting data points. A mixed model was used when deriving the overall dose-response function for soybean, and in the derivation of individual cultivar dose-response functions, with experimental study included as a random effect to account for the nonindependence of data points originating from the same study. During model fitting, the intercept was allowed to vary and was not forced through a relative yield value of 1. This decision was made to better allow for comparisons of the O 3 sensitivity of the different soybean cultivars based on their dose-response functions. Allowing the intercept to vary around 1 did not result in any systematic bias in the calculated slopes of the dose-response functions (See Appendix S2). Analysis of the effect of cultivar release date, country of study and fumigation method on O 3 sensitivity Stepwise model selection was used to determine whether the cultivar release date, country of data collection and method of O 3 fumigation were important explanatory variables in the model describing the response of soybean to O 3 . A fourth explanatory variable describing whether the O 3 concentration values had been reported as the M7 or had been converted was also included, to test for bias in the data introduced through standardization to the M7 metric. A mixed-effect model structure was used to allow experimental study to act as a random effect. Model fit was assessed using Akaike's information criterion (AIC), a goodness-of-fit parameter calculated from the number of fitted parameters in a model and the maximum likelihood estimate (Symonds & Moussalli, 2011). Cultivar release dates were taken from Specht & Williams (1984), the USDA Germplasm database (USDA, 2015) and the Indian Council of Agricultural Research Oilseed report (ICAR, 2012). Data transformation of the response variable (relative yield) was carried out before analysis by taking the base-10 logarithm, to correct for non-normality observed in model residuals.
Before beginning the analysis, a diagnostic test was carried out on the data set to test the degree of collinearity between the explanatory variables. The presence of collinearity can be a concern in multiple regression due to difficulties differentiating the separate influence of variables that are partially correlated with each other (Belsley et al., 1980). The variance inflation factor (VIF), a widely used measure of the degree collinearity of independent variables in a regression model (O'brien, 2007), was calculated for each explanatory variable (See Appendix S3 and Table S1). Calculated VIF values ranged from 1.1 to 6.1, falling well below the value of 10 considered to be a threshold above which it is recommended that measures are taken to counter the effects of collinearity (Mason & Perreault, 1991;Smith et al., 2009). The diagnostic test, however, reveals the presence of a certain degree of collinearity in the data, meaning that we cannot with complete certainty rank the explanatory variables in order of their relative importance. Nevertheless, we are able to identify which of the candidate explanatory variables are likely to be important in describing the dose-response of soybean to O 3 .
Multivariate regression analysis was stepwise and began with the simplest model (yield~O 3 ), with variables sequentially added to create a more complex model, and goodnessof-fit assessment at each step to determine whether variables should be kept or removed. The order of variable addition was determined by adding each explanatory variable individually to the simplest model, to identify the single variable which gave the greatest improvement to model fit; this model was then carried forward and the process was repeated until the best model was found. A complete list of all the model configurations tested during stepwise selection is given in Appendix S4 and Table S2 of the supporting information.
Candidate explanatory variables which were present in the 'best' model describing the response of soybean yield to O 3 were investigated further by subsequent graph plotting and separate individual regression analyses, which also used a mixed model structure.
Linear regression to determine how soybean cultivar sensitivity has changed with year of cultivar release Soybean cultivars represented in the data set by three or more data points (25 cultivars in total, 22 tested in USA and 3 tested in Indialisted in Table 2) were included in a separate linear regression analysis to determine whether cultivar sensitivity (represented by the slope of the dose-response function) was related to the year of cultivar release. The regression analysis was carried out twice, once on all cultivars and once excluding the cultivars from India, to ensure that any geographical differences in sensitivity were not biasing the observed relationship between sensitivity and year of release.

Reporting yield reductions predicted by dose-response functions
The standardization of reported yield data from the literature was achieved by scaling all data to yield at 0 ppb O 3 . However, when reporting the yield reductions predicted by our dose-response functions in the results and discussion sections of this manuscript, we reasoned that it would be more useful to express yield reductions relative to the naturally occurring background O 3 concentration. Yield reduction estimates presented in the results and discussion of this manuscript have therefore been calculated relative to pre-industrial O 3 levels in Europe, which are thought to have averaged around 20 ppb M24, or 23 ppb M7 (Vingarzan, 2004). The O 3 concentration used to represent present-day background levels was 55 ppb M7a background concentration which has been commonly exceeded in the last 20 years across different world regions (Jaffe & Ray, 2007;Wang et al., 2007;Chakraborty et al., 2015). Relative yield reduction at the present-day O 3 concentration relative to the pre-industrial concentration will hereafter be referred to as RYL c,p in this paper. A graphical representation of the method used to calculate RYL c,p is shown in Fig. 2.

Results
The overall yield response of soybean to O 3 , combined across all cultivars, regions and exposure types, is shown in Fig. 3a. Fitting quadratic and Weibull functions to the data set did not improve the model goodness of fit, suggesting that the soybean response to O 3 was linear across the range of M7 index values examined here (linear AIC = À458, quadratic AIC = À456, Weibull AIC = À453). The combined soybean doseresponse function in Fig. 3, calculated using a mixedeffect model, estimates a RYL c,p of 17.3%. For comparison with earlier studies, the response function for the same data set but using AOT40 as the O 3 metric is provided in Fig. 3b.
Of the 49 cultivars reported in the literature, 25 had three or more data points supporting their doseresponse relationship and therefore were analysed independently using linear regression. The dose-   (1986, 1991, 1983a), Heagle & Letchworth (1982), Heagle et al. (1983bHeagle et al. ( , 1987 Dwight Y = À0.0049x + Williams-79 Y = À0.0047x +  Table S2 of the supporting information. The model that performed best in describing the response of soybean to O 3 included the year of cultivar release, country of study and type of O 3 exposure as interacting variables. The AIC value for the best model shows a far greater model fit when compared to the simple model of relative yield versus O 3 concentration (delta-AIC = 42.1). It is therefore likely that the year of cultivar release, country of study and type of exposure all have some separate influence on the sensitivity of the response of soybean to O 3 . The presence of some collinearity between the candidate explanatory variables, and the observation that many of the AIC values representing different model configurations are very similar, means that we should be cautious when trying to rank the variables in order of influence. The metric conversion variable was not present in the 'best' model, and it is therefore likely that only minimal error was introduced to the data set through O 3 concentration metric conversions.
The effect of country of study on soybean sensitivity to O 3 was investigated further by fitting separate regression lines to the combined dose-response data set according to country. Dose-response data from Indian and Chinese studies were seen to exhibit a steeper decline in yield with increasing O 3 concentration than the data from the USA (Fig. 4a). The response function based on USA data alone predicts a RYL c,p of 16.5%, relative to pre-industrial levels. The Indian and Chinese functions predict a RYL c,p of 30.3% and 33.3%, respectively. The interaction between O 3 concentration and country was highly statistically significant in a separate regression analysis carried out to investigate the individual country effect (P = 0.0015, F = 6.625, df = 348). There was no significant difference between the dose-response functions for India and China (P = 0.79 when the India-O 3 and China-O 3 interactions are compared). Their data was therefore combined to produce a more robust 'Asia' function based on more data points (Fig. 4b).
The individual effect of exposure method on the observed sensitivity of soybean to O 3 was also Fig. 2 Diagram illustrating how % relative yield reduction estimates reported in the results and discussion of this manuscript were calculated. Pre-industrial yield, predicted by the dose-response function, was treated as the 100% baseline yield (relative yield = 1), relative to which yields at present-day O 3 concentrations were expressed. investigated. Data from FACE experiments were seen to exhibit a steeper dose-response relationship than data collected in OTC's (Fig. 5). A linear regression analysis to investigate the individual effect of exposure type found the interaction of exposure type and O 3 concentration to be of borderline statistical significance (P = 0.048, F = 3.93, df = 364.17). Figure 6a distinguishes the data points in the combined soybean dose-response data set by the decade of cultivar release. Modern cultivars, represented on the plot by darker colours, tend to represent the steeper side of the dose-response distribution. A separate linear regression analysis on the 25 soybean cultivars with three or more supporting data points showed that cultivar sensitivity to O 3 has increased over time (Fig. 6b). The regression analysis was carried out twice, once with and once without the Indian cultivars. The sensitivitytime function comprising data exclusively from the USA is the one that avoids the possibility of bias due to geographical differences in sensitivity. This function estimates that the average slope of the soybean doseresponse relationship would have been À0.0040 in 1960 and À0.0053 in 2000, representing an increase in the dose-response slope of 32.5%, over a period of 40 years.

Discussion
The combined dose-response function for soybean in Fig. 3 Fig. 4 (a) Subdivision of soybean dose-response data by the country in which data collection took place and (b) with the data for China and India combined into one dose-response function ('Asia'). Dose-response functions are as follows: USA, y = À0.0047x + 1.020 (df = 323, P < 0.001). India, y = À0.0079x + 1.015 (df = 9, P < 0.001). China, y = À0.0084x + 1.00 (df = 16, P < 0.001). Asia, y = À0.0081x + 1.01 (df = 26, P < 0.001)  against the year in which they were released to market. Two regression lines are shown: one which has been fitted to all cultivars (df = 23, P = 0.0019, r 2 = 0.32), and one which has been fitted to cultivars tested in the USA only (df = 20, P = 0.0271, r 2 = 0.18), excluding the data for Indian cultivars which are circled. Linear equation for all cultivars: y = À0.000058x + 0.11. USA-only linear equation: y = À0.000032x + 0.06. levels as previously published functions. RYL c,p is estimated to be 17.2% using our function, compared to 16.2% and 18.9% predicted by the functions of Mills et al. (2007) and Lesser et al. (1990), respectively. However, the dose-response relationship presented in this manuscript is linear, with 100% relative yield occurring at a theoretical background O 3 M7 value of zero. This is in contrast to the Mills et al. (2007) function which is based on the AOT40 metric and therefore assumes that O 3 concentrations below 40 ppb are not contributing to effects. The dose-response function for soybean published by Lesser et al. (1990) is in Weibull form and is therefore nonlinear, although the curve is very slight and much closer to a linear model when compared to other crop dose-response functions calculated from the NCLAN experiments (Wang & Mauzerall, 2004). Both of the previously published soybean dose-response functions are based only on data from the USA and do not include any data published after 1998. Our doseresponse function is therefore the most comprehensive published to date and predicts that some soybean yield reduction will occur even at low concentrations of ambient O 3 , consistent with the previously published Weibull function for soybean . The critical level for soybeandefined as the O 3 concentration threshold at which statistically significant yield reduction (5%) can be observed (Mills et al., 2007) is predicted using our dose-response function to be 32.3 ppb M7, when calculated relative to pre-industrial O 3 levels (M7 of 23 ppb). This is in line with the 32.4 ppb M7 critical level estimated by the function of Lesser et al. (1990) but a lower estimate than the 40.3 ppb M7 level predicted by the function in Mills et al. (2007), when both are converted to the M7 metric using the conversion functions presented in Fig. 1. The dose-response functions presented in this manuscript for India and China predict slightly lower critical levels of 28.3 ppb and 27.8 ppb M7, respectively.
Further analysis of cultivar sensitivity within the dose-response data set has revealed several important trends. The first is the significant positive correlation observed between soybean cultivar sensitivity and the year of release. Based on the sensitivity-time relationship calculated from the USA cultivars only, O 3induced RYL c,p is estimated to be on average 14.1% for cultivars released in 1960, compared to 19.3% for cultivars released in 2000. This change in cultivar sensitivity is considered to be a conservative estimate. The sensitivity-time relationship which includes the Indian cultivars estimates a greater change in cultivar sensitivity over time, with RYL c,p increasing from 13.1% in 1960 to 22.6% in 2000. However, this steeper sensitivity-time function incorporating the Asian cultivars could be artificially steep if differences in sensitivity due to geographical location are also influencing the values. The trend we have identified in cultivar sensitivity to O 3 over time is in line with the results of a number of studies conducted for wheat, which found modern wheat cultivars to have greater O 3 -sensitivity than older ones (Barnes et al., 1990;Velissariou et al., 1992;Pleijel et al., 2006;Biswas et al., 2008), although this study is to our knowledge the first evidence of this phenomenon in soybean.
The mechanism underlying this temporal trend in sensitivity is unclear, although it may be linked to varietal improvement strategies. Selective breeding across different world regions has transformed the agronomic characteristics of soybean cultivars over the last halfcentury (Morrison et al., 2000;Jin et al., 2010;Agarwal et al., 2013;Koester et al., 2014;Rincker et al., 2014). As well as having dramatically higher seed yield, modern varieties also have higher net photosynthetic rate, chlorophyll content and transpiration rate and have lower leaf area index and shorter maturation periods compared to older varieties Miladinovi c et al., 2015). It is possible that agronomic traits which have been targeted by crop breeders are mechanistically linked to physiological traits associated with O 3 sensitivity, such as a low antioxidative capacity and high g max Biswas et al., 2008). For example, selection for high yield could have simultaneously targeted a high g max to facilitate greater CO 2 fixation (Roche, 2015). This hypothesis is supported by results from a study on 24 soybean cultivars with release dates spanning 1923 to 2007, which observed an increase in g s with year of release in cultivars which also exhibited increasing instantaneous rates of carbon uptake with year of release (Koester et al., 2014). The g s of wheat cultivars has also been reported to progressively increase with their year of release and correlates positively with O 3 sensitivity (Biswas et al., 2008). Breeding for a high harvest index and rapid maturation over recent decades (Morrison et al., 2000;Jin et al., 2010) may have also played a role in the greater O 3 sensitivity of modern cultivars of soybean, by selecting for a trade-off which prioritizes vegetative and reproductive growth over antioxidant synthesis, which could be associated with a metabolic cost under O 3 enriched conditions (Huot et al., 2014;Frei, 2015).
A net increase in the yield of soybean cultivars has taken place over recent decades despite their increasing sensitivity to O 3 . The heterogeneity of O 3 concentrations temporally and geographically may explain the lack of sufficient natural selection pressure for O 3 tolerance at cultivar breeding sites (Ainsworth et al., 2008). Cultivar breeding programmes focussing on enhancing the ability of varieties to detoxify O 3 would increase tolerance and improve yield further (Frei, 2015). Another approach for breeding O 3 tolerance would be to select for reduced g max to reduce the rate of O 3 flux into the plant, and faster stomatal dynamics to allow leaves to close their stomata more rapidly in response to O 3 stress (Morgan et al., 2003). While the reduction in photosynthetic gas exchange associated with excluding O 3 could result in a small yield penalty during less polluted years, cultivars with reduced g max would likely perform better during years with high levels of air pollution, perhaps resulting in an average yield gain over time. A similar strategy in soybean with drought tolerant traits has shown early signs of success, with a 50year simulation based on US weather data showing a significant improvement in average yields, despite some of the traits being detrimental in wet years (Sinclair et al., 2010).
A second important pattern identified in the data analysis relates to the observed geographical variation in O 3 sensitivity. A steeper decline in soybean yield with increasing O 3 was observed in experimental data collected in India and China, compared to data from the USA. Unfortunately, a limited amount of doseresponse data was available for the Asian region: two studies from India and one from China met the inclusion criteria for analysis, with Asian cultivars comprising 12 of the 49 cultivars and 30 of the 379 data points included in the complete data set. Despite the small number of data points representing the Asian region in our analysis, the interaction between O 3 concentration and country of data collection exhibited a high level of statistical significance in the individual regression describing the variation in soybean yield response to O 3 (P = 0.0015), and country of study emerged as an important variable in the stepwise multiple regression.
The greater sensitivity observed in the Asian data suggests that the use of region-specific dose-response functions could potentially improve the accuracy of modelled crop loss estimates. It also highlights the urgent need for more O 3 exposure studies in India and China, which are currently significantly underrepresented in the dose-response literature compared to the USA. Historical and contemporary O 3 trends in India and China are not well documented (Cooper et al., 2014), but both countries have seen a rapid increase in emissions of O 3 precursors as a result of rapid urbanization and industrialization (Granier et al., 2011) and are likely to experience significant increases in surface O 3 concentrations by 2050 (Fiore et al., 2012). Ozone modelling in South Asia by Engardt (2008) based on emissions for the year 2000 estimated surface O 3 concentrations during the soybean growing season (September-November) to be 40-45 ppb M7 over large areas of the state of Maharashtra, which produces over 30% of India's total soybean crop (DAC, 2014). The dose-response function combined across all regions predicts relative yield reduction at this O 3 concentration to be 9.2-11.9% relative to pre-industrial levels, while the India-specific response function estimates yield reduction to be 16.2-20.9%a large discrepancy of estimation. Over large areas of the agriculturally important Indo-Gangetic plain where soybean is also grown (Singh, 2006), modelled surface O 3 exceeds 49 ppb M7, with soybean yield reduction here estimated to be 24.8% using our Indian response function. Accurate estimates of potential O 3 effects on crop yield is arguably particularly important for the South Asian region, where 21% of the population are currently undernourished, an estimated 51% of soybean cropland is reported to be experiencing stagnating or declining yields (Ray et al., 2012), and average soybean yield per hectare is less than half that in the USA (Panthee, 2010).
The mechanism underlying the differential sensitivity observed in North American and Asian doseresponse data is unclear. Interestingly, greater O 3 sensitivity of Asian cultivars compared to North American ones has been previously observed in wheat and rice . Differences between the climate and environment of the different geographical regions could be one factor driving the observed difference in sensitivity. Large areas of China and India experience a humid subtropical climate (Rubel & Kottek, 2010), which facilitates high g s and therefore high O 3 flux. Similarly, warm temperatures correlate with high g s up to a species-specific optimum temperature, above which conductance falls (Emberson et al., 2000). The cooccurrence of O 3 concentration peaks with periods of high humidity and optimum temperaturewhich follow seasonal and diurnal patterns specific to geographical regionscould therefore be a significant factor in determining the degree of crop loss. Unfortunately, the wide range of study locations, open-air experimental designs and seasonal duration of experiments included in this analysis meant that humidity and temperature could not be investigated when synthesizing the data.
Another important factor which must be considered when interpreting the Asia data is the possibility of interactions between O 3 and other ambient air pollutants. There is some evidence that the simultaneous or sequential occurrence of O 3 with SO 2 , NO 2 and NH 3 can have a greater than additive effect on the yield of crops (Fangmeier & Bender, 2002;Bender & Weigel, 2011). Two of the three experimental studies included in our analysis which took place in Asia added O 3 to nonfiltered air, and concentrations of other ambient air pollutants were not recorded during these experiments. The potential for O 3 interactions with other pollutants means we must interpret the higher sensitivity of soybean observed in the Asian studies with some caution. However, all of the data points collected in Asiaincluding those from the experiment which added O 3 to carbon-filtered air (Singh et al., 2010) and therefore removed other ambient pollutantslie below the dose-response line fitted to USA-only data (Fig. S4), suggesting that multipollutant interactions are not the sole driver of the greater sensitivity of Asian dose-response data.
As discussed earlier in relation to temporal trends, crop breeding strategies may be partly driving the observed regional differences in O 3 sensitivity. Crop breeding strategies in the USA, China and India over the last half-century have shared the common aim of increasing yield and harvest index (Karmakar & Bhatnagar, 1996;Morrison et al., 2000;Jin et al., 2010), but other breeding targets are likely to have varied by region. For example, the high sensitivity of soybean to day-length means that maturation periods are highly tailored for different latitudes (Agarwal et al., 2013). In addition, region-specific efforts to breed resistance to local diseases or pests could have increased the capacity of cultivars to upregulate antioxidants, potentially increasing their tolerance to O 3 (Bowler et al., 1992).
The third key result from the data analysis is that the sensitivity of soybean cultivars to O 3 varies widely, and varieties introduced at a similar time and from the same geographical region also exhibit a certain degree of variation in sensitivity. For example, 'Corsoy-79' and 'Hodgson'both released in the USA in the same decade (1970s)are predicted using the functions calculated in this study to experience a RYL c,p of 18.1% and 13.3%, respectively, relative to pre-industrial O 3 . A wide range of within-species variation in O 3 sensitivity has been observed before in other crop species. Quarrie et al. (2007) studied 95 wheat cultivars and observed yield reduction ranging from 0% to 56% following season-long O 3 exposure at an M7 of 91 ppb. Further evidence of differential cultivar sensitivity in wheat has come from studies on Chinese (Biswas et al., 2008) and Bangladeshi (Saitanis et al., 2014) varieties. A similar range of sensitivity has also been observed in Thai rice cultivars (Ariyaphanphitak et al., 2005). The variation in O 3 sensitivity among cultivars observed in this study suggests that there is substantial scope for breeding O 3 -tolerant soybean varieties.
The difference in the yield response observed in FACE and OTC's should be interpreted with caution, due to the marginal P-value in the individual regression (P = 0.048), and the presence of some collinearity. FACE data exhibited a marginally steeper doseresponse slope compared to data collected in OTC's. This result indicates that both methods of exposure produce dose-response data that is comparable and that the impact of the ʻchamber effectʼthe alteration of the growth environment in OTC's which can lead to heightened temperatures, altered air flow and greater vapour pressure deficit (Sanders et al., 1991;Long et al., 2005)on the soybean yield response to O 3 is only small, if it exists. More work is needed to confirm or reject the possibility that exposure method impacts the yield response of crops in O 3 exposure studies.
In conclusion, this study has revealed a large degree of intercultivar variation in soybean O 3 sensitivity and has also identified temporal and geographical patterns in sensitivity. These patterns are relevant to efforts in breeding O 3 -tolerant crop cultivars and also to those carrying out global modelling assessments of O 3 impacts on crop yield. This manuscript has discussed potential factors which might be playing a role in driving these patterns, but they are not yet fully understood. The derivation of flux-based dose-response relationships for soybean, which estimate O 3 exposure based on known relationships between climatic conditions and g s (Emberson et al., 2000), could shed light on the hypothesis that local climatic factors and particular physiological traits related to gas exchange are driving the observed regional and temporal patterns in sensitivity.