The Inﬂuence of Ozone on Net Ecosystem Production of a Ryegrass–Clover Mixture under Field Conditions

: In order to understand the effect of phytotoxic tropospheric ozone (O 3 ) on terrestrial vegetation, we quantiﬁed the impact of current O 3 concentration ([O 3 ]) on net ecosystem production (NEP) when compared to the conditions of the pre-industrial era. We compared and tested linear mixed-effects models based on [O 3 ] and stomatal O 3 ﬂux ( F sto ). The managed ryegrass–clover ( Lolium perenne and Trifolium pratense ) mixture was grown on arable land in the Czech Republic, Central Europe. Values of [O 3 ] and F sto were measured and calculated based on resistance analogy, respectively, while NEP was calculated from eddy covariance CO 2 ﬂuxes. We found the F sto -based model more precise when compared to measured NEP. High F sto was found even at low [O 3 ], while broad summer maximum of [O 3 ] was not necessarily followed by signiﬁcant NEP decline, due to low soil water content leading to a low stomatal conductivity and F sto . Comparing to low pre-industrial O 3 conditions, current levels of O 3 resulted in the reduction of cumulative NEP over the entire growing season, up to 29.7 and 13.5% when the [O 3 ]-based and F sto -based model was applied, respectively. During the growing season, an O 3 -induced reduction of NEP ranged between 13.1% in May and 26.2% in July when compared to pre-industrial F sto levels. Looking to the future, high [O 3 ] and F sto may lead to the reduction of current NEP by approximately 13.3% on average during the growing season, but may increase by up to 61–86.6% in autumn, indicating further O 3 -induced acceleration of the senescence. These ﬁndings indicate the importance of F sto and its inclusion into the models estimating O 3 effects on terrestrial vegetation. The interaction between environmental factors and stomatal conductance is therefore discussed in detail.


Introduction
Tropospheric ozone (O 3 ) is formed in UV-driven photochemical reactions from NO 2 , which are further catalyzed by precursor gases, including carbon monoxide, methane, and volatile organic compounds [1,2]. In Europe, ozone concentration ([O 3 ]) has increased since the industrial era from 10-15 ppbv to 30-40 ppbv [3]; it is expected to rise by 4.3 ppbv by et al. [2] for forest tree species. To properly assess O 3 -risk for terrestrial ecosystems, it is thus necessary to consider not only [O 3 ], but also all of the factors influencing stomatal O 3 uptake. For example, reduction in stomatal conductance to O 3 under drought conditions may lead to less plant injury [35]. Moreover, the interactive effects of O 3 with other environmental factors remain unclear [36]. While drought and elevated temperatures reduce the negative O 3 effects on plants, enhanced nitrogen depositions, particularly those applied on managed ecosystems, tend to stimulate negative O 3 impacts [2].
In this study, we evaluate the effect of natural [O 3 ] on net ecosystem production (NEP) of a managed ryegrass-clover mixture grown on arable land. The diurnal courses of NEP and total O 3 flux were continuously measured by the eddy covariance technique over the whole vegetation season from spring (April) to autumn (October) 2019 and coupled with Emberson's resistance-based model [37] to calculate stomatal O 3 flux. Consequently, measured NEP values were compared to those predicted by the linear mixed-effects model parametrized with preindustrial and current values of [O 3 ] and stomatal O 3 fluxes.

Experimental Site Description
The annual ryegrass-clover mixture (total area of 75 ha of which eddy covariance footprint area is 20 ha) is situated on arable land in close proximity to the atmospheric station, Křešín u Pacova in the Czech Republic (N 49 • 35 , E 15 • 05 , 534 m a.s.l.). The atmospheric station is equipped for the measurement of greenhouse gases, selected air pollutants, and basic meteorological parameters [38] and forms a part of European infrastructure networks ICOS [39], GMOS [40], and ACTRIS [41].
The vegetation mixture was established in spring 2018 by sowing spring barley (Hordeum vulgare), garden pea (Pisum sativum), red clover (Trifolium pratense), and perennial ryegrass (Lolium perenne). The seed mixture consisted of 75 kg spring barley, 150 kg pea, 1 kg ryegrass and 18 kg clover sown per hectare. Fertilization (200 kg per hectare) with a mixture of ammonium nitrate and limestone (27% nitrogen) was applied five days after sowing. The legume-cereal mixture was harvested in July 2018 and the clover was harvested in October 2018. In 2019, the year of study, the growth of the original understory with dominating T. pratense and L. perenne continued.
The ryegrass-clover mixture was cut three times during the investigated vegetation season (10 June, 23 July and 21 October 2019). The total harvested dry matter biomass was 8.27, 3.30, and 1.82 tons per hectare, respectively. Leaf area index (LAI) reached maximum values of 8.5, 7.0 and 8.0 m 2 m −2 before the first, second and third cut, respectively, and rapidly decreased at the beginning of September in relation with leaf senescence ( Figure 1C). The LAI was estimated from a time series of the normalized difference vegetation index (NDVI), applying the equation validated for erectophile-planophile canopies [42]. NDVI was computed from red and near-infrared bands using the Google Earth Engine built-in function using Sentinel-2 MSI: MultiSpectral Instrument, Level-2A (surface reflectance) product with a spatial resolution of 10 m × 10 m [43]. The area is characterized by a long-term (2013-2019) annual mean temperature of 8.2 • C and 620 mm of precipitation. Generally, the highest precipitation occurs in July 2019, with an average of 71 mm. February is the driest month, with 23 mm of precipitation. Usually, July is the warmest month of the year, with an average temperature of 17.8 • C, while the lowest average temperatures occur in January, with −2.1 • C. Detailed seasonal courses of selected microclimatic variables during the year 2019, the year of study, are shown in Figure 2 and presented in the Results, Section 3.1.

Measurement of Microclimatic Parameters
Air temperature (T a ) and relative air humidity (RH) were measured 2 m above the soil surface using an EMS33 Rotro ® sensor (EMS, Brno, Czech Republic). Vapour pressure deficit (VPD) was calculated based on T a and RH values [44]. Radiation sensors EMS12 (EMS, Brno, Czech Republic) and CNR4 (Kipp & Zonen, Delft, The Netherlands) were used to measure the intensity of incoming photosynthetically active radiation (PAR) and global radiation (R g ), respectively. A tipping bucket rain gauge, Model 52202 (R.M. Young Company, Traverse City, MI, USA) was used to measure precipitation. Soil water content (SWC) in a depth of 30 cm was measured by a water content reflectometer, Model CS616L (Campbell Scientific, Logan, UT, USA). The signals from all sensors were recorded at 30 s intervals and stored as 30-min averages in the data logger RailBox V32P4 DX (EMS, Brno, Czech Republic).

Eddy Covariance Flux Measurements
Measurement of CO 2 , H 2 O, and O 3 fluxes between the ecosystem and the atmosphere was performed by the eddy covariance technique. The sensors, following the recommendations of pan-European infrastructure for flux measurements (ICOS), were mounted on a micrometeorological tower at 2.6 m above the soil surface, including a Gill R3 ® ultrasonic anemometer (Gill Instruments, Hampshire, UK), LI-7200 CO 2   All high-frequency eddy covariance data were processed by EddyPro 7.0.6 software and 30-min averaged flux values were calculated as: where F is the flux, measured in nmol m −2 s −1 , w represents the vertical wind component (m s −1 ), c is the measured concentration of a given substance (CO 2 , H 2 O, O 3 ), prime ( ) refers to the fluctuations and the overbar indicates an average over the time period. Averaged flux values were evaluated for data quality. During most of the day-time, the footprint was located within the crop field. If not, the data were discarded. For the flux calculation, the OpenEddy 0.0.0.9004 software [46] was used for thorough quality-checking, applying the standardized methodology similar to Mauder et al. [47]. The screening of the data included checking the plausible range, testing for integral turbulence characteristics [48], fetch filtering accepting only the data when more than 70% of the measured flux originated from the target area, and iterative low frequency spike removal [49]. These data quality tests are used by default within the ICOS community. The REddyPro software [50] was used to additionally filter half-hourly periods with insufficient turbulent mixing (i.e., u * filtering), gap-fill NEP (CO 2 flux) data, and partition NEP into GPP and ecosystem respiration R eco . For the u * filtering of NEP, the moving point method was applied [49,51]. To gapfill the missing NEP, the marginal distribution sampling [51] was used, where fluxes within a certain time window are binned to represent measure-ments during similar meteorological conditions [50]. Estimates of GPP were derived as the residual of R eco and NEP within so-called nighttime partitioning [51]. The methods used to perform each of these steps are discussed in detail in Wutzler et al. [50].

Modelling of Stomatal O 3 Flux
Stomatal (F sto ) and total O 3 flux (F tot ) were modeled from the aerodynamic (r a ) and canopy (r c ) resistances, and the resistance of the leaf boundary layer (r b ) [37]. F tot and F sto (nmol m −2 s −1 ) were calculated as: where c is [O 3 ] in a reference height of 2.6 m (ppbv), V d is a dry deposition velocity (1/(r a + r b + r c ); m s −1 ), and V dsto is the deposition velocity to stomata and mesophyll (m s −1 ). The model was enhanced to include stability corrections and distinguish between stable, neutral and unstable atmospheric conditions within the resistance r a [52]. The degree of atmospheric stability was further incorporated into the model. See the complete set of equations used and the detailed description of the entire model in the Supplementary Materials.

Modelling and Statistical Analyses
The effects of [O 3 ] and F sto on NEP were modelled by linear mixed-effects models using the lmer function from the 'lme4' package [53] in R 4.1.0 [54]. Half-hour eddy covariance and microclimatic data of the daytime period (06:00-18:30), when R g ≥ 10 W m −2 and the photosynthesis was thus sufficiently developed, were used for further analyses. The GPP together with [O 3 ] or F sto were used as NEP predictors. We assumed pre-industrial and future [O 3 ] of 12 and 60 ppbv, respectively, as reported by Cooper et al. [3] and Dentener et al. [4]. The corresponding F sto values are 0.038 and 1.438 nmol m −2 s −1 for pre-industrial and future conditions, respectively. Relationships between predictors, their interactions and NEP were tested for linearity. Only significant linear correlations were found. Predictors were included in the model in order to minimize AIC (Akaike information criterion) and p-values. Day of observation was used as a random factor. Package 'merTools' [55] provides methods for extracting results from mixed-effect model objects fit with the 'lme4' package. We used the function to predict the interval from this package to obtain the NEP and [O 3 ] or F sto relationship with its prediction bounds.

Microclimatic Growth Conditions
During the period of measurement (1 April-21 October 2019), the mean T a was 14.5 • C, RH was 72.9%, and SWC was 11.2%. Values of R g , PAR, T a , and VPD had typical annual courses, with values highest during summer months (June-August) ( Figure 2). The highest values of soil water content (15.0-17.0%) were observed at the beginning of the growing season (April) and the lowest (7.0%) in July.
The onset of the growing season was characterized by high values of mean daylight [O 3 ] ranging between 20 (overcast, rainy days) and 60 ppbv (clear sky conditions) before the first cut of the grass mixture. During the investigated period, the maximum and minimum values of 30-min averages of [O 3 ] were 63 and 27 ppbv, which occurred in late June and late October, respectively ( Figures 1D and 2).  Figure 3). In addition, the differences between daily maxima and night minima were negligible at the end of the vegetation season. Similarly, F sto showed typical daily courses characterized by their maximal and minimal values in clear sky summer days; however, they kept more or less constant during cloudy autumn days characterized by low values of VPD. Daily maxima of F sto were up to −1.5 and −2.0 nmol m −2 s −1 in spring and summer ( Figure 3A,B), respectively, but only up to −0.5 nmol m −2 s −1 in autumn ( Figure 3C). [O 3 ] and F sto values do not necessarily correlate as shown in Figure 1D; such inconsistencies likely correspond to stomatal closure preventing the diffusion of O 3 into the inner tissues.  Figure 3C).

NEP Modelling
The linear mixed-effects models were used to determine the relations between NEP and its predictors (irradiance, temperature, VPD, etc.) along with [O 3 ] and F sto . The NEP data, as half-hour averages, were modelled for the daytime period (06:00-18:30), when R g ≥ 10 W m −2 and air turbulence was always sufficiently developed. In addition, the model was used to calculate NEP based on fixed values of gross primary production (GPP) during the investigated growing period ( Figure S1A,B). The performance of the model is shown in Table 1.  The differences between NEP values measured by eddy covariance and those modelled by linear mixed-effects models at actual [O 3 ] and F sto were less than 1% (Table 2). Both models (based on [O 3 ] and/or F sto ) confirmed a substantial reduction in NEP induced by O 3 over the whole season when measured NEP values are compared to those modelled for the pre-industrial O 3 conditions. The estimated reduction of NEP ranges between 21.1% (May) and 38.8% (July) and between 13.1% and 26.2% when models based on [O 3 ] and F sto are applied, respectively.  To predict the future O 3 -induced changes in NEP, [O 3 ] and F sto of 60 ppbv and 1.438 nmol m −2 s −1 ) were applied in the model. Although the largest percentage decreases in NEP at high O 3 conditions were found in October, the effect is comparatively lower due to a small total carbon uptake at the end of the growing season. In contrast, the largest absolute reduction in NEP occurred in June. Noticeably, reductions predicted by the F sto -based model were lower as compared to those predicted by the [O 3 ]-based model over the entire growing season (Table 2).
Finally, values of a cumulative NEP (NEP cumul ), derived from the diurnal courses over the entire growing season, are shown in Figure 5. Values of NEP cumul were calculated as sums of daily NEP totals directly estimated by eddy covariance or modelled for the pre-industrial O 3 conditions from the beginning of the vegetation season (1 April 2019). In total, approximately 8.1 t of carbon was fixed per hectare over the growing season, when calculated from NEP over day hours (6:00-18:30), when the eddy covariance technique is generally more reliable. Additionally, the yield after the 1st, 2nd, and 3rd cut was 16.57, 6.7 and 3.65 t ha −1 in the form of green fodder, which corresponded to 8.27, 3.3 and 1.82 t ha −1 when dried. That corresponds to approximately 60.5% of C in dried biomass, when night-time respiration is not taken into account. Comparing to NEP cumul measured by eddy covariance at actual O 3 conditions, the model parametrized by low, pre-industrial [O 3 ] (12 ppbv) revealed higher values of NEP cumul , as much as 29.7%; however, it was only 13.5% when the model based on pre-industrial F sto (0.038 nmol m −2 s −1 ) was used.

Discussion
Previous studies have shown that tropospheric O 3 reduces carbon uptake in the forest [56] and grassland ecosystems [27,28]. Here, we aimed to quantify reductions in NEP of a managed ryegrass-clover mixture caused by current O 3 conditions and to predict these reductions under the O 3 conditions expected in the future. Indeed, we found an annual O 3 -induced reduction of NEP by 13.5-29.7% when compared to low pre-industrial O 3 conditions, depending on the construction of the predictive linear mixed-effects model. Accordingly, the accuracy of the models based on F sto and tropospheric [O 3 ] was compared. Our findings show that F sto is a better predictor of O 3 deposition and NEP decrease as compared to [O 3 ] (Figure 1D). High [O 3 ] does not necessarily result in high penetration of O 3 into plant tissues (F sto ; Figure 1D) due to closure of stomatal aperture over the noon hours, when [O 3 ] usually peaks [32]. Prediction of negative O 3 effects may be thus overestimated in [O 3 ]-based models. On the contrary, even low [O 3 ] of 3-5 ppbv may provide sufficient F sto inducing injuries of photosynthetic apparatus when stomata remain open [57].
Previous studies at a range of ecosystems found that microclimatic conditions at cloudy or semi-cloudy skies including optimum light intensity, T a and/or low VPD lead to a substantial stomatal openness and consequently high F sto [32,58,59]. Vice versa, high T a , light intensities, and VPD, together with low SWC, result in reduced stomatal conductance and F sto [60,61]. In accordance, we observed low F sto during the summer season (June-August) when VPD, T a , and light intensities are the highest together with high [O 3 ] and low SWC (Figures 2 and 3). On the contrary, F sto was highest in spring and autumn when VPD, T a , R g and PAR were lower, and SWC had increased and [O 3 ] remained lower than in summer. Such mechanism is likely to have been responsible for the relatively low stomatal conductance to O 3 diffusion and low O 3 -induced injuries in arid and hot South-European ecosystems as compared to high injuries of photosynthetic activity, plant growth and productivity observed in cold and wet North-European ecosystems [13,62]. However, severe drought periods can induce the hydro passive opening of the guard cells in Angiosperms and Gymnosperms [63], leading consequently to exacerbated O 3 damage [64]. Moreover, actual microclimatic conditions, particularly the seasonal course of VPD, have a strong effect on phenology and overall productivity of terrestrial ecosystems [65].
Tropospheric [O 3 ] has a typical seasonal course. In western Europe, [O 3 ] reaches its maximum in Spring, while, in Central and Eastern Europe, [O 3 ] is maximal in Summer [66]. We confirm the broad Summer maximum of [O 3 ], which was, however, not followed by the biggest NEP decline. Such a paradox is the result of low soil moisture (SWC) leading to a low stomatal conductivity and F sto . We found the highest O 3 -induced reduction of NEP in the Spring to Summer transition (June-July) compared to predicted values of NEP at pre-industrial [O 3 ] (Figure 4 and Table 2), i.e., the period of convenient conditions for stomatal openness and high F sto ( Figure 1D). This finding shows a strong interactive effect of [O 3 ] and microclimatic conditions on plant injuries [2].
We estimated the reduction of NEP, relative to pre-industrial conditions, to be between 13.1% in May and 26.2% in July. Previous studies indicated that the O 3 -induced decrease in NEP is associated with impaired photosynthetic CO 2 uptake, root respiration, and stomatal conductance as well as premature senescence of plants [34]. Such changes may further lead to a reduced biomass accumulation, yield quantity and crop quality [67]. Reduced photosynthesis is particularly associated with oxidized fatty acids of bio-membranes, enhanced production of reactive oxygen species (ROS) resulting in the degradation of chloroplasts, photosynthetic pigments and morphological changes of mesophylls [12]. In a meta-analysis of O 3 -induced effects on crops, Feng [69], with the greatest reduction in summer, along with the peak of tropospheric [O 3 ]. That is in line with our observation when the largest reduction of NEP in July was found.
There is a huge variety in plant sensitivity to [O 3 ]. Such variability is likely to be the consequence of the capacity of the detoxification processes that are genetically based [70]. Studying Trifolium and Lolium species, Hayes et al. [22] reported a decrease of biomass by 33% in T. repens due to O 3 exposure, but no effect on L. perenne grown either in a mixture or monoculture (no change of canopy structure was observed when the species were grown in monoculture or mixture [22]. A mixture of T. pratense and Phleum pratense under elevated [O 3 ] showed similar results with a decrease in biomass yield of T. pratense [71]. That suggests higher capacity of L. perenne for O 3 removal and perhaps it out-competes T. repens in elevated [O 3 ], leading to changes in species composition when grown in mixtures. The possible unbalanced reduction has a significant effect on pasture communities, where those species are more-or-less co-dominant, because it allows invasion of other O 3 tolerant species. In addition, a reduction of stolon biomass has been reported in T. pratense, which results in a higher reduction in later grassland cuttings due to a carry-over effect. New leaves are produced from nodes of stolons, therefore the ability of the plant to spread and build up new biomass is reduced [22]. A further reduction has been observed due to a decrease of nodulation in T. repens when grown in monoculture [72] and mixture [73]. That further affects N-fixing capacity and diversion of assimilates away from the energetically costly process of N-fixation and nodule growth. Therefore, we assume that in our experimental site, L. perenne partially balanced the larger reduction of T. pratense biomass. A broad range of field studies has shown that precise estimation of F sto is crucial for calculation of ozone deposition velocity, and thus of F sto . Therefore, precise estimation of F sto is needed. In our case, the model error is below 1%, which is far below the varying differences in O 3 deposition rates of 13% found by Otu-Larbi et al. [74]. However, it should be noted that our model predicts pre-industrial NEP based on fixed value of [O 3 ] and F sto over the day, which has no effect on the prediction due to the linearity of the model. Nevertheless, having diurnal variation in those variables will be at least logically more relevant, although modelling the shape of the diurnal profile during the pre-industrial era will be difficult due to the lack of diurnal observations in the past. We found that the model based on [O 3 ] shows the reduction in cumulative measured NEP at low, pre-industrial [O 3 ] by as much as 29.7% against a reduction of 13.5% with the model based on F sto ( Figure 5). In line with our finding, F sto is recommended to be applied for modelling studies and estimations of vegetation damage [75]. Nevertheless, European legislation is still based on [O 3 ] instead of F sto (European Council Directive 2008/50/EC). Moreover, only [O 3 ] is monitored at most European localities [76]. Therefore, studies comparing both approaches are needed.

Conclusions
We quantified the reductions in NEP of a managed ryegrass-clover mixture caused by current O 3 conditions compared to pre-industrial O 3 conditions and we predicted these reductions under expected future O 3 conditions. Reductions of NEP were estimated by a linear mixed-effects model based on O 3 concentration ([O 3 ]) and stomatal O 3 flux (F sto ). We found the F sto -based model more precise, showing high F sto and, consequently, substantial NEP reduction even during low [O 3 ] periods. On the contrary, the broad summer maximum of [O 3 ] is not followed by the biggest NEP decline, due to low soil water content leading to a low stomatal conductivity and F sto . Compared to the lower pre-industrial O 3 conditions, current levels of O 3 result in a reduction of cumulative NEP over the entire growing season of up to 29.7 and 13.5% when the [O 3 ]-based and F sto -based model is applied, respectively. The intra-annual reduction varies between 13.1% in May and 26.2% in July in the F sto model. These findings indicate the importance of F sto and its inclusion into the models estimating O 3 effects on terrestrial vegetation.