Early Detection of Plant Physiological Responses to Different Levels of Water Stress Using Reflectance Spectroscopy

Abstract

Early detection of water stress is critical for precision farming for improving crop productivity and fruit quality. To investigate varying rootstock and irrigation interactions in an open agricultural ecosystem, different irrigation treatments were implemented in a vineyard experimental site either: (i) nonirrigated (NIR); (ii) with full replacement of evapotranspiration (FIR); or (iii) intermediate irrigation (INT, 50% replacement of evapotranspiration). In the summers 2014 and 2015, we collected leaf reflectance factor spectra of the vineyard using field spectroscopy along with grapevine physiological parameters. To comprehensively analyze the field-collected hyperspectral data, various band combinations were used to calculate the normalized difference spectral index (NDSI) along with 26 various indices from the literature. Then, the relationship between the indices and plant physiological parameters were examined and the strongest relationships were determined. We found that newly-identified NDSIs always performed better than the indices from the literature, and stomatal conductance (G_s) was the plant physiological parameter that showed the highest correlation with NDSI(R₆₀₃,R₅₅₈) calculated using leaf reflectance factor spectra (R² = 0.720). Additionally, the best NDSI(R₆₈₅,R₄₁₅) for non-photochemical quenching (NPQ) was determined (R² = 0.681). G_s resulted in being a proxy of water stress. Therefore, the partial least squares regression (PLSR) method was utilized to develop a predictive model for G_s. Our results showed that the PLSR model was inferior to the NDSI in G_s estimation (R² = 0.680). The variable importance in the projection (VIP) was then employed to investigate the most important wavelengths that were most effective in determining G_s. The VIP analysis confirmed that the yellow band improves the prediction ability of hyperspectral reflectance factor data in G_s estimation. The findings of this study demonstrate the potential of hyperspectral spectroscopy data in motoring plant stress response.

1. Introduction

One of the significant impacts of climate change due to rising temperatures is alteration in water supply and demand for the plant communities, which are a critical source of food [1]. It has been long recognized that the lack of water availability has an adverse effect on the plant growth process, thus decreasing crop productivity and yield [2]. Therefore, studies focusing on early detection and warning of plant water stress are of considerable significance.

Grapevine (Vitis spp.) is considered the most commercially-important berry crop in the world [3]. North American Vitis species play a vital role in the global grape industry by imparting important pest resistance through their use as rootstocks (the below-ground part makes up the lower stem and roots), their contributions to hybrid scions (the above-ground part of the plant; the scion produces the stem, leaves, flowers and berries) or through planting as an ungrafted plant. Today, grape growing is emerging as a more important part of rural agriculture in the Midwestern United States. For instance, in Missouri, the economic impact of grapes and wine has grown to be worth $1.6 billion with a 16% annual growth rate [4]. However, abiotic and biotic stresses are restrictive factors impacting cultivation of even the grafted V. vinifera ssp. Vinifera and hybrid scions derived from crosses between V. vinifera ssp. vinifera and one of the native North American Vitis species are dominant cultivated grapevines in the Midwest. Moreover, fruit quality and yield will be significantly affected by prolonged drought, which is projected in the Midwestern United States [5].

The severity and duration of a plant’s dehydration determines the impact of water stress on the photosynthetic performance of the plant and whether or not a plant can recover from stress damage when irrigated. In the short term, water stress causes stomata closure, which in turn not only leads to leaf temperature increase by reducing the transpiration rate, but also creates a reduction in CO₂ concentration [6,7]. If there is an imbalance between the absorbed light energy and the energy requirement for carbon fixation, promoted by stomata closure, this may cause over-excitations and subsequent the photodamage to photosystem II (PSII) reaction centers [8]. As a result, the maximum quantum efficiency of PSII will decline in response to the onset of water stress [9]. When plant stress severity or duration exceeds a critical threshold, plants manifest biochemical and morphological symptoms as adaption strategies. Some adaptation strategies include reduction in chlorophyll content, leaf area, premature leaf senescence and stunted growth [10,11]. In general, if the stressor is removed before the damage is visible to the naked eye, plants revive and develop a new physiological standard [12].

Compared to traditional field measurements, remote sensing can provide timely and reliable information about the current plant physiology in a cost-effective and timely manner [13]. In particular, optical remote sensing exploits reflected radiation in the visible (VIS, 400–700 nm), the near-infrared (NIR, 700–1200 nm) and short-wave-infrared (SWIR, 1300–2500 nm) regions of electromagnetic spectrum recorded with use of ground, air- and space-borne sensors. Generally, these regions are strongly correlated with leaf pigment concentration, cell structure and water content, respectively [14,15,16,17,18]. The reflectance data recorded with remote sensors are often used to calculate vegetation indices (VIs), which are a mathematical combination of several bands within the visible and NIR spectral regions. VIs have been proven to be a very simple yet effective approach in estimating biophysical variables [19,20,21]. Narrow band hyperspectral remote sensors, unlike broad band multispectral remote sensors, record reflected radiant energy from the objects at a high number of wavebands, and VIs calculated from hyperspectral data provided new insights into the early detection of plant stress [22,23,24,25]. The photochemical reflectance index (PRI), originally designed for epoxidation state of the xanthophyll cycle pigments, is an example promising narrow band index that can detect non-photochemical quenching, one of the two photo-protective mechanisms [24]. Sun-induced fluorescence (SIF) emission is another alternative photo-protective mechanism that minimizes the damages to the photosystem by re-emitting the excess absorbed light energy as fluorescence [26]. In recent years hyperspectral sensors have demonstrated the feasibility of SIF retrieval at multiple scales and confirmed further possibilities to detect plant stress before it is visible to the naked eye [27,28,29,30,31,32].

Recent studies demonstrated that normalized difference spectral indices (NDSI) using all possible combination of two bands outperformed previously-published indices for prediction of biochemical, biophysical and structural plant parameters [33,34,35,36,37,38,39,40]. Additionally, two-dimensional visualization of coefficient of determination (R²) between NDSIs and biophysical variables provide a clear overview of effective wavebands and spectral regions for determining optimal normalized indices to predict various parameters under study [35]. On the other hand, partial least squares regression (PLSR) has been frequently used for band selection analysis at the leaf and canopy levels [41,42]. PLSR reduces collinearity that commonly exists in hyperspectral data by producing non-correlated, statistically independent latent factors that are a linear combination of the original spectral bands. Both the spectral information and the dependent variables are given equal consideration when latent vectors are regressed against the dependent variable using a cross-validated linear model. The advantages of PLSR approach include: (1) it utilizes the continuous spectrum rather than in a band-by-band type of analysis [43,44]; and (2) statistical over-fitting can be minimized by defining the number of orthogonal latent vectors that produces the smallest root mean square error in cross-validation [45].

The primary objectives of this study were to: (1) investigate the potential of field spectroscopy for characterizing the physiological status of grapevines exposed to different levels of water stress based on in situ measurements; and (2) identify the most effective indices and predictive models for early detection of plant response to water stress using the NDSI and PLSR approaches.

2. Materials and Methods

2.1. Study Site

Field data collection was carried out in a vineyard situated in Mount Vernon, MO, USA (37°4′27.17′′N, 93°52′46.70′′W, altitude 376 m) during the growing seasons of 2014 and 2015 (Figure 1). The region has a continental climate with an average annual temperature of 15.6 °C and mean annual rainfall of 1066.8 mm. The vineyard was established on 25 June 2008 to investigate varying rootstocks and irrigation interactions in tandem. At establishment, six irrigation zones were installed allowing for randomization of blocks for both three different irrigation regimes and four different rootstocks. Chambourcin vines, either own-rooted or grafted onto 1103 Paulsen, 3309 Couderc and Selectioin Oppenheim 4 (SO4) were planted with varying irrigation patterns of either: (i) nonirrigated (NIR); (ii) full replacement of evapotranspiration (FIR); or (iii) intermediate irrigation (INT, 50% replacement of potential evapotranspiration). The vine density was 504 vines ha⁻¹ with 3 m × 3 m row spacing, including 25 rows and 1034 vines in total. The soil texture was a combination of sandy loam, silt loam and loam, with an average pH of 6. Each vine row was oriented in the east to west direction, and the vines were trained with a high wire cordon trellis and spur pruned. To avoid soil erosion, grass was sown between the rows with a weed free strip maintained just below the vines. Vine rows were numbered from the north and plants were numbered from the west. Originally, the entire field received irrigation at a rate necessary to replace evapotranspiration (ET) until the 2014 season. In the 2014 season, the treatments were initiated, and data collection on the vine water status and fruit quality began.

This study’s field measurements involved 9 total rows, Numbers 8–16, where the different irrigation treatments were applied. Each row consisted of 8 plots and 4 vines in each plot with the same rootstocks. For measurement purposes, the third vine from the east side of the plot was chosen. FIR and INT rows were irrigated using a drip irrigation system. To maintain the different treatments during the growing season, both timing and amount of water were determined based on ET calculated from the weather data obtained from a weather station adjacent to the vineyard. The start dates of the treatments were 14 July 2014 (Day of Year: DOY 195) and 1 September 2015 (DOY 244). The end dates were 3 October 2014 (DOY 276) and 9 October 2015 (DOY 282), respectively. Daily average precipitation, maximum and minimum air temperature and precipitation are presented in Figure 2 together with field data collection dates for two study years.

2.2. Field Data Collection

In the years 2014 and 2015, we collected plant physiological and nutritional variables of the vineyard along with leaf spectra using field spectroscopy. The data collection dates were during early fruit set stage (18 June, DOY 169) and late veraison stage (19 August, DOY 231) in 2014 and during berry touch (10 July, DOY 191) and fruit ripening stage (21 September, DOY 264) in 2015. On these measurement days, no clear damages caused by water stress were identified on grapevine leaves.

2.2.1. Plant Physiological Measurements

Midday leaf physiological status was determined using an LI-6400XT Portable Photosynthesis system coupled with a pulse amplitude-modulated (PAM) leaf chamber fluorometer at incident photosynthetic photon flux density (PPFD) level of 1000 μmol m⁻² s⁻¹ generated by a red LED array, with an additional 10% blue light to maximize stomatal opening (Li-Cor, Inc., Lincoln, NE, USA). Measurements were made on a fully-expanded sunlit leaf on the south-facing side of the vine in each plot for gas exchange and fluorescence variables: stomatal conductance (G_s), photosynthetic CO₂ assimilation rate (A_i), chlorophyll fluorescence, electron transport rate (ETR), photochemical (qP) and non-photochemical quenching (NPQ) during midday (1000–1500 h) under full sun conditions. Inside the leaf chamber, the CO₂ concentration was set at 400 μmol CO₂ mol⁻¹ air in the cuvette, and the relative humidity of the incoming air ranged between 40 and 60%; temperature and water vapor pressure deficit (VPD) were not controlled. The chlorophyll fluorescence (ΔF/F_m′, also known as instantaneous photochemical efficiency of PSII) was calculated as (F_m′-F_s)/F_m′, where F_m′ is the fluorescence yield and F_s the steady-state fluorescence of the light-acclimated leaf [46].

2.2.2. Hyperspectral Reflectance Measurements

Reflectance factor data of vine leaves, specifically hemispheric conical reflectance factor (HCRF; [47]), were obtained using a Spectral Evolution portable spectroradiometer PSR-3500 (Spectral Evolution, Inc., Lawrence, MA, USA). The spectroradiometer records spectral information of a target in full wavelength range (350–2500 nm) with a resolution of 3.5 nm in the 350–1000 nm range, 10 nm in the 1000–1900 nm range and 7 nm in the 1900–2500 nm range. The spectroradiometer was equipped with a specifically developed leaf clip for the bifurcated fiber-optic connected to both the device and a 5-watt tungsten halogen lamp light source. With the leaf clip, we collected the reflectance factor of the same leaf used for photosynthetic measurements. This combination is useful for acquiring the leaf hemispherical reflectance factor with a spectrally black background after a spectrally white background on the opposite side of the clip was measured for the reference spectrum. The black background measurements were taken three times at different points avoiding leaf veins, and the averaged leaf reflectance factor was used for further analysis.

We used Getac^® PS336 PDA preloaded with DARWin software (Compact V.1.2.4903, Spectral Evolution, Inc., Lawrence, MA, USA) to manipulate the spectroradiometer and collect data efficiently in the field. The spectroradiometer was configured to average 40 spectra automatically per sampling, and the raw spectra bandwidth was interpolated to 1 nm. This resulted in 2151 individual spectral bands.

In addition to plant physiological and reflectance measurements, manual harvesting was carried out on 5 October 2014 (DOY 278) and 11 October 2015 (DOY 284), respectively, and the berry weight measurements were determined on-site for each individual vine within the plots.

2.3. Methods

2.3.1. Selection of Spectral Indices and Sensitive Bands

One of the approaches to explore the significant relationships between plant physiological parameters and hyperspectral data is conducting a comparative analysis of simple normalized difference spectral indices (NDSI) calculated from narrow band reflectance factor spectra. We identify the wavelengths or normalized indices (Equation (1)) that are capable of estimating physiological parameters such as G_s.

We applied NDSI to identify optimal wavelengths and/or indices. The NDSI is defined as:

$$NDSI\left(i,j\right)=\frac{R_i-R_j}{R_i+R_j},$$

(1)

where R is the reflectance factor value, and the subscripts are wavelengths in nanometer (nm).

NDSIs were calculated for the measured leaf hyperspectral reflectance factor spectra using all possible combinations of available bands (i and j nm) in the full spectral region (350–2500 nm), excluding the 1350–1450 nm, 1800–2000 nm and 2300–2500 nm regions due to strong atmospheric H₂O and CO₂ absorption. Here, we examine the linear relationship between in situ physiological parameters and NDSIs and 2-dimensional maps of coefficient of determination (R²). The R² maps allow the evaluation of the different band combinations and the selection of a sensitive NDSI for each physiological parameter under study [35,39,40]. The performance of identified NDSIs was then compared with previously-published indices from the relevant literature as shown in

Table 1. Spectral indices used in this study.

Reflectance Index	Acronym	Equation	References
Leaf pigment
Anthocyanin (Gitelson)	AntGitelson	AntGitelson = (1/R550 − 1/R700) × R780	[48]
Carotenoid Reflectance Index	CRI1	CRI1 = 1/R510 − 1/R550	[49]
Carotenoid Reflectance Index	CRI2	CRI2 = 1/R510 − 1/R700	[49]
Chlorophyll Index	CI	CI = (R750 − R705)/(R750 + R705)	[50]
Optimized Soil-Adjusted Vegetation Index	OSAVI	OSAVI = (1 + 0.16) × (R800 − R670)/(R800 + R670 + 0.16)	[51]
Red Green Index	RGI	RGI = R690/R550	[52]
Structure Intensive Pigment Index	SIPI	SIPI = (R800 − R450)/(R800 + R650)	[53]
Transformed Chlorophyll Absorption in Reflectance Index	TCARI	TCARI = 3 × ((R700 − R670) − 0.2 × (R700 − R550) × (R700/R670))	[54]
	TCARI/OSAVI	TCARI/OSAVI	[54]
Normalized Pigment Chlorophyll Index	NPCI	NPCI = (R680 − R430)/(R680 + R430)	[55]
Greenness
Enhanced Vegetation Index EVI	EVI	(2.5(R782 − R 675)/(R782 + 6 × R675 − 7.5 × R445 + 1))	[56]
Normalized Difference Vegetation Index	NDVI	NDVI = (R800 − R670)/(R800 + R670)	[57]
Greenness Index	GI	GI = R554/R677	[52]
Green NDVI	GNDVI	GNDVI = (R750 − R540 + R570)/(R750 + R540 − R570)	[58]
Red Edge Inflection Point	REIP	REIP = 700 + 40 × {[(R670 + R780)/2 − R700]/(R740 − R700)}	[59]
Simple Ratio	SR	SR = R900/R680	[57]
Triangular Vegetation Index	TVI	TVI = 0.5 × (120 × (R750 − R550) − 200 × (R670 − R550))	[60]
Stress
Fluorescence Ratio Index 1	FRI1	FRI1 = R690/R600	[61]
Fluorescence Ratio Indices 2	FRI2	FRI2 = R740/R800	[61]
Modified Red Edge Simple Ratio Index	mRESR	mRESR = (R750 − R445)/(R705 + R445)	[20]
Normalized Phaeophytinization Index	NPQI	NPQI = (R415 − R435)/(R415 + R435)	[62]
Photochemical Reflectance Index	PRI	PRI = (R531 − R570)/(R531 + R570)	[24]
Plant Senescence Reflectance Index	PSRI	PSRI = (R680−R500)/R750	[63]
Red-Edge Vegetation Stress Index	RVSI	0.5(R722 + R763) − R733	[64]
Simple Ratio Pigment Index	SRPI	SRPI = R430/R680	[65]
Water
Water Index	WI	WI = R900/R970	[66]

. The most effective NDSIs and the indices from the literature were determined using R² and root mean square error of calibration (RMSE_cal) on the randomly-selected calibration dataset (80% of samples, n = 169). The predictive ability the indices was evaluated using R² and RMSE_val on the independent validation dataset (20% of samples, n = 42).

2.3.2. Multivariate Method

The partial least squares regression (PLSR) is a multivariate regression method that specifies a linear relationship between a set of dependent (response) variables, Y, and a set of predictor variables, X [67]. It is a powerful tool specifically designed to deal with the data consisting of many independent variables and is used to reduce collinearity within the data to non-correlated latent variables or factors [68,69,70]. To select the optimum number of factors and avoid overfitting, we calibrated the model by an iterative leave-one-out cross-validation criterion called the minimum predicted residual sum of squares (PRESS) RMSE_cal. RMSE_cal is minimized by iteratively leaving one sample out of the calibration dataset and calibrating the model from the remaining dataset. To further evaluate the significance of each wavelength for model prediction, variable importance of projection (VIP) values were calculated [71]. The higher the VIP value of a wavelength is, the greater its contribution to the model becomes. Thus, the wavelengths with VIP-values greater than 1 are the most influential predictors in a model. The predictive ability of the best selected PLSR model was assessed using the R² and RMSE_val on the independent validation dataset.

2.3.3. Sensitivity Analysis for Early Stress Detection

As mentioned in the Introduction section, water stress symptoms in the early stage are not visible. To avoid yield loss, the stress factor must be removed before irreversible damage occurs. Therefore, it is important to test the capability of newly-identified NDSI and the PLSR model for early stress detection. To evaluate the sensitivity of the best selected NDSI model and the PLSR model to the induced water stress based on field measured physiological data, a one-way analysis of variance (ANOVA) followed by an honest significant difference (HSD) Tukey test (α = 0.05) was performed to determine significant differences between predicted plant physiological parameters corresponding to the observed different stress levels. Since there were three irrigation treatments (NIR, INT and FIR) intended to induce different levels of water stress within the vineyard, we represented physiological status of grapevine treated with INT irrigation as the early stage of water stress.

2.3.4. Statistical Analyses

One-way analysis of variance (ANOVA) was used to investigate the effect of irrigation treatments on the field measured plant physiological parameters. Significant differences between treatments were assessed with HSD Tukey tests (α = 0.05). Pearson correlation coefficients were used to explore the significant relationships between plant physiological parameters. All data analyses were performed using IBM SPSS Version 24.0 (SPSS Inc., Chicago, IL, USA).

3. Results

3.1. Vineyard Weather Condition and Physiological Responses of Grapevines to Induced Water Stress

The average air temperature was higher by 1.7 °C in 2014 for the period covering the stress initiation through the second field measurement date (DOY 195–DOY 231) compared to the same period in 2015 (DOY 244–DOY 264). There was only one >20-mm rain event that occurred during this deficit irrigation period in 2014. In 2015, however, there were two >20-mm events. In both years, there were nine no-rain days prior to the second field measurement date.

In 2014, the mean values of leaf-level physiological parameters were generally lower on 19 August (DOY 231) than the mean values on 18 June (DOY 169), although there were rainfall events during the deficit irrigation period (

Table 2. Grapevine physiological parameters as a function of the irrigation treatments in 2014. Values indicate the mean and pooled RMSE of the ANOVA test. G_s: stomatal conductance (mol H₂O m⁻² s⁻¹); A_i: photosynthetic CO₂ assimilation rate (µmol CO₂ m⁻² s⁻¹); F_s: steady-state fluorescence; F_m': maximum fluorescence; ΔF/F_m': fluorescence yield; ETR: electron transport rate (µmol m⁻² s⁻¹); qP: photochemical; NPQ: non-photochemical quenching; RMSE is pooled root mean square error from ANOVA. Different letters in the column indicate significant differences among treatments according to Tukey’s test (p < 0.05). * p < 0.05, ** p < 0.01, *** p < 0.001; ns., not significant.

18 June 2014	Gs	Ai	Fs	Fm'	ΔF/Fm'	ETR	qP	NPQ
FIR (n = 12)	0.37	20.3	811	1221	0.33	145	0.64	5.19
INT (n = 12)	0.31	19.7	800	1183	0.32	140	0.62	5.53
NIR (n = 8)	0.39	20.1	784	1211	0.35	153	0.67	5.18
RMSE	0.08	3.38	84	145	0.04	17	0.04	0.90
Significance level	ns.	ns.	ns.	ns.	ns.	ns.	ns.	ns.
19 August 2014
FIR (n = 12)	0.17a	21.1a	827a	1213a	0.32a	139a	0.63	5.00a
INT (n = 12)	0.06b	16.0b	768b	1067b	0.28b	121b	0.59	5.89b
NIR (n = 12)	0.04b	15.4b	721b	983b	0.26b	115b	0.57	6.46b
RMSE	0.04	3.4	54	94	0.04	20	0.07	0.59
Significance level	***	**	***	***	*	*	ns	***
10 July 2015
FIR (n = 24)	0.30	7.88	694	962	0.27	117	0.57	1.88
INT (n = 24)	0.33	7.35	711	1006	0.28	124	0.59	1.93
NIR (n = 24)	0.33	9.45	694	967	0.27	117	0.57	1.87
RMSE	0.08	3.85	83	180	0.07	32	0.10	0.2
Significance level	ns.	ns.	ns.	ns.	ns.	ns.	ns.	ns.
21 September 2015
FIR (n = 24)	0.19	17.4	955	1359	0.29	128	0.56	2.10
INT (n = 24)	0.18	17.0	1008	1408	0.29	125	0.54	2.12
NIR (n = 24)	0.18	17.3	921	1369	0.32	140	0.61	2.13
RMSE	0.03	3.8	172	219	0.07	32	0.11	0.21
Significance level	ns.	ns.	ns.	ns.	ns.	ns.	ns.	ns.

). On 19 August, as expected, INT and NIR treatments showed significantly lower photosynthetic activities than the FIR treatment except for the qP, but the photosynthetic differences between INT and NIR treatments were not significant. Regarding the 2015 data, G_s decreased for all treatments on 21 September (DOY 264) compared to G_s on 10 July (DOY 191), yet the differences between treatments were not significant. Despite the irrigation treatment and nine no rain days prior to the measurement date, INT and NIR treatments closely followed the photosynthesis levels of the FIR treatment on 21 September (Table 2). The relationship between A_i and G_s was consistent in both study years, and the relationship was stronger during the deficit irrigation period (r = 0.82, p < 0.01 for 19 August 2014; r = 0.76, p < 0.01 for 21 September 2015). Physiological parameters that showed significant correlation with A_i were F_s, F_m′, ΔF/F_m', ETR and NPQ, but these relationships were not consistent in both years (

Table 3. Pearson correlation coefficients (r) between grapevine physiological parameters measured in 2014 and 2015. G_s: stomatal conductance (mol H₂O m⁻² s⁻¹); A_i: photosynthetic CO₂ assimilation rate (µmol CO₂ m⁻² s⁻¹); F_s: steady-state fluorescence; F_m': maximum fluorescence; ΔF/F_m': fluorescence yield; ETR: electron transport rate (µmol m⁻² s⁻¹); qP: photochemical; NPQ: non-photochemical quenching; * p < 0.05, ** p < 0.01, *** p < 0.001.

	Gs	Ai	Fs	Fm'	ΔF/Fm'	ETR	qP	NPQ
18 June 2014
Gs	1.00
Ai	0.75 **	1.00
Fs	0.05	0.21	1.00
Fm'	0.27	0.37 **	0.85 **	1.00
ΔF/Fm'	0.39 *	0.34	0.05	0.56 **	1.00
ETR	0.39 *	0.34	0.05	0.56 **	0.99 **	1.00
qP	0.26	0.06	−0.32	0.16 *	0.85 **	0.85 **	1.00
NPQ	−0.23	−0.31	−0.82**	−0.97 **	−0.59 **	−0.59 **	−0.26	1.00
19 August 2014
Gs	1.00
Ai	0.82 **	1.00
Fs	0.47 **	0.29	1.00
Fm'	0.50 **	0.34 *	0.82 **	1.00
ΔF/Fm'	0.26	0.22	0.15	0.68 **	1.00
ETR	0.14	0.22	0.15	0.69 **	0.99 **	1.00
qP	0.13	0.14	−0.05	0.51 **	0.96 **	0.96 **	1.00
NPQ	−0.50	−0.31	−0.81 **	−0.99 **	−0.69 **	−0.69 **	−0.53 **	1.00
10 July 2015
Gs	1.00
Ai	0.31 **	1.00
Fs	−0.01	−0.08	1.00
Fm'	−0.10	0.00	0.87 **	1.00
ΔF/Fm'	−0.18	−0.10	0.40 **	0.84 *	1.00
ETR	0.19	−0.10	0.48 **	0.84 **	0.99 **	1.00
qP	−0.22	−0.14	0.24 *	0.66 **	0.95 **	0.95 **	1.00
NPQ	−0.12	−0.04	0.74 **	0.96 **	0.91 **	0.91 **	0.75 **	1.00
21 September 2015
Gs	1.00
Ai	0.76 **	1.00
Fs	0.11	−0.10	1.00
Fm'	0.29 *	0.10	0.79 **	1.00
ΔF/Fm'	0.26 *	0.31 **	−0.39 **	0.25 *	1.00
ETR	0.26 *	0.31 **	−0.39 **	0.25 **	1.00 **	1.00
qP	0.18	0.29 *	−0.65 **	−0.06 *	0.93 **	0.93 **	1.00
NPQ	0.31 **	0.20	0.37 **	0.82 **	0.63 **	0.63 **	0.31 **	1.00

). Failing to induce different levels of water stress in 2015 did not affect the results of this work, as the purpose of the study was to find the best indicator of water stress, not differentiating levels of water stress. In particular, we selected G_s as the water stress indicator because it is considered an important parameter to assess water stress by Flexas et al. (2002a) [72], and stomata respond to water stress before there is a detectable change in the leaf water potential and/or leaf water content [73,74].

Across the seasons, berry yield was higher in 2014 (14.51 kg/vine) than in 2015 (14.21 kg/vine), but the difference was not statistically significant (p > 0.05). Irrigation only affected the yield in 2014 and there was a significant (p < 0.05) difference between the FIR (13.82 kg/vine) and INT (15.70 kg/vine) treatments.

3.2. Complete-Combination Indices Analysis of the Hyperspectral Reflectance Factor Data

Figure 3a shows one of the two-dimensional maps of the coefficient of determination (R²) calculated using a full spectral range (350–2500 nm) and G_s. The NDSIs with the highest R² and the lowest p-values (R² > 0.7 and p < 0.001) were found in the VNIR (visible and near infrared) spectral region. In most cases, the SWIR region (1300–2500 nm) correlated relatively well with G_s when only combined with the VNIR region (R² ≤ 0.69). The similar results were obtained for other physiological parameters, yet with very weak correlations (data not shown). The VNIR bands have a higher signal to noise ratio than SWIR bands. Besides, the SWIR bands of imaging spectroscopies tend to have lower spatial resolution than the VNIR bands due to the physical limitations. Using indices formed with the combination of VNIR and SWIR would be another source of error caused by resampling to the same spatial resolution [75]. Moreover, longer wavelengths tend to be affected by atmospheric water vapor, as the water absorption coefficient becomes higher, and it is hard to completely remove the effects of water vapor, even with a good atmospheric correction [76]. Therefore, we focused on the 400–1100 nm spectral region corresponding to VNIR bands for the rest of this study. In

Table 4. Maximum values of coefficients of determination (R²) and root mean square error of calibration (RMSE_cal) between the grapevine physiological parameters and selected NDSIs.

Spectral Indices	Gs	Ai	Fs	Fm'	NPQ
	R2 (RMSEcal)	R2 (RMSEcal)	R2 (RMSEcal)	R2 (RMSEcal)	R2 (RMSEcal)
NDSI(603,558)	0.720 (0.063)	-	-	-	-
NDSI(728,525)	0.711 (0.066)	-	-	-	-
NDSI(830,525)	0.703 (0.068)	-	-	-	-
NDSI(1000,525)	0.694 (0.068)	-	-	-	-
NDSI(715,620)	0.715 (0.066)	-	-	-	-
NDSI(818,620)	0.707 (0.067)	-	-	-	-
NDSI(1000,620)	0.709 (0.067)	-	-	-	-
NDSI(726,630)	0.714 (0.068)	-	-	-	-
NDSI(778,635)	0.707 (0.069)	-	-	-	-
NDSI(1000,635)	0.707 (0.070)	-	-	-	-
NDSI(705,535)	-	0.256 (5.268)	-	-	-
	-
NDSI(704,540)	-	-	0.275 (140.281)	-	-
NDSI(704,540)	-	-	-	0.284 (208.748)	-
NDSI(685,415)	-	-	-	-	0.681 (0.992)

Notes: The NDSIs presented here were selected according to the R² maps, a part of which is shown in Figure 3b–e. For G_s and NPQ, NDSIs with R² values higher than 0.7 and 0.6 are presented in this table, respectively. All relationships are significant at p < 0.001 level.

, the strength of the relationship between selected NDSIs and photosynthetic parameters is presented. Figure 3b–e shows a representative part of the NDSI analysis that was conducted for all leaf physiological parameters under study. The coefficient of determination (R²) for the relationship between the leaf reflectance factor spectra and G_s is presented in Figure 3b. The most highly correlated was NDSI(R₆₀₃,R₅₅₈) (R² = 0.720; RMSE_cal = 0.063) with G_s is in the yellow (570–630 nm) and green (530–580 nm) spectral region. NDSIs between the red (630–680 nm) and red-edge (690–750 nm) region also resulted in significant and high R² (>0.7) values. Additionally, a broad spectral region with the combination of red and red-edge over the NIR region correlated well with G_s. In general, except NDSI(R₆₀₃,R₅₅₈) mentioned previously, G_s was also strongly correlated with other NDSIs (Table 4).

The strength of the relationship between NDSIs was moderate for A_i, F_s and F_m', and negligible for ΔF/F_m', ETR and qP (Table 4). A_i is best correlated with NDSI(R₇₀₅,R₅₃₅) (R² = 0.256; RMSE_cal = 5.268) and the most significant region was narrow between 520 and 550 nm. The distribution of significant spectral regions identified for the F_s and F_m’ are very similar, as shown in Figure 3c,d. In the NDSI map (Figure 3c), F_s has a maximum R² = 0.275 at (R₇₀₄,R₅₄₀). This significant region was narrow (approximately 25 nm) along 700 nm (R_i), but relatively wide over 515–595 nm (R_j). Similarly, for F_m', NDSI(R₇₀₄,R₅₄₀) with the highest R² = 0.284 value was found in the 700–725 nm and 520–590 nm spectral regions. For both parameters, the significant spectral region extended from around 700 nm toward longer wavelengths up to 740 nm, covering the far-red chlorophyll fluorescence emission region.

The spectral region formed with blue (400–450 nm) and red (670–690 nm) regions showed the highest correlation with NPQ (R² > 0.6, Figure 3e). Accordingly, NDSI(R₆₈₁,R₄₁₅) within this region had a maximum R² (0.663). NDSIs between the green (510–530 nm) and blue (400–500 nm) region also resulted in significant and high R² values (>0.4 and <0.6). Furthermore, a broad spectral region with the combination of green and red over the NIR region showed relatively good correlation with NPQ.

Figure 4a,b shows the predictive ability of the models for the best NDSI(R₆₀₃,R₅₅₈) and NDSI(R₆₈₁,R₄₁₅) in G_s and NPQ estimation using the independent validation dataset. The significant spectral regions identified for G_s and NPQ were broad and useful in sensor applications. However, regarding the consistent and significant correlations between G_s and A_i, we only focused on G_s for the rest of the study.

3.3. The Relationship between the Grapevine Water Stress Response and Hyperspectral Reflectance Indices from the Literature

Table 5. Coefficients of determination (R²) and root mean square error of calibration (RMSEcal) between the grapevine physiological parameters and indices from the literature. * p < 0.05, ** p < 0.01, *** p < 0.001.

	Gs	NPQ
	R2 (RMSEcal)	R2 (RMSEcal)
Leaf Pigment
AntGitelson	0.567 (0.100) ***	0.269 (1.509) ***
CRI1	0.663 (0.087) ***	0.332 (1.441) ***
CRI2	0.648 (0.090) ***	0.207 (1.571) ***
CI	0.648 (0.090) ***	0.529 (1.211) ***
OSAVI	0.662 (0.088) ***	0.485 (1.261) ***
RGI	0.592 (0.097) ***	0.364 (1.407) ***
SIPI	0.632 (0.092) ***	0.376 (1.394) ***
TCARI	0.537 (0.103) ***	0.368 (1.403) ***
TCARI/OSAVI	0.587 (0.098) ***	0.433 (1.327) ***
NPCI	0.518 (0.106) ***	0.544 (1.191) ***
Greenness
EVI	0.156 (0.140) **	0.038 (1.731) **
NDVI	0.662 (0.088) ***	0.529 (1.211) ***
GI	0.277 (0.130) ***	0.057 (1.714)**
GNDVI	0.218 (0.135) ***	0.045 (1.724)**
REIP	0.540 (0.104) ***	0.456 (1.301) ***
SR	0.680 (0.086) ***	0.531 (1.196) ***
TVI	0.257 (0.131) ***	0.086 (1.687) ***
Stress
FRI1	-	0.046(1.723)**
FRI2	0.314 (0.126) ***	0.356 (1.416) ***
mRESR	0.656 (0.089) ***	0.383 (1.386) ***
NPQI	0.222 (0.134) ***	0.094 (1.680) ***
PRI	0.022 (0.151)	0.043 (1.725) **
PSRI	0.278 (0.130) ***	-
RVSI	0.263 (0.131)	-
SRPI	0.530 (0.105) ***	0.450 (1.309) ***

summarizes the R² and corresponding RMSE_cal for the previously published indices in estimation of the leaf level photosynthetic parameters. G_s showed moderate to strong correlations with mostly pigment and greenness-based indices. The best correlations were obtained with the carotenoid reflectance Index 1 (CRI1), optimized soil-adjusted vegetation index (OSAVI), NDVI and the simple ratio (SR). The modified red edge simple ratio index (mRESR) was the only index among stress based indices that was well correlated with G_s. The highest and more frequent correlations were observed for NPQ with all of the indices except for WI. The Normalized Pigment Chlorophyll Index (NPCI) appeared to be the best index for estimating the NPQ as it had the highest R² and lowest RMSE_cal. Traditional indices (e.g., NDVI and SR) outperformed the improved once (EVI and triangular vegetation index (TVI)). It is worth noting that several stress based indices (e.g., PRI, plant senescence reflectance index (PSRI) and Red-Edge Vegetation Stress Index (RVSI)) showed very low to non-correlation with NPQ (R² < 0.2). As opposed to G_s and NPQ, the rest of the photosynthetic parameters did not correlate well with all of the published indices (R² < 0.17). Lastly, it is important to note that stress-based indices (e.g., FRI1 and FRI2), designed for chlorophyll fluorescence estimation, showed a very weak correlation with fluorescence-related photosynthetic parameters (data not shown).

3.4. PLSR Analysis

The PLSR analysis of the hyperspectral reflectance factor data and the G_s showed that the highest R² (0.680) and the lowest RMSE_CV (0.065) value were found with three factors (which explains 98.27% of the variance in the predictive variables). This PLSR model with three factors was selected based on the rule that the addition of another factor should reduce the RMSE_CV by more than 2% [77,78]. The influence of each wavelength in the PLSR model is illustrated in Figure 5 with corresponding VIP values. The VIP method revealed the importance of the 400–720 nm region for G_s. Particularly, the local maximum VIP values were found with 522, 604, and 700 nm, while 604 had the highest VIP value of 1.29. Figure 6 shows the predictive ability of the best PSLR model in G_s estimation using the independent validation dataset. Overall, however, the PLSR models did not result in any improvement in terms of variance explained compared to the NDSI band selection methods for G_s estimation.

3.5. Feasibility of Early Stress Detection

Stomatal closure induced by water stress is associated with changes in other photosynthetic parameters. Based on the significant physiological differences and absence of visible stress symptoms between treatments in August 2014, we tested the sensitivity of the best selected NDSI(R₆₀₃,R₅₅₈) and the PLSR model in discriminating the difference observed in G_s. Figure 7 shows the ANOVA performed on the predicted G_s using the leaf reflectance factor data collected on 19 August 2014. NDSI(R₆₀₃,R₅₅₇) (Figure 7a) was able to differentiate NIR and INT treatments from FIR treatment with a significance of p < 0.05 compared to FIR treatment. The PLSR model (Figure 7b) showed similar patterns for the treatments; however, the difference is less significant (p > 0.05) than that of Figure 7a.

4. Discussion

Pure leaf reflectance factor spectra were used in this study to assess the plant physiological response to the induced water stress using data collected on four measurement dates in 2014 and 2015.

Even though there were nine no rain days prior to the second field measurement dates in both years, the induced water stress was obvious only in 2014. This could be attributed to the drier and high temperature during the growing season of 2014 compared to 2015. After stress initiation in 2015, a partial stomatal closure was detected, but other physiological parameters were higher than before stress initiation. This could be explained by nine no rain days prior to the measurement date and lower temperature on that specific measurement date. Furthermore, the reason for this could be optimized G_s by the plants to maximize photosynthesis and minimize water loss [79].

Reduction in photosynthesis is a commonly-observed response of plants to water stress [80]. Depending on the intensity of water stress, both stomatal and non-stomatal limitations are responsible for photosynthetic decline [81,82]. The consistent and significantly positive relationship between G_s and A_i indicates that a main limiting factor for after treatment initiation could be stomatal closure. However, in a study by Zarco-Tejada et al. [83], 2013, the best relationship for A_i was found with F_s consistently for two years. This is because when grapevines are exposed to mild to moderate water stress, stomatal closure acts as a dominate photosynthetic limiting factor until G_s reaches below 0.1 mol H₂O m⁻² s⁻¹, if a further decrease occurs in G_s due to prolonged stress leading to noticeable change in F_s as the non-stomatal factor becomes dominant [84].

Despite the fact that there was no significant difference between treatments in most of the instances, G_s was still found to be the plant physiological parameter that showed the highest correlation with NDSIs calculated using leaf reflectance factor spectra (the highest R² = 0.720). This result is in agreement with the findings of Sellers et al. [85], Myneni et al. [86], Verma et al. [87] and Carter [88], who reported strong linear or non-linear relationships between vegetation indices calculated using VIS and NIR bands, and G_s. G_s is an integrative stress indicator by responding to all external (soil water availability and vapor pressure deficit) and internal (abscisic acid, xylem conductivity, chlorophyll content and leaf water status) influences caused by water stress [89,90,91,92]. This may explain the broad spectral regions with the combination of different wavelengths that correlated well with G_s. NDSI(R₇₂₈,R₅₂₅), NDSI(R₇₁₅,R₆₂₀) and NDSI(R₇₂₆,R₆₃₀) for G_s and NDSI(R₇₀₅,R₅₃₅) for A_i emphasized the importance the red-edge region (695–730 nm). These results agree with the findings of Carter [88], who reported that the red-edge region around 701 nm is the most suitable for G_s and A_i estimation because of its sensitivity to subtle changes in chlorophyll. Evidently, the red-edge region is sensitive to chlorophyll concentration, but this region is also sensitive to changes in cell structure as it is closer to the NIR region. Therefore, the red-edge region is more likely to respond to changes in leaf cell structure separate from pigmentation [93,94]. This could explain why NDSI(603, 558) was superior, where green and yellow bands are less sensitive to changes in leaf structure and water content [95]. Furthermore, both green and yellow bands are known to be sensitive to subtle changes in photosynthetic pigments [58,95,96]. In addition, the yellow band highlights the negative change in slope from the green peak around 550 nm [97]. Figure 8 depicts the effects of irrigation treatments on leaf reflectance factor spectra. Apparently, the effect of treatments on 557 nm is minimum, whereas 603 nm is highly effected by treatments. It is also worth noting that 603 nm is separate from the strong chlorophyll absorption band 680 nm, which is less sensitive to more subtle changes in photosynthesis pigment concentrations [98,99].

A very similar and significant combination of spectral regions (R² is highest at 704 and 540 nm) was found for F_s and F_m', and these results somewhat agree with the findings of Stratoulias et al. [40], who reported the combination of similar spectral regions (including red and far-red fluorescence emission regions) correlated with F_s and F_m' for deep water reed plants with lower chlorophyll concentration, but the shapes of the R² maps from their study showed broad significant spectral regions with diffuse edges compared to the ones in our study (Figure 4b,c). Although plants only emit 2–5% of the absorbed sun light energy as fluorescence, the fluorescence emission spectrum spans a broad spectral region of red and far-red (600–800 nm) with two distinct peaks at around 685–690 nm and 730–740 nm [100,101]. The red part of the fluorescence emission is subject to reabsorption because this part of the fluorescence emission spectrum overlaps with the chlorophyll absorption spectrum, while the far-red part is minimally affected [102,103]. Therefore, red fluorescence emission is low at the sensor especially when chlorophyll is not markedly damaged, and this explains why our identified significant spectral regions were located in far-red region for both F_s and F_m'.

NDSIs calculated with leaf reflectance factor spectra were weakly correlated with ΔF/F_m', ETR and qP. This could be explained by the high level heterogeneity of grapevines photosynthetic activity affected by irrigation treatments, even though the differences in grapevine physiology were not consistently significant. In contrast, there was a strong correlation between NDSIs and NPQ, and the strength of the correlation was comparable to that of G_s, suggesting that grapevines were more subject to photosynthetic down-regulation via NPQ, probably related to the xanthophyll cycle [104,105]. In general, NDSI(R₆₈₅,R₄₁₅) showed the highest correlation for NPQ and selected spectral bands corresponding to the absorption maxima of chlorophyll and carotenoids. The latter compose the xanthophyll cycle to protect the PSII from photodamage by thermal dissipation [104].

G_s response to a subtle leaf internal structure and pigment changes could explain the moderate to low (yet frequent) correlation of greenness and pigment-based indices from the literature with G_s. The NDSIs identified for G_s showed R² higher than 0.7, while the best performed indices from the literature had the highest R² of 0.68. These results are somewhat expected regarding RGI and NDVI because these indices were the specific cases of NDSIs, and their correlations with G_s were already accounted for. Similarly, for NPQ, the NPCI band combination of (R₆₈₀,R₄₃₀) showed the highest R² among the indices from the literature, while the best identified NDSI(R₆₈₅,R₄₁₅) had the R² of 0.681. Moreover, the spectral bands used for these indices can be found within the same spectral regions, e.g., red and blue regions, and therefore, similar outputs are expected. The PRI, an indicator of the epoxidation state of xanthophyll, turns out to be one of the worst performing indices. This is because PRI captures short-term plant photosynthetic performance in response to environmental conditions; for the long term; however, variation in plant structure and foliar pigment concentration act as confounding factors for elucidating the temporal photosynthetic status by using PRI [106,107,108,109,110].

Compared with the NDSI method, the PLSR method did not improve the predictive ability of the model. However, VIP values indicated the importance of the yellow wavelength region. This is consistent to some extent with a recent study conducted by Zengeya et al. [97], who found that the WorldVidw-2 yellow band improved the detection of low nitrogen level within savanna grasses especially at the beginning of the dry season, while the red-edge band was weakly correlated with low nitrogen levels. Up to 75% of the leaf nitrogen is in the chloroplasts that contains chlorophyll [111], and the sensitivity of the yellow band to a low level of nitrogen could be attributed to subtle changes in chlorophyll concentration at the beginning of the dry season, thus G_s. Additionally, Inoue et al. [112] reported the inferior performance of the multivariate models compared to index-based models in canopy chlorophyll content estimation.

Field-measured physiological data in our study presented a relatively wide range of stress by showing significant physiological differences between treatments, and this was especially true for the data collected on the second measurement date (DOY 231) in 2014. Such a difference allowed us to test the capability of effective NDSI(R₆₀₃,R₅₅₈) in early stress detection.

The two-dimensional maps of R² provide a clear overview of effective wavebands and spectral regions for determining optimal normalized indices to predict various parameters under study [35]. Alternatively, multivariate regression methods, machine learning-based support vector machine and non-linear machine learning methods can be employed for leveraging spectral data [15,113,114]. However, these methods are dependent on a large number of high quality training datasets, and their operational applicability is limited [113,115,116].

5. Conclusions

The objectives of this contribution were to (1) investigate the potential of field spectroscopy for characterizing the physiological status of grapevines exposed to different levels of water stress based on in situ measurements and (2) identify the most effective indices and predictive models to evaluate the potential of visible and near-infrared band combinations in the early detection of plant response to water stress. To this end, we utilized leaf reflectance factor spectra and plant physiological parameters measured in a vineyard treated with three different irrigation treatments throughout two consecutive growing seasons. Consistently significant correlation in both experimental years was only found between in situ measured leaf level A_i and G_s, suggesting that stomatal closure acted as a dominate photosynthetic limiting factor. Strong correlations between leaf reflectance factor spectra and plant physiological status demonstrated the potential of ground-based spectroscopy in determining plant photosynthetic performance regardless of different irrigation treatments, growth stages and other environmental factors. In particular, we found that NDSI(R₆₀₃,R₅₅₈) showed the highest R² (0.720) for G_s and NDSI(R₆₈₅,R₄₁₅) for NPQ (R² = 0.681). These indices are spectrally wide and applicable for broad satellite bands. The PLSR model performance was inferior to the NDSI-based method. However, the VIP method provided insights into the importance of the yellow region, confirming that the yellow band improves the prediction ability of hyperspectral reflectance factor data in G_s within the vineyard. Nevertheless, this study only relied on ground-based field spectroscopy. Therefore, the applicability of our findings as an operational tool for predicting plant early stress indicators requires validation using the air- and space-borne data acquired over different plant species exposed to water stress. Scaling up these findings to satellite observations by canopy radiative transfer simulations, unmanned aerial systems-based imaging and current and future satellites systems is the future work of this research.