Neural Network Analysis to Evaluate Ozone Damage to Vegetation Under Different Climatic Conditions

Tropospheric ozone (O3) is probably the air pollutant most damaging to vegetation. Understanding how plants respond to O3 pollution under different climate conditions is of central importance for predicting the interactions between climate change, ozone impact and vegetation. This work analyses the effect of O3 fluxes on net ecosystem productivity (NEP), measured directly at the ecosystem level with the eddy covariance (EC) technique. The relationship was explored with artificial neural networks (ANNs), which were used to model NEP using environmental and phenological variables as inputs in addition to stomatal O3 uptake in Spring and Summer, when O3 pollution is expected to be highest. A sensitivity analysis allowed us to isolate the effect of O3, visualize the shape of the O3-NEP functional relationship and explore how climatic variables affect NEP response to O3. This approach has been applied to eleven ecosystems covering a range of climatic areas. The analysis highlighted that O3 effects over NEP are highly non-linear and site-specific. A significant but small NEP reduction was found during Spring in a Scottish shrubland (−0.67%), in two Italian forests (up to −1.37%) and during Summer in a Californian orange orchard (−1.25%). Although the overall seasonal effect of O3 on NEP was not found to be negative for the other sites, with episodic O3 detrimental effect still identified. These episodes were correlated with meteorological variables showing that O3 damage depends on weather conditions. By identifying O3 damage under field conditions and the environmental factors influencing to that damage, this work provides an insight into O3 pollution, climate and weather conditions.


INTRODUCTION
Tropospheric ozone (O 3 ) is a harmful air pollutant which affects human health (Ainsworth et al., 2012), damages vegetation, including natural ecosystems and crops (The Royal Society, 2008), and contributes to climate change, being a greenhouse gas with a radiative forcing of 0.35-0.37 W m −2 (Shindell et al., 2009). It is a secondary pollutant, mainly produced through photochemical reactions of methane, carbon monoxide and volatile organic compounds in the presence of nitrogen oxides (Monks et al., 2015).
Although O 3 is a natural component of the troposphere, its concentration has been increasing since the pre-industrial era as a result of anthropogenic emission of its precursors (Ainsworth et al., 2012). Nowadays, the background O 3 mixing ratio of the northern hemisphere is 30 to 40 ppb (Parrish et al., 2012), although large regional differences are recorded due to the strong influence of weather, which promotes O 3 formation in warm, dry and sunny conditions (see Table 2 for mean O 3 mixing ratio at our study sites).
Following chemical destruction, the second most important sink of tropospheric O 3 is the dry deposition on land surfaces, primarily controlled by vegetation, which contributes to O 3 removal through stomatal uptake, deposition onto the surface, and in-canopy chemistry (Stevenson et al., 2006). Major O 3 uptake occurs at leaf level, controlled by stomatal absorption (Cieslik, 2004;Fowler et al., 2009). Entering the leaves through stomata, O 3 sets off a chain of oxidative reactions within the apoplast, damaging cell metabolism (Wohlgemuth et al., 2002). The main detrimental effect is a reduction in carbon assimilation, which represents the first evidence of O 3 impact over vegetation, before the occurrence of visible injuries.
Manipulation experiments have been widely used to assess the impact of O 3 over vegetation carbon assimilation capacity (Karlsson et al., 2000;Manning, 2005;Gerosa et al., 2015). While this approach has been useful in understanding vegetation behavior in standard conditions and to derive dose-response functions, it may be inadequate to provide the complete picture, since plants are often exposed to unrealistic concentrations, the approach is often limited to young plants, and the experimental facilities alter the microclimate.
An ecosystem approach is thus of primary importance for understanding how O 3 pollution affects CO 2 uptake by vegetation and to assess the validity of extrapolating the effect (Sitch et al., 2007). Eddy covariance (EC) towers, from which the carbon flux is measured with a wide range of meteorological variables at high temporal resolution, provide large datasets which can be used to extrapolate information about ecosystem responses to O 3 pollution (Fares et al., 2018).
A useful tool for investigating functional relationships between site characteristics and environmental factors such as climate and other atmospheric conditions is Artificial Neural Network (ANN) modeling (Aitkenhead and Coull, 2016). ANNs are very powerful in analyzing and modeling non-linear relationships owing to their capacity to learn from examples and generalize, allowing them to explore relationships without making assumptions about the shape of these relations (such as are made by other approaches such as multiple linear regression) (Olden and Jackson, 2002;Moffat et al., 2010). Although ANNs are primarily used in the building of predictive models, methods for quantifying the independent variable contributions within networks have also been developed (Olden et al., 2004), allowing researchers to use them to understand how climate variables drive ecosystem responses (Moffat et al., 2010).
In this work, feed-forward ANNs were used to test the hypothesis that current O 3 concentration affect vegetation photosynthetic CO 2 assimilation under field condition by isolating the effect of O 3 on the net ecosystem productivity (NEP) of eleven ecosystems, characterized by different climatic condition and O 3 concentration, taking into account the influence of other NEP climatic drivers (solar radiation, air temperature, vapor pressure deficit, soil water content) and stomatal conductance. This approach is fully empirical and avoids a priori assumption on the functional relationships between the study variables, which are measured directly. The analysis was conducted using daytime eddy covariance (EC) data directly measured over: eight northern hemisphere open tree canopies; one moorland; one grassland and one cropland. We had the following aims: (1)

Study Sites
Eleven sites from three eddy covariance flux measurement networks were selected to test the effect of O 3 pollution on NEP different type of vegetation: four semi-natural forests, three planted forests, one orange orchard, one moorland, one grassland and one cropland. The study sites are spread across five different Köppen climatic areas: Oceanic (Auchencorth Moss, Grignon, Lochristi, Speulderbos), Semi-arid-Continental (Bugac), Subartic (Hyytiälä), Humid-Subtropical (Bosco Fontana, Ispra) and Mediterranean (Castelporziano, Blodgett, Lindcove). Information about site location and ecosystem type can be found in Table 1 and a map showing the ECLAIRE site locations can be found in Fares et al. (2018).

Datasets
Data used in the development of the ANNs was recorded continuously from EC towers, at half-hour intervals, from January 2012 to December 2013 for all sites, except for: (1) Blodgett where data was collected from January 2001 to December 2007; (2) Lindcove which operated from 2009 to 2010; (3) Castelporziano for which data were collected from 2013 to 2015 and (4) Grignon, where only the dataset with rapeseed was used (31 August 2012 to 15 August 2013). Only relevant data for detecting O 3 effect over NEP were considered in the study. Last column indicates references where a detailed description of the sites is available.
Since damage occurs due to O 3 stomatal absorption (Reich and Amundson, 1985;Biswas et al., 2007;Broschè et al., 2010), we reduced the dataset to daytime data (10:00 -18:00 UTC time) from the Spring and Summer seasons, when stomata are open and O 3 levels in the atmosphere are high. A summary of data used in this study is given in Table 2.
Although data coverage was incomplete over the measured period at each site, interpolation of large gaps was avoided in order to make sure that the functional relationships captured by ANNs were unbiased. Small gaps (less of 50% of missing data over 10:00-18:00 period in a day) were replaced by the mean of correspondent half-hour data of adjacent days (Falge et al., 2001).
Stomatal conductance to H 2 O (G st , m s −1 ) was calculated as the inverse of stomatal resistance (R s ), derived from measured latent heat flux (E, kg m −2 s −1 ) using the evaporative/resistance method (Monteith, 1981): where cp is the specific heat capacity of air (J kg −1 K −1 ), ρ is the density of the dry air (kg m −3 ) q a is the vapor pressure at measurement height (Pa), q s is the saturation mass fraction (Pa) of H 2 O at air temperature and roughness length z 0 , γ is the psychrometric constant (67 Pa K −1 ) and λ is the vaporization heat for H 2 O (2.5 × 10 6 J kg −1 ). The use of E to calculate R s is valid only if transpiration is the only significant source of water vapor from the ecosystem and thus only data recorded during dry-daylight conditions were used. Data were discarded if they met any of the following criteria: net solar radiation <20 W m −2 , relative humidity > 80%, rainy days (daily rainfall > 2 mm day −1 ) or the day after a rain event. Discarded data are mainly located when R s is also large, hence little influence on overall dose is expected. Percentages of data discarded by this procedure are as follow: 25% for Auchencorth Moss (Au), 11% for Bugac (Bu), 24% for Grignon (Gr), 18% for Hyytiälä (Hy), 19% for Speulderbos (Sp), 31% for Lochristi (Lo), 15% for Bosco Fontana (BF), 18% for Ispra (Is), 15% for Castelporziano (CPZ), 6% for Blodgett (BL) and 11% for Lindcove (Ci).
O 3 stomatal uptake (FO 3sto , nmol m −2 s −1 ) was calculated as: where 0.61 is the ratio of diffusivity between O 3 and H 2 O (Marrero and Mason, 1972) and [O 3 ] canopy is the O 3 concentration at the canopy level. The latter was calculated following the standard resistance analogy (Hicks et al., 1987). A detailed explanation of the calculation can be found in Supplementary Appendix 1.

Artificial Neural Network Modeling
Artificial neural networks (ANNs) were used to model the NEP. An ANN can be defined as a large series of simultaneous equations with each variable equivalent to a simple processing element (node) connected to each other by connection weights. Appropriate values within the connection weights provide the network with the ability to store knowledge about some modeled system. A supervised learning algorithm (i.e. with predefine input and output values within the training data) adjusts the connection weights, randomly assigned at the beginning, to approximate relationships that are present in the data. Three ANN model runs were conducted, trained using three different groups of input variables. The three cases were used to test the hypothesis that O 3 influences NEP: Case 1. ANNs were trained using solar radiation (St, W m −2 ), air temperature (T, • C), vapor pressure deficit (vpd, kPa), soil water content measured between 10 and 30 cm depth (swc,%, not available for BL) and stomatal conductance to H 2 O (G st , m s −1 ). Case 2. Included all Case 1 input variables plus O 3 stomatal uptake (FO 3sto , nmol m −2 s −1 ). The latter was included in the analysis under the assumption that, if O 3 has a detrimental effect on vegetation (and thus on NEP), it would be caused by the O 3 entering the leaves. If O 3 absorbed through stomata affects NEP, the ANN model's ability to predict NEP would be improved. Case 3. consisted of Case 1 input variables plus the O 3 dose absorbed through stomata integrated over 3 h before the measure time (FO 3cum , µmol s m −2 ). This case was used to test if accumulated O 3 entering the stomata was a better predictor of NEP damage than instantaneous O 3 stomatal uptake, under the assumption that antioxidants are consumed by O 3 during the day and a longer exposition to high O 3 level reduces the leaf capacity to detoxify O 3 entering the stomata due to fast scavenging of antioxidant defense in the intercellular spaces. The 3-h interval was chosen because it was the longest interval which permitted us to not include the night-time data (i.e. at 10 am, FO 3cum integrated measurements between 7 am and 9 am).
In Case 1, O 3 variables were not considered so that, if the model performance was better using Case 1 rather than Case 2 or 3, O 3 had no effect on NEP.
All input variables and NEP values within the dataset were normalized by scaling between 0 and 1, to ensure that no variables had an inherently greater effect than others. The dataset was split into two subsets, Spring (from 21st of March to 20th of June) and Summer (from 21st of June to 22nd of September), and ANNs were trained separately for each subset, with the aim of highlighting the seasonal variability of the response of NEP to O 3 . Feed forward ANNs with a sigmoid activation function were used.
The feed-forward ANN was made of three layers: one input layer, a single hidden layer and an output layer. In a feed-forward ANN the information flows only in a forward direction, from the input to the output through the hidden layer. Layers are composed of nodes, with nodes in adjacent layers fully interconnected by weights which are determined by a supervised learning algorithm appropriate for non-linear regression (backpropagation algorithm, Rumelhart, 1986). In this work, the hidden layer consisted of 8 nodes. The number of nodes of the hidden layer war chosen by comparing the performance of different networks, with 1 to 10 hidden nodes, and choosing the number that produced the best network performance (Gevrey et al., 2003;Olden et al., 2004).
For each group, ANN training was repeated 100 times, because different ANNs trained with the same dataset may return different connection weights, depending on the training procedure and initially randomized connection weights. A common criticism of ANN modeling is "overfitting, " which is the case that ANN memorizes the training data but may fail to fit new data (Chan et al., 2006). Overfitting occurs when the model is parameterized to give the best possible fit to the training data, rather than to the "global dataset" possible from all possible examples of the system being studied. While this is a risk of all data mining or statistical regression approaches, the same solution can be applied as here: datasets were split randomly into three subsets: training (70% of dataset), test (15% of dataset) and validation (15% of dataset). The training subset was used to compute the weights of the network's nodes and the test subset for stopping the training process and checking the model generalization ability. The validation subset was used to validate the model and prove the ANN's ability to generalize beyond the training dataset.
Artificial neural network development and training was carried out using Neural Network Toolbox (Matlab 2010, Natick, MA, United States).

Performance of ANNs
One-way analysis of variance (ANOVA, confidence interval 95%) was used to determine whether there were any statistically significant differences between the means of the original measured NEP and the 100 modeled NEP values derived from ANN simulations of each case. If a statistical difference was found, a post hoc test was performed to detect which specific simulation differed from measured NEP, in order to discard that simulation and train the ANN again. The coefficient of determination (r 2 ) was used as a measure of goodness of fit, and as an indicator to evaluate if the inclusion of O 3 parameters into ANN models improved the model ability to simulate NEP, thus suggesting an effect of O 3 over NEP.
Artificial neural network model was compared with a linear statistical approach, Multiple Linear Regression (MLR). MLR model is used to explore the relationship between a dependent variable and independent variables, under the assumption that each independent variable has a linear relationship with the dependent variable (Civelekoglu et al., 2008). In this work, MLR was used to model the linear relationship between NEP and the three groups of input variables (case 1, case 2 and case 3, see above) which were the same input variables of ANN modeling. The MLR r 2 was calculated and compared with ANN r 2 , in order to evaluate if the ANN approach better performed than the linear approach in predicting NEP behavior.

Analysis Tools for Quantifying O 3 Contributions as NEP Driver in ANN Modeling
The integrated information gathered from ANNs can be decomposed to disentangle the effects of different inputs on the output values, to improve understanding of how each input variable affects the predictions. Gevrey et al. (2003) and Olden et al. (2004) provided a comparison of the different existing methods for estimating variables importance in ANN applications. In this work, the partial derivative method (Dimopoulos et al., 1995) was used to isolate the effect of O 3 over NEP estimated by ANN modeling.
The partial derivative method produces a profile of the output variations for unit change of selected input variable. The link between the modification of the input, x j , and the variation of the output, y j = f(x j ), is the partial derivative of each activation function with respect to its input (d j ), with j = 1,. . .,N and N the total number of observations.
Given an ANN with n inputs i (i = 1,. . .,n), one hidden layer with m h nodes h (h = 1,. . .,m h ) where the logistic sigmoid function is used for activation, the partial derivative of y j with respect to x i is d ji (Dimopoulos et al., 1999;Gevrey et al., 2003): where S j is the derivative of the output with respect to its input, I hj is the response of the h hidden node, w ho is the weight between the output node and h, w ih is the weight between h and the input node (n i ). Partial derivatives were calculated for each of the ANN runs and averaged to calculate the mean absolute change of NEP associated with O 3 . To get information about positive and negative change of NEP, positive and negative fraction of the partial derivatives were averaged separately (Moffat et al., 2010).
The weather influence on O 3 down-regulating effect was tested using Spearman partial correlation. This is a nonparametric measure of rank correlation that assesses monotonic relationships of two variables whilst controlling for other, potentially confounding variables. The negative fraction of partial derivative associated with O 3 (FO 3st or FO 3cum ) was correlated with environmental factors such as solar radiation, air temperature, vapor pressure deficit, soil water content and O 3 concentration at canopy level. The latter was included in the analysis with the aim of controlling the confounding effect it may have on correlation coefficients, since O 3 concentration strongly depends on weather (Monks et al., 2015). All environmental factors were transformed between 0 and 1 to avoid scale effects.

Performance of ANNs
The ANOVA test highlighted that there are no statistically significant differences between the means of the original measured NEP and the 100 modeled NEP derived from ANN

Do Current O 3 Levels Affect NEP?
The inclusion of O 3 stomatal uptake in the ANN simulation did not change ANN performance (i.e. no impact of O 3 ) in some cases, and in others it improved the ANN performance (i.e. O 3 had an impact). In particular, the results suggest that O 3 damage does not occur in most of the northern sites which are less exposed to O 3 pollution, while in the other ecosystems a limited effect was observed.
Artificial neural networks trained with Spring data ( Since the 3-h O 3 dose was calculated from cumulated instantaneous O 3 stomatal fluxes, r 2 between the two parameters was calculated at each site to verify the degree of independence of the two variables. They were shown to be fully independent (r 2 = 0.00) for Lindcove, Blodgett, Castelporziano, Bosco Fontana, Speulderbos and Hyytiälä. An r 2 value of 0.00 was also found for Auchencorth Moss and Ispra in Summer, while in Spring r 2 were 0.10 and 0.25, respectively. Higher r 2 were found for Grignon (0.45 and 0.73 for Spring and Summer, respectively), Bugac (0.37 and 0.42 for Spring and Summer, respectively), and Lochristi (0.25 and 0.48 for Spring and Summer, respectively).

Assessing the Sensitivity of NEP to Current O 3 Levels
Partial derivatives represent NEP rate of change with respect to O 3 (FO 3st or FO 3cum ). If this rate is negative, then the NEP will tend to decrease as FO 3st or FO 3cum increases while if the rate is positive, NEP will tend to increase. Partial derivatives were calculated for each half hour observation. Partial derivatives calculated with respect to one predictor can be positive for some half hours and negative for other half hours. This means that a predictor has not always a positive or negative effect on NEP, and that it depends on the combination of all ANN predictors values occurring at that time.
Averaging separately positive and negative partial derivatives helps to discern when O 3 has a negative effect over NEP. Figure 1 shows that, although the O 3 variables were significant factors for predicting NEP in the ANN model runs, they did not always lead to a reduction in NEP. Reduction of NEP related to O 3 was detected at Auchencorth Moss, Bugac, Grignon, Bosco Fontana, Castelporziano, Blodgett and Lindcove in Spring, Grignon, Bosco Fontana and Lindcove in Summer (Figure 1), with values ranging from 0.15 to 2.64% average NEP loss due to O 3 ( Table 5).
In a few cases the response is only negative, indicating that current O 3 concentration level consistently reduces NEP during Spring or Summer, for all combinations of the other predictors. This is the case of Bosco Fontana in Spring and Lindcove in Summer. In other cases, the ANN did not always determine O 3 to be a damaging factor for NEP. To highlight how the rate of change of NEP responds to different levels of O 3 entering the stomata, a profile of the NEP partial derivative versus O 3 input (FO 3st or FO 3cum depending on the case) was plotted (Figures 2, 3).
The relationships shown in Figures 2, 3 are highly non-linear and present different behaviors at each site: at Auchencorth Moss site (Figure 2A) a down-regulating effect was found in Spring for FO 3sto below the 40th percentile (−0.27 nmol m 2 s −1 ), peaking around 23rd percentile (−0.20 nmol m 2 s −1 ). The same trend was observed at the Castelporziano site ( Figure 2C) during Spring, where FO 3sto below the 33rd percentile (−0.13 nmol m 2 s −1 ) negatively affected NEP. At the Bugac site (Figure 2B), the FO 3sto effect on NEP was almost linear, and the effect turned negative above the 52nd percentile (−0.19 nmol m 2 s −1 ). At the Blodgett site, during Spring (Figure 2E), FO 3cum negatively affected NEP in a range between 29th (68.33 µmol m −2 ) and 49th (91.23 µmol m −2 percentile, peaking at the 40th (80.67 µmol m −2 ) percentile.
At the Bosco Fontana site the ANN analysis predicted a consistently negative effect of FO 3cum flux for the Spring data (Figure 2D), although the largest effect was recorded for the smallest doses (below 20th percentile, 43.68 µmol m −2 ). During Summer (Figure 3C), a depressing effect of FO 3cum over NEP was observed only below the 26th percentile (70.72 µmol m −2 ). At the Lindcove site, during Spring (Figure 2F), only low doses of FO 3cum had a negative effect on NEP (below 7th, 45.42 µmol m −2 , and above 95th percentile, 151.02 µmol m −2 ), whilst during Summer (Figure 3A), all FO 3st values induced a decrease in NEP, peaking at 88th percentile (0.19 nmol m −2 s −1 ). At the Grignon site, the  (Figures 4, 5). During the Spring season, the average dNEP(O 3 ) followed a pronounced bell-shape curve at the Lindcove site, where the maximum effect of O 3 was observed during the middle hours of the day, coinciding with FO 3cum peak ( Figure 4F). The same pattern was found for Bugac ( Figure 4B) sites, although the bell-shape was slightly accentuated. At Auchencorth ( Figure 4A) and Castelporziano ( Figure 4C) larger effects of FO 3st were recorded at the end of the afternoon. The damaging effect of instantaneous O 3 followed an exponential decrease during the day at Bosco Fontana site (Figure 4D), opposite to the shape of the response to FO 3cum . The damaging effect of O 3 peaked at 10:00 h, when FO 3cum was still low. At Blodgett (Figure 4E) no significant variations were found during the day.
For the Summer season, average dNEP(O 3 ) showed a bellshape curve which follows the FO 3st trend at Grignon and Lindcove (Figures 5A,B, respectively), while at Bosco Fontana ( Figure 5C) the dNEP (FO 3cum ) presents the same pattern as during the Spring season.

How O 3 Reduces NEP According to Other Environmental Factors
For the sites where a negative effect of O 3 over NEP was found (Au, Bu, Gr, BF, CPZ, Bl and Ci), the correlation between the negative fraction of dNEP(FO 3st ) or dNEP(FO 3cum ) and the other environmental variables was tested through Spearman partial correlation (

DISCUSSION
How O 3 pollution alters vegetation carbon sequestration capacity is considered an important component of global change (Ashmore and Bell, 1991), but few studies have quantified its impact over ecosystems. Some of them confirm a detrimental effect of O 3 over vegetation occurring in sites where ambient O 3 concentrations are typically high (Zapletal et al., 2011;Fares et al., 2013), while others did not find any effect of high levels of tropospheric ozone concentrations (Zona et al., 2014;Verryckt et al., 2017). Results of these studies are difficult to interpret in the context of ozone/plant interactions because of the great variability among site characteristics, vegetation type and methodological approaches (Cailleret et al., 2018).
This work demonstrates that ANN modeling is a useful tool to understand O 3 -NEP correlation considering other co-varying environmental factors. r 2 values produced by ANN were found higher than r 2 values produced by MLR, indicating that a non-linear statistical data modeling approach as ANN is more appropriate in modeling complex relationships such the dependence of NEP from co-varying environmental factors. Our results are in line with other ecological studies in literature which compared the two methods and found ANN models more accurate than MLR (e.g. Lek et al., 1996;Paruelo and Tomasel, 1997;Brion et al., 2005). The strength of ANN lies in its fully inductive approach, which allows multidimensional relationships to be investigated without a priori knowledge of the shape of these relations. Coupling ANN analysis with EC provided a picture of the current status of O 3 pollution effects over ecosystems.
Although overfitting was controlled in this work, high r 2 values calculated between simulations and measured data ( Table 3) indicates that the possibility of overfitting from ANN exists.
We did not carry out an assessment of whether the model accuracy was significantly better for the training or test subsets than for the final validation subset in each case, and so cannot provide an indication of whether or not the model was actually overfitted to the training data. However, in this work, the power and flexibility of ANN in fitting the data is an advantage and not a limitation, indeed the aim was not to find a general model for NEP, but to evaluate the influence of a variable (O 3 ) in each single study site.  Artificial neural network performance analysis found that both FO 3st and FO 3cum are suitable indicators for predicting NEP reduction, depending on season and type of vegetation. Since FO 3cum is obtained from cumulated values of FO 3st , an analysis of correlation between the two variables was performed with the aims of verifying the degree of dependence of the two variables. None or low correlation was found for most of the study sites. That may suggest that the O 3 dose entering the stomata is discontinuous and that vegetation is subjected to O 3 pulses rather than a constant flux. However, we did not find any relation between the degree of correlation between the two variables and the selection of ANN model best predictor, which means that the most suitable predictor was selected for each site and season regardless of how discontinuous the O 3 fluxes were, and it probably depends on vegetation type and climate conditions.
Of the eleven sites tested, four ecosystems were free from O 3 damage. These sites were the Finnish Scots pine forest (Hyytiälä), the Dutch Douglas fir plantation (Speulderbos), the Belgian poplar plantation (Lochristi) and the Italian  mixed forest (Ispra). These results are in line with Zona et al. (2014), who did not find a negative relationship between O 3 and net ecosystem exchange at Lochristi. All of these ecosystems are located in northern areas with the exception of Ispra, where O 3 concentrations are not particularly high compared with typical high ozone-prone Mediterranean sites ( Table 2) and this may support the hypothesis that low to moderate ozone concentrations and therefore lower stomatal ozone fluxes may generate cumulative exposure to ozone far below possible critical levels in northern ecosystems. Seven ecosystems showed a significant but limited NEP loss due to O 3 entering the stomata: the United Kingdom shrubland (Auchencorth Moss), the Hungarian grassland (Bugac), the Italian Holm oak forest (Castelporziano), the Californian pine plantation (Blodgett) during the Spring season, the French cropland (Grignon) during Summer, the Italian mixed forest (Bosco Fontana) and the Californian Citrus orchard (Lindcove) during both Spring and Summer seasons. Mean NEP loss was estimated at between 0.15 and 2.64%.
These values are low compared with other studies: for example Fares et al. (2013) adopted more traditional statistical methods based on step-wise regression analysis and multivariate analysis and found up to 12-19% of the carbon assimilation reduction in Blodgett and Lindcove sites explained by O 3 entering the stomata. Such results suggest that either our approach is extremely conservative and does not appropriately attribute O 3 effect, or statistical methods adopted in earlier studies may have overestimated O 3 effects by including the effects of covariates in the predictive model of NEP. It should be noted, however, that both approaches focus on the quantification of the instantaneous or nearinstantaneous effect of O 3 on NEP, and capture neither the effect of this NEP reduction on biomass reduction which may further reduce NEP in the future nor the long-term effects of leaf injury. Among the northern sites affected by O 3 , Auchencorth Moss showed higher sensitivity to O 3 damage. This was despite nonstomatal deposition being the principal sink of O 3 at this site, representing 70% of the overall O 3 flux (Fowler et al., 2001), and where moorland mosses like Sphagnum are relatively tolerant to elevated O 3 concentrations (Rinnan et al., 2003). The reason of the depression of NEP linked to O 3 entering the stomata can be attributed to the increase of the plant respiration rate, as already observed by Niemi et al. (2002) for moorland vegetation, as a result of the plants repairing O 3 damaged tissues (Williamson et al., 2015).
Partial derivative results indicate that the O 3 effects on NEP are highly non-linear and site-specific. In almost all study sites, positive relationships between stomatal O 3 flux and NEP were found. Both are controlled by stomatal conductance and thus, in the absence of O 3 damage, stomatal O 3 fluxes positively correlated with NEP, as already observed by Proietti et al. (2016). This also implies that a limited effect of ozone on stomatal closure may still take place, as this matches with a moderate reduction in NEP.
Episodes in which O 3 detrimental effect occurred were identified from the models using partial derivative analysis. These episodes were correlated to climatic variables showing that O 3 damage dependence on weather varies with the climate. O 3 damage occurred primarily during Spring, especially for those sites where stomatal conductance decreases in Summer as affected by water availability (see Table 2). Partial correlation analysis showed that swc decreases the O 3 negative effect for those sites where no drought stress occurred (Blodgett and Bosco Fontana).
In Mediterranean regions, drought periods (which coincide with high O 3 levels) limit stomatal conductance, protecting vegetation from O 3 oxidative stress (Paoletti, 2006). At Lindcove, a well irrigated Mediterranean citrus plantation as previously reported by Fares et al. (2012), O 3 detrimental effects were observed both during Spring and Summer and positively correlate with swc. In these warm periods, with O 3 concentrations often exceeding 80 ppb, the mean stomatal O 3 fluxes were 3.29 and 3.19 nmol m −2 s −1 during the central hours of the day during Spring and Summer, respectively. swc was not a significant predictor of O 3 damage at Castelporziano, a Mediterranean Holm oak forest, where the high water table (Bucci, 2006) protected trees from water stress, although during Spring high stomatal fluxes were associated with high levels of precipitations .
The inclusion of O 3 concentration at canopy level in the partial correlation analysis showed that NEP damage does not always occur at peak O 3 concentrations, as the case of the Italian sites of Bosco Fontana and Castelporziano. It must be noted that being photochemically produced, O 3 concentrations tend to peak when solar radiation is high and that BVOC emitted by some tree species can contribute to ozone formation (Monks et al., 2015). When radiation increases, vpd also increases, causing stomatal closure and leading to a protective effect against O 3 entering the leaves (Mereu et al., 2009;Fares et al., 2010bFares et al., , 2014. For the same sites, partial correlation analysis highlights that when solar radiation and vpd increase, O 3 impact on NEP decreases. We do not exclude that to some minor extend high reactive terpenoid emissions at Castelporziano and Bosco Fontana in the central hours of the day [documented by Fares et al. (2013) and Acton et al. (2016)] may be responsible for ozone scavenging in the gas phase, thus reducing the amount of ozone entering stomata. This phenomenon of mid-day exclusion of O 3 damage does not happen in the Mediterranean Ponderosa pine plantation at Blodgett, this species being relatively insensitive to vpd when drought is not a limiting factor (Panek and Goldstein, 2001). At this site, air temperature was found to be a limiting factor for O 3 damage during Spring, when high temperature constrained gas exchange (Panek and Goldstein, 2001) and helped reduce O 3 oxidative stress.
At Auchencorth Moss, the relation of O 3 damage to weather remains unclear, although a significant negative effect of O 3 on NEP was found. O 3 damage estimates were averaged over the course of the day, with the aim of highlighting hourly patterns. For some sites (Lindcove, Bugac and Grignon), the detrimental effect of O 3 followed a bell-shaped curve, thus suggesting that most of damage to photosynthetic apparatus occurred rapidly during hours of maximum O 3 absorption. Interestingly, for Auchencorth Moss and Castelporziano, major damage was observed at the end of the afternoon, indicating that: (1) high vpd reduces stomatal conductance, and therefore O 3 damage, during the central hours of the day; and (2) to some extent plants may be able to detoxify O 3 during hours of maximum exposure to the pollutant, while at the end of the day detoxification capacity of leaves decreases. While there is evidence of a mid-day depression of stomatal conductance in those sites (Fares et al., 2010a, the second hypothesis is highly speculative and deserves further investigation. However, the possibility of changes in reducing power during the day has been previously described by Dizengremel et al. (2008) who showed decreasing foliar level of antioxidants during the afternoon hours in response to oxidative stress. Conversely, O 3 damaging effect and O 3 absorption were completely decoupled at Bosco Fontana, where the O 3 damaging effect peaks in the morning, under low O 3 concentrations, suggesting the occurrence of species-specific acclimatization phenomena along the day which we cannot explain in this study.

CONCLUSION
This work clearly suggests that long-term datasets are required to identify O 3 damage to vegetation under field conditions. We found that O 3 has a detrimental effect on NEP, although damage can be sporadic and is driven by specific weather conditions and in general, has lower magnitude compared with observations carried out through manipulative experiments or in the field using traditional statistical methods. Our results suggest that vegetation response to O 3 depends not only on pollution level but also on how the ecosystems respond to climate variables. Future climate changes may therefore either expose ecosystems to further O 3 damage by increasing temperatures or rather lead to a reduction in ozone damage in drought-prone ecosystems.

DATA AVAILABILITY STATEMENT
The datasets generated for this study are available on request to the corresponding author.

AUTHOR CONTRIBUTIONS
All authors contributed providing experimental data and/or supporting the preparation of the manuscript.

FUNDING
We want to thank the ÉCLAIRE project (Effects of Climate Change on Air Pollution and Response Strategies for European Ecosystems) funded by the EU's Seventh Framework Program for Research and Technological Development (FP7), which supported the flux measurements at the sites Au, Bu, Gr, Hy, Sp, Lo, BF, and Is. The CREA team wants to acknowledge the LIFE financial instruments of the European Union (LIFE15 ENV/IT/000183 -LIFE18 PRE IT 003) in the framework of the projects MOTTLES ("Monitoring ozone injury for setting new critical levels") and VEG-GAP ("Vegetation for Urban Green Air Quality Plans") LIFE18 PRE IT 003 and the General Secretariat of the Presidency of Italian Republic for financing the CASTEL4 project and the Directorate of Castelporziano Estate. The CEH measurement sites were supported by National Capability funding from the United Kingdom Natural Environment Research Council. The ECOSYS team acknowledge INGOS (grant agreement 284274), the French ANR project ANAEE, and ICOS France.