Early Detection of Sage (Salvia officinalis L.) Responses to Ozone Using Reflectance Spectroscopy

Advancements in techniques to rapidly and non-destructively detect the impact of tropospheric ozone (O3) on crops are required. This study demonstrates the capability of full-range (350–2500 nm) reflectance spectroscopy to characterize responses of asymptomatic sage leaves under an acute O3 exposure (200 ppb for 5 h). Using partial least squares regression, spectral models were developed for the estimation of several traits related to photosynthesis, the oxidative pressure induced by O3, and the antioxidant mechanisms adopted by plants to cope with the pollutant. Physiological traits were well predicted by spectroscopic models (average model goodness-of-fit for validation (R2): 0.65–0.90), whereas lower prediction performances were found for biochemical traits (R2: 0.42–0.71). Furthermore, even in the absence of visible symptoms, comparing the full-range spectral profiles, it was possible to distinguish with accuracy plants exposed to charcoal-filtered air from those exposed to O3. An O3 effect on sage spectra was detectable from 1 to 5 h from the beginning of the exposure, but ozonated plants quickly recovered after the fumigation. This O3-tolerance was confirmed by trends of vegetation indices and leaf traits derived from spectra, further highlighting the capability of reflectance spectroscopy to early detect the responses of crops to O3.


Introduction
Tropospheric ozone (O 3 ), produced by a variety of precursors, such as nitrogen oxides (NO x ) and volatile organic compounds (VOCs), under light conditions, is a major phytotoxic air pollutant, with deleterious effects also on animal health [1][2][3]. Although numerous attempts have been made to attenuate the emission of its precursors, levels of O 3 are still high in many areas of the world (e.g., North America, Europe, and South and East Asia) and are expected to increase further due to the presence of climatic changes and to anthropogenic activities [4], especially in hot-spot areas such as the Mediterranean basin [5]. However, decreases in surface O 3 concentrations have been also reported in some Mediterranean sites [6]. The deleterious effects of excessive O 3 uptake on plants include reduction of photosynthesis and growth, partial stomatal closure, cell dehydration, excessive excitation energy, accelerated leaf senescence, and appearance of chlorotic/necrotic leaf injuries, that overall result in reduction of crop production and huge economic loss [1]. Although O 3 is commonly debated as a major air pollutant with detrimental effects on plants, some studies have also proposed exposure to an properties. In particular, sage contains antioxidant molecules that play an important role against human diseases [32,33]. This is a species with good level of tolerance against salinity [34] and drought [35]. Recent studies have also demonstrated that sage is able to activate biochemical compounds, including phytohormones, signaling molecules, and antioxidants, in order to overcome O 3 -induced oxidative pressure, and for that it can be considered as O 3 tolerant [36,37].
Here, we test the capability of reflectance spectroscopy to rapidly and non-destructively monitor the responses of asymptomatic sage plants exposed to an acute O 3 exposure. Specifically, the aims of this work were (i) to develop spectroscopic models for the estimation of a variety of traits related to photosynthesis, the oxidative pressure induced by O 3 , and the antioxidant mechanisms adopted by plants to cope with the pollutant, (ii) to evaluate the potential of spectral phenotyping to accurately detect and predict O 3 stress, and (iii) to elucidate the variations of both vegetation indices and the abovementioned leaf traits predicted from spectra by PLSR-models, under O 3 exposure.
Predictive models very accurately characterized gas exchange traits with average model goodness-of-fit (R 2 ) for validation of 0.65 for A, 0.80 for E, 0.70 for g s , 0.72 for C i , 0.90 for WUE i , 0.76 for WUE in , 0.67 for k, and 0.80 for T l . The RMSEs for validation data were from 1.4 to 3.6 fold higher than RMSEs for calibration data, while larger differences were for WUE i and WUE in (6.2 and 16 fold). Lower prediction performance metrics were found for biochemical traits with average R 2 for validation of 0.65 for MDA, 0.71 for ORAC, 0.50 for HORAC, 0.45 for DHA, 0.63 for DHA/ASA TOT , 0.58 for GSH, 0.67 for GSH TOT , 0.53 for Chl a, 0.42 for Chl TOT , 0.59 for Car, 0.61 for Phen. The RMSEs for validation data were from 1.2 to 2.6 fold higher than RMSEs for calibration data, while larger differences were for DHA/ASA TOT , Car and Phen (6.0, 4.0, and 4.6 fold). Other performance statistics of PLSR-models are shown in Table 1 and  Standardized coefficient and variable important to the projection (VIP) statistic profiles highlighted important wavelengths for prediction from 450 to 750 nm for A and k, around 700, 1450, and 1900 nm for E, WUE i , WUE in , and T l , and around 1400 and 1900 nm for g s and C i (Figures 1 and 2). For biochemical traits, standardized coefficients and VIP were most pronounced around 700 and 1900 nm for MDA, around 1400 and 1900 nm for ORAC, around 700, 1400, and 1900 nm for HORAC, around 400 and 700 nm for GSH and GSH TOT , and around 700 nm for Chl a and Chl TOT . Given the large number of peaks observed in standardized coefficients and VIP profiles of DHA, DHA/ASA TOT , Car and Phen, it was not possible to highlight specific wavelengths more important than others for predictions of these traits (Figures 3-5). Figure 1. (a,c,e,g) Observed versus partial least squares regression (PLSR)-predicted values of CO2 assimilation rate (A), transpiration (E), stomatal conductance (gs), and intercellular CO2 concentration (Ci) in sage; error bars for predicted values represent the standard deviations generated from 500 simulated models; dashed line is 1:1 relationship; model goodness-of-fit (R 2 ), root mean square error (RMSE), bias and %RMSE for validation data generated using 80% of the data for calibration and 20% for validation are reported. (b,d,f,h) Mean (solid), 5th and 95th percentile (dotted) of standardized coefficients (black), and variable importance for projection values (VIP, blue) by wavelengths for PLSR-models predicting A, E, gs, and Ci in sage. Figure 1. (a,c,e,g) Observed versus partial least squares regression (PLSR)-predicted values of CO 2 assimilation rate (A), transpiration (E), stomatal conductance (g s ), and intercellular CO 2 concentration (C i ) in sage; error bars for predicted values represent the standard deviations generated from 500 simulated models; dashed line is 1:1 relationship; model goodness-of-fit (R 2 ), root mean square error (RMSE), bias and %RMSE for validation data generated using 80% of the data for calibration and 20% for validation are reported. (b,d,f,h) Mean (solid), 5th and 95th percentile (dotted) of standardized coefficients (black), and variable importance for projection values (VIP, blue) by wavelengths for PLSR-models predicting A, E, g s , and C i in sage. Figure 2. (a,c,e,g) Observed versus partial least squares regression (PLSR)-predicted values of intrinsic water use efficiency (WUEi), instantaneous water use efficiency (WUEin), CO2 assimilation rate:intercellular CO2 concentration ratio (k) and temperature of adaxial leaf surface (Tl) in sage; error bars for predicted values represent the standard deviations generated from 500 simulated models; dashed line is 1:1 relationship; model goodness-of-fit (R 2 ), root mean square error (RMSE), bias and %RMSE for validation data generated using 80% of the data for calibration and 20% for validation are reported. (b,d,f,h) Mean (solid), 5th and 95th percentile (dotted) of standardized coefficients (black), and variable importance for projection values (VIP, blue) by wavelengths for PLSR-models predicting WUEi, WUEin, k, and Tl in sage. (a,c,e,g) Observed versus partial least squares regression (PLSR)-predicted values of intrinsic water use efficiency (WUE i ), instantaneous water use efficiency (WUE in ), CO 2 assimilation rate:intercellular CO 2 concentration ratio (k) and temperature of adaxial leaf surface (T l ) in sage; error bars for predicted values represent the standard deviations generated from 500 simulated models; dashed line is 1:1 relationship; model goodness-of-fit (R 2 ), root mean square error (RMSE), bias and %RMSE for validation data generated using 80% of the data for calibration and 20% for validation are reported. (b,d,f,h) Mean (solid), 5th and 95th percentile (dotted) of standardized coefficients (black), and variable importance for projection values (VIP, blue) by wavelengths for PLSR-models predicting WUE i , WUE in , k, and T l in sage.       Table 1. Range, number of latent variables (LV), and model goodness-of-fit (R 2 ), root mean square error (RMSE), and bias for calibration (cal) and validation (val; i.e., root-mean square error for prediction) data generated using 500 random permutations of the data with 80% used for cal and 20% used for val for PLSR-models predicting leaf traits from sage spectra. Data are shown as mean ± standard deviation. Bias for cal is not shown as it was always <0.001. Trait abbreviations: A, CO 2 assimilation rate (µmol m −2 s −1 ); E, transpiration (mmol m −2 s −1 ); g s , stomatal conductance (mol m −2 s −1 ); C i , intercellular CO 2 concentration (µmol mol −1 ); WUE i , instantaneous water use efficiency (µmol mmol −1 ); WUE in , intrinsic water use efficiency (µmol mol −1 ); k, instantaneous carboxylation efficiency (µmol m −2 s −1 Pa −1 ); T l , temperature of adaxial leaf surface ( • C); MDA, malondialdehyde (µmol g −1 DW); ORAC, oxygen radical absorption capacity (µmol TE g −1 DW); HORAC, hydroxyl radical antioxidant capacity (µmol GAE g −1 DW); DHA, oxidized ascorbate (mg g −1 DW); DHA/ASA TOT , oxidized:total ascorbate ratio; GSH, reduced glutathione (µmol g −1 DW); GSH TOT , total glutathione (µmol g −1 DW); Chl a, chlorophyll a (mg g −1 DW); Chl TOT , total chlorophyll (mg g −1 DW); Car, carotenoids (mg g −1 DW); Phen, total phenols (mg GA g −1 DW). An external validation of PLSR-models is suggested to further test their accuracies in estimation of leaf traits.

Analyses of Spectral Signatures
A number of different spectral ranges (Table S2) were examined to optimize the statistical outputs (i.e., to get highest significances for the tested effects on reflectance profiles) of the permutational multivariate analysis of variance (PERMANOVA), but the best results were obtained using the full range (i.e., 400-2400 nm). Final PERMANOVA revealed that only O 3 affected reflectance profile of sage, whereas no significant effects were found for time and O 3 × time (Table 2). However, keeping times of analysis separated, PERMANOVA showed significant O 3 effects only at 1, 2, and 5 h from the beginning of exposure (FBE), as shown by Figure 6 related to the principal coordinates analysis (PCoA). By partial least squares discriminant analysis (PLS-DA), the best classification of control versus O 3 plant groups from spectra (higher kappa) was found with a 80:20 ratio for calibration:validation data using nine components (accuracy: 0.88; 95% confidential interval of accuracy: 0.59-0.98; kappa: 0.74 ± 0.17). According to PERMANOVA, keeping times of analysis separated, good classifications between groups were found at 1, 2, and 5 h FBE (average accuracy: 0.83, 0.72, and 0.75; kappa: 0.69, 0.40, 0.48, respectively), whereas lower classification outputs were found at 0, 8, and 24 h FBE (average accuracy: 0.52, 0.64 and 0.43; kappa: 0.04, 0.16, −0.13, respectively) ( Table 3). Table 2. F-and p-values of two-way permutational multivariate analysis of variance (PERMANOVA) for the effects of ozone, time, and their interaction on full-range (400-2400 nm) reflectance profiles of sage leaves. df represents the degrees of freedom. *** p ≤ 0.001, ns: p > 0.05.   Table 3) are shown on the bottom-left corners of panels. Table 4. F-and p-values of one-way repeated measures analysis of variance (ANOVA) for the effects of ozone (between factor), time (within factor), and their interaction on spectral indices and leaf traits derived from sage spectra (only leaf traits from well performing PLSR-models were used: validation R 2 > 0.70 and 0.60 for gas exchange and biochemical traits, respectively). df represents the degrees of freedom. *** p ≤ 0.001, ** p ≤ 0.01, * p ≤ 0.05; ns: p > 0.05. See Table 1 Table 3) are shown on the bottom-left corners of panels. Table 3. Latent variables (LV), average accuracy (Acc), 95% confidential interval (CI) of accuracy, and kappa (K; mean ± standard deviation) for validation data generated using partial least squares discriminant analysis (PLS-DA), using 80% of the data for calibration and 20% for validation, for the classification of control versus ozonated plant groups from full-range (400-2400 nm) sage spectra, at 0, 1, 2, 5, 8, and 24 h from the beginning of exposure.

Variations of Vegetation Indices and Leaf Traits Derived from Spectra
The effects of O 3 , time, and their interaction on vegetation indices and selected leaf traits derived from spectra are shown in Table 4. Only leaf traits derived by well-performing PLSR-models were selected (validation R 2 > 0.71 and 0.60 for gas exchange and biochemical traits, respectively). The interaction O 3 × time and time alone were significant on PRI, WUE i , WUE in , MDA, DHA/ASA TOT , and GSH TOT . A significant effect of time alone was only found for C i and Phen, while both O 3 and time alone were significant on ORAC, but not their interaction. Normalized Difference Vegetation Index, CI, and sPSRI did not show any significance.
Both WUE i and WUE in were reduced by O 3 exposure at 2 and 5 h FBE (around 33 and 28% in comparison to controls, respectively), while PRI were reduced by the pollutant only at 5 h FBE (−46%); but these traits then recovered at the following times of analysis (Figure 7a

Discussion
The present study shows an accurate and non-destructive approach by which plant responses to O3 can be rapidly detected and monitored using reflectance spectroscopy. Firstly, by combining reflectance measurements, standard physiological and biochemical analyses, and robust statistical modelling, this study demonstrated the potential to concomitantly predict numerous widely used leaf traits related to crop-O3 interaction from spectral data.
Photosynthesis is a major plant process that is severely affected by O3, primarily through stomatal and mesophyll limitations [1]. Standard measurements of light-and CO2-saturated photosynthesis, however, are sometimes logistically challenging, potentially taking several minutes per leaf as the leaf has to acclimate to the cuvette conditions. Spectral approaches present a valid alternative to standard reference measurements and have been used previously to estimate photosynthetic activity in plants both indirectly, via xanthophyll cycling (PRI, [12,38]) and directly, by predicting specific traits from spectral data, mainly RuBP, Vcmax, and Jmax (e.g., [39][40][41][42]). Interestingly, we found excellent prediction performance for all the modeled gas-exchange traits (A, E, gs, Ci, WUEi, WUEin, k, and Tl), with validation R 2 ranging from 0.65 to 0.90 (an external validation of these PLSR-models would be suggested to further test their estimation accuracy, especially for those with larger differences between RMSEs for calibration and validation data). According to the coefficient and VIP profiles, pigment-(450-750 nm) and water-related wavelengths (around 1400 and 1900 nm) were particularly important for these accurate estimations. However, best prediction accuracies (R 2 > 0.75) were found for models using the full range (400-2400; i.e., E, WUEi, WUEin and Tl), followed by those using a narrower range (950-2400 nm; i.e., gs and Ci), and then by those using  Table 4. F-and p-values of one-way repeated measures analysis of variance (ANOVA) for the effects of ozone (between factor), time (within factor), and their interaction on spectral indices and leaf traits derived from sage spectra (only leaf traits from well performing PLSR-models were used: validation R 2 > 0.70 and 0.60 for gas exchange and biochemical traits, respectively). df represents the degrees of freedom. *** p ≤ 0.001, ** p ≤ 0.01, * p ≤ 0.05; ns: p > 0.05. See Table 1

Discussion
The present study shows an accurate and non-destructive approach by which plant responses to O 3 can be rapidly detected and monitored using reflectance spectroscopy. Firstly, by combining reflectance measurements, standard physiological and biochemical analyses, and robust statistical modelling, this study demonstrated the potential to concomitantly predict numerous widely used leaf traits related to crop-O 3 interaction from spectral data.
Photosynthesis is a major plant process that is severely affected by O 3 , primarily through stomatal and mesophyll limitations [1]. Standard measurements of light-and CO 2 -saturated photosynthesis, however, are sometimes logistically challenging, potentially taking several minutes per leaf as the leaf has to acclimate to the cuvette conditions. Spectral approaches present a valid alternative to standard reference measurements and have been used previously to estimate photosynthetic activity in plants both indirectly, via xanthophyll cycling (PRI, [12,38]) and directly, by predicting specific traits from spectral data, mainly RuBP, V cmax , and J max (e.g., [39][40][41][42]). Interestingly, we found excellent prediction performance for all the modeled gas-exchange traits (A, E, g s , C i , WUE i , WUE in , k, and T l ), with validation R 2 ranging from 0.65 to 0.90 (an external validation of these PLSR-models would be suggested to further test their estimation accuracy, especially for those with larger differences between RMSEs for calibration and validation data). According to the coefficient and VIP profiles, pigment-(450-750 nm) and water-related wavelengths (around 1400 and 1900 nm) were particularly important for these accurate estimations. However, best prediction accuracies (R 2 > 0.75) were found for models using the full range (400-2400; i.e., E, WUE i , WUE in and T l ), followed by those using a narrower range (950-2400 nm; i.e., g s and C i ), and then by those using an even smaller region (400-1000 nm; i.e., A and k). This outcome suggests that the use of wider ranges is a good practice to obtain better predictions of leaf traits from spectra, since (i) using larger ranges means to incorporate more signals of chemical, morphological, and physiological properties of leaves included in the spectra, and (ii) using more prediction variables allows to increase the number of components adopted for modelling without overfit the models.
However, the use of narrower ranges containing only specific absorption wavelengths for the trait to estimate sometimes leads to better predictions than using the full range, since the incorporation of other spectral regions may reduce the prediction ability of those trait-specific wavelengths [27]. Thus, it is not surprising that the best predictions for A and k were found using only the wavelengths from 400 to 1000 nm, with coefficient and VIP profiles highlighting the importance of the region 450-750 nm, since this range includes leaf pigment absorption features [43], as well as the red-edge (700-750 nm; [44]). The importance of these spectral features in the assessment of the photosynthetic processes has been previously reported for several plant species (e.g., [12,39,42,43]), and numerous studies have also shown that the shape of the red-edge is dependent on chlorophyll content (e.g., [45,46]) and stress conditions (e.g., [10,44]). The spectral region best predicting g s and C i was instead in the NIR-SWIR, from 950 to 2400 nm, which is dominated by water content and outside of wavelengths commonly associated with pigments [27]. Thus, the ability to estimate stomatal behavior from spectral data was likely due to the sensitivity of spectra to g s -related water regulations. The ability to estimate C i was instead likely ascribable to specific CO 2 -absorption features included in the NIR-SWIR [47], although also this trait could be indirectly associated to water regulations since g s and C i are usually strongly related. These interpretations are supported by the profiles of coefficients and VIP statistics, showing that wavelengths around 1400 and 1900 nm were extremely important for predictions of g s and C i . The best ranges found in the present study for predicting gas-exchange traits are not in agreement with our previous findings on maize under drought conditions (best predictions occurred using the 500-900 nm spectral range [48]), suggesting that these outcomes are likely species-and environmental-specific.
Stomatal closure is considered the first barrier against O 3 since it limits uptake and preserves plants from oxidative damage, but plants have also evolved enzymatic (e.g., superoxide dismutase, catalase, glutathione peroxidase) and non-enzymatic (e.g., ascorbate, glutathione, phenols, carotenoids) antioxidant systems to cope with considerable amount of reactive oxygen species (ROS) generated by the reaction of O 3 with the leaf cell apoplast [1,49]. The characterization of these plant-O 3 interactions by traditional biochemical analyses may be precise, but has several limitations since these analyses are destructive, time consuming and expensive, all aspects that make these methods logistically impractical for monitoring a large number of individual plants. Here, we developed PLSR-models to predict from spectral data a number of widely used leaf traits related to the lipid peroxidation induced by O 3 (MDA) as well as the antioxidant capacity (ORAC and HORAC) and the main non-enzymatic antioxidants of plants (DHA, DHA/ASA TOT , GSH, GSH TOT , Car, and Phen), in addition to Chl a and Chl TOT having already been largely investigated. However, prediction accuracies for biochemical traits were lower than those reported for gas-exchange parameters (validation R 2 ranging from 0.42 to 0.71), especially considering the higher amount of outliers removed before running final models: only PLSR-models for MDA, ORAC, DHA/ASA TOT , GSH TOT , and Phen showed acceptable performances (validation R 2 > 0.60), while PLSR-models developed to estimate HORAC, DHA, GSH, and Chl a and Chl TOT showed even lower accuracies (validation R 2 < 0.60). Outputs for ORAC and chlorophyll traits, the few parameters already tested by previous studies, were not in accordance with previous reports since Yendrek et al. [42] found low accuracy for a PLSR-model developed to predict ORAC from maize spectra, whereas the ability to estimate chlorophyll contents from spectral data has been demonstrated since long [43]. The inconsistencies between the present study and previous investigations confirm the requirement to develop different models for specific species and environments, as well as to improve the protocols adopted for modeling. Although different ranges were used to best predict the biochemical traits (400-2400 nm for HORAC, GSH TOT , and Chl TOT ; 400-750 + 1400-2400 for MDA; 1400-2400 nm for ORAC; 1100-1800 nm for DHA and DHA/ASA TOT ; 400-900 nm for GSH; 400-700 + 1600-1800 nm for Chl a; 1100-1400 nm for Car; and 1100-1600 nm for Phen), the regions from 400-700 nm and around 1400 and 1900 nm resulted the most important for estimations, as shown by coefficient and VIP profiles. This is in accordance with previous reports about remote sensing of foliar chemistry [10,18,47,50]. Although the approach here proposed may have limitations discriminating fine scale differences in biochemical traits (and we thus encourage caution when interpreting results from a narrow range of values) and an external validation of these PLSR-models (as well as of those for ecophysiological traits) would be suggested to further test their estimation accuracy, the present study also demonstrates the potential to expand prediction capabilities of spectral data for key leaf biochemical features involved in the response of plants to O 3 pressure.
The prediction of these outcomes, in combination with analyses of spectral signatures, has the potential to provide multiple layers of stress-specific information to growers, including the characterization of the underlying physiological and biochemical responses to O 3 but also the rapid identification of stress conditions, that can increase the efficiency of management practices. Effectively, the second major achievement of the present study was the demonstration that reflectance profiles of sage plants are sensitive to acute O 3 exposure. In the absence of visible symptoms, we were able to distinguish with high accuracy sage plants exposed to charcoal-filtered air from those exposed to O 3 using the full range spectral region (i.e., 400-2400). The utilization of full range spectral data (i.e., a phenotypic expression of the overall physiological, morphological, and biochemical characteristics of leaves under specific environmental conditions; [28]) may prevent loss of information and represent a more powerful approach than direct measurements of leaf traits to monitor crop status. However, the two-way O 3 × time interactive effect on spectra profiles was not detected. Better outputs might be reached by increasing the experimental/plant replications; this is especially true in field settings, where growing conditions can be highly variable. However, keeping times of analysis separated, we were able to distinguish an O 3 effect on sage already at 1 h FBE (when none of the traits reported below showed a significant effect), as well as the next two times of analysis, until the end of the O 3 -exposure (i.e., 2 and 5 h FBE). Also, these discriminations resulted in highly accurate results, highlighting the capability of this approach to detect early O 3 stress, as previously reported for several abiotic and biotic stressors (e.g., [51,52]).
By a phytotoxicologic point of view, the analyses of spectral signatures suggested that O 3 induced some physiological changes in sage leaves during the exposure (from 1 to 5 h FBE), but plants quickly recovered after the fumigation, confirming the tolerance of this species to the pollutant [36,37]. These responses were confirmed by variations of the investigated vegetation indices and leaf traits derived from spectra (again, only traits from best performing PLSR-models were used), further highlighting the potential of reflectance spectroscopy in monitoring the responses of plants to O 3 . The photosynthetic performance was impaired by the pollutant at 2 and 5 h FBE but recovered to optimal levels after the end of exposure, as confirmed by the trends in WUE i , WUE in , and PRI (C i was not affected by O 3 ). The unchanged values of NDVI, CI, and sPSRI suggest that a senescing process was not induced. Indeed, slightly higher levels of lipid peroxidation (i.e., increased MDA) were only observed at the end of the exposure, and the antioxidant capacity was overall increased by the pollutant (i.e., increased ORAC). The antioxidant response seemed to be finely regulated by the activation/suppression of specific antioxidants at the different times of analysis, both during and after the exposure. A key role in the antioxidant response seemed to be played by the Halliwell-Asada cycle [53], given the O 3 -induced reduction of the DHA/ASA TOT ratio observed at 8 h FBE and the induction of GSH TOT at 5 h FBE; whereas the phenolic metabolism was not likely involved in the response. However, other enzymatic and non-enzymatic antioxidants should be investigated to better highlight how sage cope the O 3 -induced oxidative pressure [49]. Overall, the responses showed by sage to tolerate the acute O 3 exposure are in accordance with several previous findings (e.g., [37,54,55]).
In conclusion, this study confirms that reflectance spectroscopy can provide a rapid, non-destructive, and relatively inexpensive approach to accurately and concomitantly quantify a number of physiological and biochemical traits in crops associated to the oxidative pressure induced by O 3 using a single spectral measurement. In addition, it shows that spectral information can successfully identify stress conditions in crops exposed to O 3 in absence of visible symptoms and using an O 3 -tolerant species, as confirmed by the investigated physiological and biochemical leaf traits. Outcomes and approaches presented in the current study could be of application in a number of scientific fields (e.g., precision agriculture, high-throughput plant phenotyping) with benefits not only for further plant science research, but also for growers to achieve greater crop yield and quality (and incomes), with lower environmental impact. Furthermore, this approach may help to monitor plant function over large geographic regions given its potential capability of being scaled to remote collection from air-or space-borne platforms [10].

Plant Material and Experimental Design
Experimental activities were performed at the field-station of San Piero a Grado (Pisa, Italy; 43 • 40 48 N, 10 • 20 46 E, 2 m a.sl.) owned by the Department of Agriculture, Food and Environment of the University of Pisa. On September 2018, thirty-six sage seedlings (4 month-old, 10-15 stems per seedling), grown under field conditions in plastic pots (3.7 L volume) containing a mix of peat and steam-sterilized soil (1:1, v/v), were selected for uniformity in size (approximately 25 cm tall), equally distributed among four controlled environment fumigation chambers placed inside a walk-in growth chamber, and kept under 22 ± 1 • C of temperature, 85 ± 5% of relative humidity (RH), 500 µmol m −2 s −1 of photosynthetic active radiation (PAR) provided by incandescent lamps during a 12 h photoperiod, and charcoal-filtered air. The fumigation system was continuously ventilated (two complete air changes per min) with charcoal-filtered air. After two weeks, half of the plants were exposed to 200 ppb of O 3 (1 ppb = 1.96 µg m −3 , at 25 • C and 101.325 kPa) for 5 h from 10:00 to 13:00, while the other half of the plants were maintained under charcoal-filtered air (controls, O 3 concentration < 5 ppb; two chambers per O 3 -treatment). Ozone was generated by a Fisher 500 air-cooled generator (Fisher America Inc., Houston, TX, USA), supplied with pure oxygen, and mixed with the inlet air entering the fumigation chamber; its concentration in the fumigation chambers was continuously monitored with a Serinus 10 analyzer (Ecotech Acoem Group, Milan, Italy). The O 3 concentration in fumigation chambers used for the high O 3 treatment was <5 ppb after the end of the exposure. The O 3 exposure was performed according to Cotrozzi et al. [56] Investigations were performed at 0, 1, 2, 5, 8, and 24 h from the beginning of exposure (FBE). Six control plants and six to ten O 3 plants were dedicated to the analyses of spectral signatures and to the final estimations of leaf traits by PLSR-models developed (see below), thus repeatedly measured, within few minutes and at each time of analysis, only for leaf reflectance and never harvested. The remaining plants were instead dedicated to the development of PLSR-models and, thus, as soon as the previous measures (i.e., those dedicated to the analyses of spectral signatures) were completed, they were consecutively measured for leaf-gas exchange and reflectance and finally harvested (see below). Some of these plants were processed (i.e., gas-exchange-reflectance-harvest) more times (to a maximum of four consecutive times), in order to obtain a total number of 81 gas-exchange measurements and samples for biochemical analyses, coupled with reflectance collections. A summary of the measurement design is reported in Table S1. Gas exchange measurements were performed on one of the second mature and fully expanded leaves (counting from the top of stems), whereas similar and contemporary leaves (three per plant) were harvested, stored at −80 • C, and later freeze-dried for biochemical analyses. Gas exchange analyses made on plants addressed to the development of PLSR-models and measured more times were performed on distinct leaves at each time of analysis. The onset of visible foliar symptoms was constantly monitored.

Collection of Leaf Spectra
Reflectance profiles of leaves were collected using a full range (350-2500 nm) ASD FieldSpec 4 spectroradiometer (Analytical Spectral Devices, Boulder, CO, USA), equipped with a leaf-clip including an internal halogen light source attached to a plant probe. Measurements were made on two areas of the adaxial surface of each leaf, with one measurement per area, and collections were combined to produce an average leaf spectrum. The relative reflectance of each leaf was determined from the measurement of leaf radiance divided by the radiance of a white reference panel included in the leaf-clip, measured every 16 spectral collections.

Foliar Gas-Exchange
The CO 2 assimilation rate (A), transpiration (E), stomatal conductance (g s ), intercellular CO 2 concentration (C i ), and leaf temperature (T l ) were determined using a LI-6400 portable photosynthesis system equipped with a 2 × 3 cm chamber and a 6400-02B LED light source (Li-COR Inc., Lincoln, NE, USA), operating at 400 ppm CO 2 concentration, 25 ± 2 • C of leaf temperature, 45 ± 5 % of RH, 1.8 ± 0.2 kPa of VPD and saturating light conditions (1500 µmol m −2 s −1 PAR). Instantaneous (WUE i ) and intrinsic water use efficiencies (WUE in ) were calculated as A/E and A/g s , respectively. Instantaneous carboxylation efficiency (k) was calculated as A/C i .

Biochemical Analyses
Lipid peroxidation was measured by the TBARS (thiobarbituric acid reactive substances) method, according to Guidi et al. [57]. Briefly, 30 mg of leaf samples were extracted with 750 µL of 80% ethanol, sonicated three times for 10 min and centrifuged at 13,000× g for 10 min at 4 • C. Then, 100 µL of each sample supernatant were mixed with 400 µL of 20% trichloroacetic acid (TCA) and 0.5% thiobarbituric acid (TBA). Samples were incubated at 95 • C for 30 min, and centrifuged at 12,000× g for 10 min at The antioxidant properties were tested by measuring the oxygen radical absorption capacity (ORAC) and hydroxyl radical antioxidant capacity (HORAC) according to Ou et al. [58] and Ou et al. [59], respectively. For both analyses, 10 mg of leaf samples, were extracted with 1 mL of 100% methanol, sonicated three times for 10 min, incubated over-night at 4 • C and centrifuged at 13,000× g for 15 min at 4 • C. Then, 10 µL of each sample supernatant were mixed with 170 µL of 200 µM fluorescein and incubated at 37 • C for 20 min. The antioxidant scavenging activity was induced by 2,2 -azobis-(2-amidino-propane) dihydrochloride (AAPH) and Co(II) complex, respectively. The antioxidant activity was quantified with excitation at 480 nm and emission at 530 nm, using a Victor3 1420 Multilabel Counter (Perkin Elmer, Waltham, MA, USA). Fluorescence/absorbance was converted to concentration based on an antioxidant standard curve expressed in Trolox (T) equivalents for ORAC and gallic acid equivalents for HORAC.
The total amount of ascorbic acid (ASA TOT ) and its oxidized form (AsA) were determined according to Kampfenkel et al. [60]. Around 10 mg of leaf samples were extracted with 1 mL of 6% TCA, sonicated three times for 10 min and centrifuged at 12,000× g for 10 min at 4 • C. This assay is based on the reduction of ferric ion (Fe 3+ ) to ferrous ion (Fe 2+ ) with AsA in acid solution, followed by formation of a red chelate between Fe 2+ and 2,2-dipiridil. Samples were finally read for absorbance at 525 nm with the same spectrophotometer reported above. Absorbances were converted to concentration based on an ascorbic acid standard curve. Concentrations of oxidized ascorbate (DHA) were calculated as ASA TOT − AsA. Oxidized:total ascorbate ratio (DHA/ASA TOT ) was also calculated.
Total amount of total glutathione (GSH TOT ) and its oxidized form (GSSG) were determined according to Sgherri and Navari-Izzo [61]. Around 10 mg of leaf samples were extracted with 1 mL of 5% TCA, sonicated three times for 10 min and centrifuged at 12,000× g for 10 min at 4 • C. This assay is based on an enzymatic recycling procedure in which glutathione is sequentially oxidized by 5,5 -dithiobis-2-nitrobenzoic acid and reduced by NADPH in the presence of glutathione reductase. The GSSG was determined after removal of GSH from the sample extract by derivatization with 4-vinylpyridine. Absorbances were measured at 412 nm with same spectrophotometer reported above and converted to concentration with a GSH standard curve. Concentrations of the GSH were calculated as GSH TOT − GSSG.
Photosynthetic leaf pigments were measured according to Lichtenthaler [62]. Around 10 mg of leaf material were extracted with 1 mL of 100% acetone, sonicated three times for 10 min and centrifuged at 13,000× g for 10 min at 4 • C. The absorbances of sample supernatants were determined at 470, 645, and 662 nm, using the same spectrophotometer reported above. The amounts of chlorophyll a (Chl a) and b (Chl b) and Car were calculated as reported by Lichtentaler and Buschmann [63]. Total chlorophyll (Chl TOT ) was calculated as Chl a + Chl b.
Total phenols (Phen) were determined according to Ainsworth and Gillepsie [64]. Concisely, 5 mg of leaf samples were extracted with 1.9 mL of 95% methanol, sonicated three times for 10 min, incubated at room temperature for 48 h in the dark and centrifuged at 13,000× g for 10 min at 4 • C. Then, 100 µL of sample supernatant were collected and mixed with 200 µL 10% Folin-Ciocalteu reagent and 800 µL of 700 mM Na 2 CO 3 . Samples were incubated at room temperature for 2 h and then measured for absorbance at 765 nm with the same spectrophotometer reported above. Absorbances were converted to concentration with a gallic acid standard curve.

PLSR-Model Calibration and Validation
The PLSR [16] models were generated from untransformed reflectance profiles to predict A, g s , C i , E, WUE i , WUE in , k, T l , MDA, ORAC, HORAC, DHA, DHA/ASA TOT , GSH, GSH TOT , Chl a, Chl TOT , Car, and Phen. To avoid potential overfitting, the PLSR models, the numbers of latent variables to use were selected based on reduction of the predicted residual sum of squares (PRESS) statistics [65] using leave-one-out cross-validation. Finally, the selected sets of extracted components were combined into linear models predicting leaf traits based on leaf spectral profiles.
Model performance was evaluated by conducting 500 randomized permutations of the datasets using 80% of the data for calibration and the remaining 20% for validation (internal validation). For each permutation, we calculated statistics to assess model performance when applied to the calibration and the validation data sets: the R 2 , the overall error rate (RMSEC, root mean square error; the RMSE for the validation dataset is known as root mean square error for prediction, RMSEP), the percentage of error over the data range (%RMSE), and bias. These randomized analyses generated a distribution of fit statistics allowing for the assessment of model stability as well as uncertainty in model predictions. The strength contribution of PLSR loadings by individual wavelengths was also determined using the variable important to the projection (VIP) statistics. These statistics highlight the importance of individual wavelengths in explaining the variation in both the response and predictor variables, with larger weightings confer greater value to contribution of individual wavelengths to the predictive model [16,66].
Before developing the final modelling, we tested preliminary models to identify poorly predicted outliers. Prediction residuals were used to identify potential outliers, following Couture et al. [18]. Spectral data of outliers were further examined for errors, which included elevated reflectance in the VIS region, misaligned detector splicing producing jumps, or concave spectral shape at the red-edge peak, all a result of the leaf clip is not being fully closed during the reference or target measurements. Standard data of outliers were also examined for extremes in the data distribution. Outliers removed accounted for approximately 10% and 20% of the initial data for gas exchange (less steps and complexity in determinations) and biochemical traits (more steps and complexity in determinations), respectively. The modelling approach and data analyses were performed using the "pls" package in R (www.r-project.org).
Best practices suggest to further external validate the developed PLSR-models (i.e., test their accuracy in prediction on other independent samples not used in model development). This operation was not performed in the current study because of the limitation in sample numerousness. Thus, although outputs for validation are usually in agreement with external validation (e.g., [42,48]), we encourage an external validation before using coefficients from PLSR-models here reported.
Other leaf traits were calculated by applying the coefficients of the best performing PLSR-models. Both spectral indices and spectra-derived leaf traits were determined from spectra averaged per plant, using plants only dedicated to spectral measurements (as stated above).

Analyses of Spectral Signatures
The influence of O 3 , time and their interactions on the reflectance profiles of sage (averaged per plant) was determined by permutational analysis of variance (PERMANOVA) [67], employing Euclidian measurements of dissimilarity and 10,000 permutations. The same approach was used to determine the influence of O 3 at each time of analysis (i.e., 0, 1, 2, 5, 8, 24 h FBE). Spectral responses at these times were visualized using Principal Coordinates Analysis (PCoA) on the same spectral data utilized for PERMANOVA, using the "vegan" package in R (www.r-project.org; [68]). Principal Coordinates Analysis uses a distance of uncorrelated variables, or principal coordinates, reducing the dimensionality of the data.
Partial least squares discriminant analysis (PLS-DA) [69] was additionally used to determine the ability of spectral data to classify experimental groups that showed statistical significance by PERMANOVA. Similar to PERMANOVA, this approach was firstly used on the whole dataset and then on subsets related to separated times of analysis. The analyses were applied 500 times by splitting observations into different groups of calibration (training) and validation (testing) sets. We used the number of correct classifications both in the calibration and the validation sets across the 500 iterations to evaluate the accuracy of the tested model. The calibration:validation data ratio and the number of components called to obtain the models that would give the best fit to the data were determined by iteratively running the PLS-DA models with different calibration:validation data ratios (i.e., 50:50, 70:30, 80:20) and numbers of components and was based on the highest kappa values returned for the validation models. The PLS-DA was performed using the "caret" and "vegan" packages in R (www.r-project.org; [68,70]).

Other Statistical Analyses
Normal distribution of leaf traits derived from spectra (by both spectral indices and PLSR-models) was preliminary analyzed following the Shapiro-Wilk test. The effects of O 3 (between factor), time (within factor), and their interaction on these leaf traits were tested using a one-way repeated measures analysis of variance (ANOVA). Tukey's test was used as the post-hoc test. Effects with p ≤ 0.05 were considered statistically significant. These univariate statistical analyses were performed in JMP 13.2.0 (SAS Institute Inc. Cary, NC, USA).