Combining Genetic Analysis and Multivariate Modeling to Evaluate Spectral Reﬂectance Indices as Indirect Selection Tools in Wheat Breeding under Water Deﬁcit Stress Conditions

: Progress in high-throughput tools has enabled plant breeders to increase the rate of genetic gain through multidimensional assessment of previously intractable traits in a fast and nondestructive manner. This study investigates the potential use of spectral reﬂectance indices (SRIs; 15 vegetation-SRIs; 15 water-SRIs) as alternative selection tools for destructively measured traits in wheat breeding programs. The genetic variability, heritability (h 2 ), genetic gain (GG), and expected genetic advances (GA) of these indices were compared with those of destructively measured traits in 43 F 7-8 recombinant inbred lines (RILs) grown under limited water conditions. The performance of SRIs to estimate the destructively measured traits directly was also evaluated using the partial least squares regression (PLSR) and stepwise multiple linear regression (SMLR) models. Most vegetation-SRIs exhibited high genotypic variation, similar to the measured traits, and phenotypic correlations with these traits, compared with the water-SRIs. Most vegetation-SRIs presented comparable values for h 2 ( > 60%) and GG ( > 20%) as intermediate traits, while about half of water-SRIs exhibited a high h 2 ( > 60%), but low GG ( < 20%). Principle component analysis revealed that most vegetation-SRIs and seven of 15 water-SRIs were grouped together in a positive direction, had a moderate to strong relationship with measured traits, and could identify the drought-tolerant parent Sakha 93 and several RILs. The PLSR model based on all SRIs as a single index showed moderate to high R 2 in calibration (0.53–0.75) and validation (0.46–0.72) datasets, with strong relationships between observed and predicted values of measured traits. The SMLR models identiﬁed four and three SRIs from vegetation-SRIs and water-SRIs, respectively, to explain 63–86% of the total variability in measured traits among genotypes. These results demonstrated that vegetation-SRIs can be used individually or combined with water-SRIs as alternative breeding tools to increase genetic gains and selection accuracy in spring wheat breeding.


Introduction
The combination of water shortage and global climate change will increase the occurrence of drought events in arid and semiarid regions, which will become a critical issue for the sustainability of field crop production in such regions. Improving agronomic practices and cultivation of drought-tolerant cultivars are the most feasible strategies to maintain the sustainability of crop production under limited water supply conditions. These strategies aim to maximize crop production per unit of water applied rather than per unit of area [1,2].
Breeding strategies to develop drought-tolerant genotypes rely on the evaluation of several segregating genotypes or crossing lines using multidimensional assessment of phenotypic differences in several plant traits. In-depth multidimensional descriptions of different plant traits, especially under environmental stress, are urgently needed to improve the understanding of the mechanisms underlying adaptive behavior for stress conditions and to provide plant breeders with comprehensive information about the reliability of various components of plant traits as screening criteria for drought tolerance and the ability to incorporate these traits into breeding programs [3][4][5][6]. The first step to develop drought tolerance in breeding programs is to monitor several morpho-physiological plant traits related to drought tolerance for a large number of crossing lines. However, integrative assessments of several plant traits in a large number of genotypes based on plant sampling and traditional measurement techniques is destructive and inefficient. This emphasizes the urgent need to apply the practical, efficient, and economic high-throughput phenotyping tools, which help breeders to assess key plant traits relevant to stress tolerance for a sufficient number of genotypes in a rapid and non-destructive manner comprehensively. Hyperspectral spectroradiometers contain sensors that can detect small changes in biophysical and biochemical characteristics of the canopy by monitoring changes in canopy spectral signatures at the visible (VIS), near-infrared (NIR), and shortwave-infrared (SWIR) parts of the spectrum. They represent an important phenotyping tool in breeding programs and are used to accelerate the development of genotypes for different environmental stress conditions by performing fast, in-depth phenotypic assessments of various phenological plant traits in multiple genotypes [7][8][9][10].
The concept behind this phenotyping tool is that the behavior of spectral reflectance from the canopy at certain wavelengths in the VIS-SWIR domains of the spectrum is strongly affected by substantial changes that take place in specific plant characteristics. These include changes in photosynthetic efficiency, aboveground biomass, internal leaf structures, leaf pigment content, green leaf area index (LAI), plant water content, and other characteristics related to the growth and physiological states of the canopy. For example, changes in the LAI and fresh and dry aboveground biomass impact the spectral reflectance at 760-1300 nm in spring and winter wheat crop [11][12][13]. Changes in the internal leaf structures due to decreases in leaf relative water content impact the spectral reflectance in the red edge (680-740 nm) and NIR (740-940 nm) regions and increase spectral reflectance in the VIS-NIR regions (400-1300 nm) in wheat and maize crops [14][15][16]. Cochavi et al. [17] reported that changes in photosynthetic pigments and function significantly changed the reflectance of the 680 and 740 nm wavelengths in sunflower crop. The spectral reflectance at several specific wavelengths in the NIR and SWIR regions (i.e., wavelengths centered at approximately 970, 1100, 1200, 1240, 1400, 1450, 1665, 1730, 1950, 2100, 2250, and 2500 nm) was strongly associated with changes in leaf water status in wheat, cotton, and grass crops [18][19][20][21][22][23]. Consequently, the close association between canopy spectral reflectance signatures at specific wavelengths in the VIS, NIR, and SWIR parts of the spectrum and different plant characteristics indicates that the plant spectral reflectance tool can be exploited for integrative assessment of several plant traits of a large number of genotypes in an inexpensive, fast, and nondestructive manner. However, associations between spectral reflectance characteristics and different morpho-physiological properties of plants are usually determined through different specific spectral reflectance indices (SRIs), which are computed from the reflectance at given wavelengths using simple mathematical equations, such as ratios, differences, and sums.
Several SRIs have proven useful in the assessment of different morpho-physiological and agronomic parameters, such as LAI, green area index (GAI), aboveground biomass, photosynthetic capacities, Remote Sens. 2020, 12, 1480 3 of 23 chlorophyll content, grain yield, and different traits related to canopy water content in stressed and healthy plants [24][25][26][27][28][29][30]. Among these, the most commonly used indices to assess plant traits related to growth and yield status are the normalized difference vegetation index (NDVI = [R NIR − R Red ]/[R NIR + R Red ]) and the simple ratio (SR = [R NIR /R Red ]. The simple water index (SWI) and different normalized water indices (NWIs), which focus on the minor water absorption bands within the NIR spectrum (970 nm), have been used to assess parameters related to plant water status at both the leaf and canopy levels. Because the detection of canopy spectral signatures by a spectroradiometer is limited by different soil background conditions, such as soil color and moisture, and is very sensitive to some vegetation biophysical parameters, such as LAI, some SRIs have been developed and modified to reduce the impact on the accuracy of this phenotypic tool. Examples of such indices are the optimized soil-adjusted vegetation index (OSAVI, [31]), enhanced vegetation index (EVI, [32]), modified transformed vegetation index (MTVI, [33]), and modified chlorophyll absorption ratio index (MCARI, [34]). Thus, these indices are able to estimate several traits related to growth and grain yield for different crops and conditions.
The ability of different SRIs to assess various phenological plant traits indicates that these indices may be used as indirect selection tools, replacing destructive and direct traditional phenological selection criteria in breeding trials related to drought studies. In addition, the close relationship between SRIs and phenological plant traits may also identify a true genetic relationship and may therefore help plant breeders reliably evaluate thousands of genotypes in a fast and nondestructive manner after successful incorporation in breeding programs [5,9,24,35,36].
Although it is unclear how different SRIs can estimate the performance of different genotypes based on destructively measured phenological plant traits, particularly under abiotic stress conditions, the existence of a statistical and genetic relationship between the indices and destructively measured plant traits is a proof of concept for the adequacy of these indices as an indirect selection tool. Using SRIs as an indirect selection tool in breeding programs is acceptable if the genetic correlation between these indices and the destructively measured plant traits is very high, and if the heritability of these indices is generally greater than that of the destructively measured traits, especially grain yield [9,24,35,36]. These two features of genetic and other analyses enable plant breeders to determine the efficiency of SRIs as alternative indirect selection tools for destructively measured traits among breeding trials in the future. The different SRIs will be more efficient as indirect selection tools in breeding programs if the genetic analyses for SRIs are conducted on genetic populations closely representing the area where the improvement is being targeted. A comprehensive study is thus urgently required to determine the usefulness of SRIs as indirect selection tools to increase the chance of recognizing wheat genotypes that are well adapted to water shortage under arid conditions. SRIs focus on only two-three wavelengths, reducing the accuracy of indirect estimates of destructively measured traits by a single SRI. Therefore, several studies have focused on multivariate analysis techniques, such as partial least squares regression (PLSR), which deal with a large number of SRIs as a single independent index, to provide more flexibility in the direct estimation of destructively measured traits [28,[37][38][39]. Importantly, because the PLSR tool can create reliable models when the number of SRIs (independent variables) greatly exceeds the number of measured traits (dependent variables), several studies have applied this tool to estimate different measured traits accurately, such as LAI, aboveground biomass, photosynthetic capacities, pigment contents, grain yield, and water relation traits [5,29,40,41].
Therefore, the main objectives of this study were: to (1) evaluate the potential use of different SRIs as alternative indirect breeding tools for breeding purposes under limited water irrigation conditions based on their genetic variability, heritability, genetic gain, and expected genetic advances in advanced breeding wheat lines; (2) evaluate the performance of different PLSR models based on all SRIs as a single index for the accurate estimation of measured traits; and (3) select the most influential SRIs that explain most variability in each destructively measured trait across genotypes via a stepwise multiple linear regression analysis.  Figure S1). The soil texture of the experimental area was sandy loam with a bulk density of 1.50 g cm −3 , along with a field water capacity and wilting point of 0.220 and 0.101 m 3 m −3 , respectively. The soil was characterized as being slightly alkaline (pH 7.9) with low organic matter content (7.8 g kg −1 ). The average monthly temperature during the growing period (December, January, February, March, and April), was 16.4, 14.6, 17.5, 21.5, and 26.9 • C, respectively, and the average monthly precipitation was 12.0, 11.9, 10.0, 11.0, and 3.8 mm, respectively.
Each genotype was planted in 10 row plots at a seeding rate of 150 kg ha −1 with a plot size 1.5 m in width and 4 m in length. The seeds of each genotype were planted in the middle of December in each season. The genotypes were distributed randomly in the plots and arranged in a randomized complete block design with a total of 135 plots. The amount of phosphorus, nitrogen, and potassium fertilizers in each plot were calculated based on the rates of 31, 180, and 120 kg ha −1 of P 2 O 5 , N, and K 2 O, respectively. The entire dose of phosphorus with 20% of nitrogen was applied 2 weeks after sowing, while 40% of nitrogen and half the dose of potassium were applied in the middle of tillering. The remaining dose of nitrogen and potassium was applied at the beginning of the booting growth stage.
Limited water irrigation was applied at a rate of 0.60 of the estimated crop evapotranspiration (ETc). The ETc was estimated based on the crop coefficient (Kc) of spring wheat and the daily reference evapotranspiration (ETo). The Kc values for spring wheat as recommended by FAO-56 were used after adjusting with the actual values of relative humidity and wind speed for the study area. The ETo was calculated by applying the daily meteorological data to the FAO-56 Penman-Monteith equation given by Allen et al. [42]. Based on this calculation, the cumulative amount of irrigation was approximately 3600 m 3 ha −1 , being equivalent to 360 mm ha −1 . About 45 and 40 mm ha −1 of water were applied at sowing and 15 days after sowing, respectively, to ensure full germination and complete establishment of seedlings. Limited water irrigation treatment was started 4 weeks after sowing. The amount of water was applied at the Zadoks growth stages (ZS) of ZS 24, ZS 32, ZS 37, ZS 49, ZS 59, ZS 68, ZS 73, and ZS 83, which referred to the tillering, stem elongation, flag leaf emergence, booting, heading, complete emergence of inflorescence, grain milk, and grain dough of the wheat growth stage, respectively [43]. The irrigation water was applied using a low-pressure surface irrigation system, which had one water emitting tube per plot. To control the amount of irrigation water delivered to each plot, the sub-main hose at each plot was equipped with a manual control valve. The amount of water delivered for each plot was monitored by a flow meter, which connected at the location where the main line encountered the main water source (Supplementary Figure S1).

Measurements of Plant Traits
At the middle anthesis growth stage (ZS 65), thirty plants from the second and ninth row of each plot were randomly selected to determine the green leaf area (GLA), total dry weight (TDW), and leaf area index (LAI). All the green leaves from 30 plants were separated and run through an area meter (LI 3100; LI-COR Inc., Lincoln, NE, USA) to calculate the surface GLA. Thereafter, all parts of the 30 plants (stem, green and complete brown leaves, and spike) were oven-dried at 70 • C for 72 h and then weighed to obtain TDW. Data of the total GLA (A), the number of plants in 1 m (B), the number of plants sampled (C = 30, in this study), and the row spacing (D = 15 cm, in this study) were applied using the following equation to calculate the LAI [27]: Stomatal conductance (Gs) and relative water content (RWC) were measured using the flag leaf of 10 randomly selected plants per plot. The Gs trait was measured at 10:00-13:30 AM using an AP4 diffusion porometer (Delta-T Devices, Cambridge, U.K.). An area of approximately 7-10 cm 2 from the leaves used to measure Gs was excised, and the fresh weight (FW) was immediately recorded. Then, the leaf samples were rehydrated in deionized water at 25 • C until they were fully turgid and rubbed carefully between paper towels to remove external water from the leaf surface before they were weighed to record the turgid weight (TW). Finally, the dry weight (DW) of the leaf samples was recorded after being oven-dried at 70 • C for 72 h. The data of FW, TW, and DW were applied to the following equation to calculate the RWC: After the plants reached maturity, fifty spikes from the third row of each plot were randomly selected and air-dried until the water content reached approximately 15-16% to determine the spike weight (SPW). The five internal rows in each plot, each 3 m in length, were harvested and air-dried for 5 days and then weighed to record the total biological yield (BY). Next, all spikes were separated and threshed, and the final grain yield (GY) was determined. BY and GY values were applied to the following equation to calculate the harvest index (HI) after BY and GY were expressed as ton ha −1 :

Spectral Reflectance Measurements
The spectral reflectance of the whole canopy was measured at ZS 65, simultaneous to the GLA and TDW measurements using a portable high-spectral-resolution spectroradiometer (Fieldspec4, Analytical Spectral Devices Inc., Boulder, CO, USA). This spectroradiometer is able to collect reflectance from the canopy at 350-2500 nm of the light spectrum with finally 1.0 nm band intervals (Supplementary Figure S1). The reflectance of the canopy was determined under cloudless and windless conditions between 10:00 and 14:00. The optical fiber of the sensor was held approximately 80 cm above the plant in the nadir direction with a 25 • field of view to sense a sufficient area of the wheat canopy (~23.4 cm in diameter). Before and as needed during the spectral reflectance measurements, the incoming reflectance from the plant was calibrated using the Spectralon white reference panel (Labsphere, Inc., North Sutton, NH, USA) coated with BaSO 4 . Reflectance was measured from the five internal rows in each plot at four random places, each representing an average of 10 scans, while the average of four sequential measurements and 10 scans per plot was used to estimate 30 spectral reflectance indices (SRIs). These SRIs covered 15 vegetation-SRIs related to the status of growth vigor and photosynthetic efficiency and 15 water-SRIs related to canopy water status, internal leaf structure, and leaf biochemical compounds. The full name, abbreviation, formula, and references of these SRIs are presented in Table 1. Table 1. Full name and abbreviation of the vegetation-and water-spectral reflectance indices (vegetation-SRIs and water-SRIs) used in this study and their algorithms and references.

Data Analysis
Combined analysis of variance for different destructively measured traits and SRIs were performed across seasons based on a randomized complete block design. Combined ANOVA was performed to analyze the genotypic differences across seasons using Bartlett's test and the Shapiro-Wilk test for homogeneity of variances and the normality distribution of the residuals, respectively. Based on this test, the results indicated homogenous variance across seasons for both destructively measured traits and SRIs, and therefore, the data of two seasons were combined and analyzed using the CoStat system for Windows, Version 6.311 (CoHort software, Berkeley, CA 94701). The PROC MIXED model for an across seasons combined analysis of variance was conducted as follows: where y is the phenotypic observation of the genotype for the trait of interest, µ is the overall mean, s i is the effect of the ith season, g j is the effect of the jth genotype, sg ij is the interaction effect between season i and genotype j, r k(i) is the effect of the kth replication nested within season, and ε ijk is the residual error and was assumed to be normally and independently distributed. All effects except for genotypes were considered as random.
Phenotypic correlations between each destructively measured trait and SRI were performed using Pearson's correlation matrix (r).
Genotypic (GCV) and phenotypic coefficients of variability (PCV) were estimated using the following formulae [66]: where genotypic (σ 2 g ) and phenotypic (σ 2 p ) variance were calculated using the formulae given by Fehr [67]. Based on the mean squared values of genotype (GMS) and error (EMS), as well as the number of replicates (r), σ 2 g was calculated as the difference between GMS and EMS divided by r, where σ 2 p was calculated as the summation of σ 2 g and EMS. GCV and PCV values of more than 20% were regarded as high, whereas values less than 10% were considered to be low, and values between 10 and 20% were considered as medium. The high values of GCV and PCV for the traits indicated the possibility of improving these traits through selection. Importantly, if the PCV value of the trait was higher than the corresponding GCV, this indicated that the influence of environmental factors on the expression of this trait was considerable.
Broad-sense heritability (h 2 ), which measured the proportion of the phenotypic variance that was derived from genetic effects, was calculated across two seasons using the following formula [67]: where σ 2 g , σ 2 gy , and σ 2 e are the variance components of the genotype, genotype × season interaction, and residual, respectively, y is the number of seasons, and r is the number of replications.
Genetic advance (GA) and genetic gain (GG) were estimated using the following formulae: where k is the selection differential (the value of which was 2.06 at a 5% selection intensity) and X is the phenotypic mean for destructively measured traits and SRIs.
To better view the interrelationships between all destructively measured traits and SRIs, to classify the contribution of different SRIs rather than each SRI individually, and to classify a large number of variables into major components, the mean values of all destructively measured traits and SRIs for all genotypes over two seasons were analyzed using principal component analysis (PCA). The PCA was performed based on the data of the correlation matrix to reduce the dimensions of the data space. The first two components of PCA explained always as much of the variability in the data as possible. The XLSTAT statistical package software (vers. 2019.1, Excel Add-ins soft SARL, New York, NY, USA) was used to analyze and draw a biplot of the PCA.
Partial least squares regression (PLSR) coupled with leave-one-out cross-validation (LOOCV) were used to relate the input variables (all SRIs) to the target variables (destructively measured traits). PLSR is a powerful tool that can effectively handle data when the number of input variables is much larger than that of the target variables, and the collinearity and noise in the input variables are strong [68,69]. Therefore, this tool maximizes the covariance between the scores of the input variable and target variable because it projects both variables into a low-dimensional space.
An important step in PLSR analysis is to select the optimal number of latent factors (ONLFs) in order to represent the calibration data without over-fitting or under-fitting. In this study, the ONLFs in PLSR analysis for each destructively measured trait were decided via LOOCV. The ONLFs were selected when the predicted residual error sum of squares (PRESS) statistic was most minimized [68,70]. To increase the robustness of the results, random 12-fold cross-validation was conducted on the data from two seasons, and the maximum number of latent factors was set at 7, which was suggested by the software (Unscrambler X software Version 10.2 (CAMO Software AS, Oslo). In this study, the ONLFs, which minimized the residual error sum of squares (PRESS) statistic, were found with four ONLFs for all destructively measured traits. During model training, twenty-three (25%) out of 90 datasets were used for validation (Val.), while the remaining 67 datasets were included in the calibration (Cal.). Then, the predicative performance of the PLSR models for the Cal. and Val. datasets was evaluated using the root mean squared error (RMSE), the coefficient of determination (R 2 ), and the relative percent deviation (RPD). The best model for both Cal. and Val. was selected based on the minimum value of RMSE along with the maximum value for R 2 . A value of RPD, used to study the analytical error with variations in the dataset, more than 1.4 indicated an acceptable predicative ability of the model. Finally, the calibration models derived from two seasons' dataset of cross-validation based on all SRIs were applied to predict each destructively measured trait under each season. Linear regressions between the observed and predicted values for each destructively measured trait were performed using Sigma Plot (Sigma Plot 11.0).
To identify the most influential SRIs accounting for the most variability in each destructively measured trait, the different SRIs were further analyzed using the stepwise multiple linear regression (SMLR) method. This method combines a forward selection and a backward elimination [71] and selects the input variable (different SRIs) for the regression equation step-by-step according to the importance of the input variables. The input variable is important to the regression equation based on the influence, distribution, and the significance of a target variable (measured traits) as affected by input variables. At each step, only SRIs that remained significant at p-values of an F-statistic ≤0.05 were retained in the models. Furthermore, the SRIs were eliminated with a significant probability level set to 0.01 in backward elimination. The stepwise forward and backward process terminated if no further SRIs could be added to the model or if the SRI just entered into the model was the only SRI removed in the subsequent backward elimination [71,72].

Phenotypic Variability of Different Destructively Measured Traits and Spectral Reflectance Indices Across Genotypes
Based on ANOVA, all destructively measured traits presented highly significant differences (p < 0.001) across genotypes. Significant genotypic differences (p < 0.05, 0.01, or 0.001) were also found for most SRIs. The only exception was for the red edge inflection point (REIP) and structure insensitive pigment index (SIPI) among the vegetation-SRIs and the normalized difference moisture index-1 (NDMI-2), normalized difference water index-1240 (NDWI 1240 ), normalized difference water index-D (NDWI-D), shortwave infrared water index (SWWI), normalized difference spectral index for relative water content (NRWC), dry matter content index-D (DMCI-D), and ligno-cellulose absorption index (LCA) among the water-SRIs, for which the genotypic main effect was non-significant for theses indices ( Table 2). The season main effect and the genotype × season interaction were non-significant for all SRIs, or for all destructively measured traits except for Gs for season main effect and GLA, Gs, and HI for genotype × season interaction (Table 2).
Interestingly, over the two seasons, there was a wide range between the maximum and minimum values for either destructively measured traits or SRIs, with a few exceptions for SRIs such as REIP, SIPI, WI, and R 1000/1100 , which presented a narrow range between the maximum and minimum values among genotypes. For all destructively measured traits and most SRIs, the maximum values were two-three-times higher than the minimum values (Table 3), indicating broad genotypic differences for destructively measured traits and SRIs.   Table 1. *, **, ***, Significant at the 0.05, 0.01, and 0.001 probability levels, respectively, and ns, not significant.  Figure 1 shows the correlation matrix between different destructively measured traits and SRIs across the genotypes and seasons. Generally, the different intermediate destructively measured traits (GLA, TDW, LAI, WC, Gs, SPW, and HI) showed a high positive correlation with GY; that is, correlation coefficients (r) ranged from 0.78 to 0.94. The vegetation-SRIs were better correlated with GY and intermediate traits than the water-SRIs. The index (R 1000/1100 ) showed non-significant correlation only with TDW, SPW, and HI. In addition, the vegetation-SRIs were always positively correlated with all destructively measured traits (except REIP and SIPI) and vice versa for water-SRIs (except NDWI 1240 , SWWI, NRWC, development water balance index (WBI-D), cellulose absorption index (CAI), and normalized difference lignin index (NDLI)). Importantly, Pearson's correlation coefficients of the vegetation-SRIs that were positively correlated with GY were comparable to those of the intermediate traits correlated with GY. Furthermore, these vegetation-SRIs were also highly correlated with intermediate traits (r ranged from 0.55 to 0.85). The water-SRIs showed moderate correlation coefficients with GY and intermediate traits, with the exception of normalized difference moisture index-1 (NDMI-1), normalized difference moisture index-D (NDMI-D), NRWC, and DMCI-D, which were highly correlated with TDW; that is, r ranged from 0.72 to 0.76 for these three water-SRIs (Figure 1).

Traits
Generally, PCV was higher than the corresponding GCV for destructively measured traits and SRIs. However, the PCV was clearly higher than the GCV for GY, GS, SPW, and several SRIs, especially the water-SRIs. Furthermore, the vegetation-SRIs had higher PCV and GCV than the water-SRIs ( Table 4). The highest values for PCV (49.79%) and GCV (45.93%) were observed for PRI. Among the destructively measured traits, GY, GLA, and LAI exhibited the highest values for PCV (24.17, 29.47, and 25.39%) and GCV (21.23, 29.01, and 25.32%), respectively. The PCV and GCV for the red-edge model index (RM) and modified chlorophyll absorption ratio index (MCARI) (vegetation-SRIs) and NDWI 1240 (a water-SRI) were comparable with those of these three destructively measured traits (GY, GLA, and LAI). Interestingly, the values of PCV and GCV for several vegetation-SRIs ((normalized difference vegetation index (NDVI), red normalized difference vegetation index (RNDVI), green chlorophyll index Green chlorophyll index (GCI), simple ratio pigment index (SRPI), TVI, modified transformed vegetation index (MTVI), optimized soil adjusted vegetation index (OSAVI), and enhanced vegetation index-2 (EVI-2)) were two-three-fold higher than those of WC and HI, and some presented comparable PCV and GCV values to the destructively measured traits TDW, GS, and SPW (Table 4).
Genetic advance (GA) and genetic gain (GG) showed high variability for destructively measured traits and SRIs ( Table 4). The GA and GG for destructively measured traits ranged from 0.08 to 117.37 and from 12.41 to 58.78%, respectively ( Table 4). The GA and GG were higher for the vegetation-SRIs (GA ranged from 0.02 to 35.57 and GG ranged from 4.62 to 53.19%) than for the water-SRIs (GA ranged from 0.002 to 0.06 and GG ranged from 1.16 to 18.79%) ( Table 4).

PCA of Different Destructively Measured Traits and SRIs
A PCA was conducted for all destructively measured traits, SRIs, genotypes, and seasons simultaneously (Figure 2). The first and second PCA explained 72.04 and 10.47% of the total variation between all traits, respectively. All destructively measured traits and vegetation-SRIs (except REIP and SIPI), as well as seven out of 15 water-SRIs were grouped together in a positive direction, and most identified the drought-tolerant parent Sakha 93. Eight out of 15 water-SRIs were spread in a negative direction, and the drought-sensitive parent Sakah 61 was favored by these water-SRIs. Importantly, the angle between the vectors of the destructively measured traits and all vegetation-SRIs (except REIP and SIPI) was acute (less than 90 • ), which indicate positively associated, while the angle between the vector of the destructively measured traits and most water-SRIs (eight out of 15 water-SRIs) was higher than 90 • , which indicated negatively associated (Figure 2). Remote Sens. 2020, 12, x FOR PEER REVIEW 13 of 24 most identified the drought-tolerant parent Sakha 93. Eight out of 15 water-SRIs were spread in a negative direction, and the drought-sensitive parent Sakah 61 was favored by these water-SRIs. Importantly, the angle between the vectors of the destructively measured traits and all vegetation-SRIs (except REIP and SIPI) was acute (less than 90°), which indicate positively associated, while the angle between the vector of the destructively measured traits and most water-SRIs (eight out of 15 water-SRIs) was higher than 90°, which indicated negatively associated (Figure 2).  Table 5 summarizes the R 2 and corresponding RMSE and RPD of the calibration (Cal.) and validation (Val.) datasets of the PLSR analysis to estimate destructively measured traits based on all SRIs. Figure 3 also shows a 1:1 scatter plot of the measured and predicted values of different measured traits for PLSR analysis under each season. The PLSR model provided accurate estimations of GY, TDW, and LAI in both the Cal. and Val. datasets, with R 2 ranging from 0.68 to 0.75 in the Cal. and from 0.65 to 0.72 in the Val. datasets ( Table 5). The PLSR model generated a moderate estimation performance of GLA, RWC, Gs, SPW, and HI, with R 2 cal/val = 0.54/0.46, 0.61/0.53 0.57/0.51, 0.63/0.60, and 0.53/0.47 for GLA, RWC, Gs, SPW, and HI, respectively ( Table 5). The RPD values for all destructively measured traits in Cal. and Val. were higher than 1.45, with the exception of GLA, GS, and HI, the RPD of which was less than 1.45 in Val. only ( Table 5).

Multivariate Statistical Analysis to Estimate Different Destructively Measured Traits
The vegetation-SRIs, water-SRIs, or both SRIs were analyzed by stepwise multiple linear regression (SMLR) in order to identify the most influential SRIs to better estimate the most variation  Table 5 summarizes the R 2 and corresponding RMSE and RPD of the calibration (Cal.) and validation (Val.) datasets of the PLSR analysis to estimate destructively measured traits based on all SRIs. Figure 3 also shows a 1:1 scatter plot of the measured and predicted values of different measured traits for PLSR analysis under each season. The PLSR model provided accurate estimations of GY, TDW, and LAI in both the Cal. and Val. datasets, with R 2 ranging from 0.68 to 0.75 in the Cal. and from 0.65 to 0.72 in the Val. datasets ( Table 5). The PLSR model generated a moderate estimation performance of GLA, RWC, Gs, SPW, and HI, with R 2 cal/val = 0.54/0.46, 0.61/0.53 0.57/0.51, 0.63/0.60, and 0.53/0.47 for GLA, RWC, Gs, SPW, and HI, respectively ( Table 5). The RPD values for all destructively Remote Sens. 2020, 12, 1480 13 of 23 measured traits in Cal. and Val. were higher than 1.45, with the exception of GLA, GS, and HI, the RPD of which was less than 1.45 in Val. only (Table 5).

Multivariate Statistical Analysis to Estimate Different Destructively Measured Traits
The vegetation-SRIs, water-SRIs, or both SRIs were analyzed by stepwise multiple linear regression (SMLR) in order to identify the most influential SRIs to better estimate the most variation for each destructively measured trait. In general, the SMLR model based on vegetation-SRIs or all SRIs (R 2 ranged from 0.63 to 0.86 for both models) performed better for estimating the destructively measured traits than those based on water-SRIs (R 2 ranged from 0.41 to 0.75) ( Table 6). The NDVI performed better for estimating the destructively measured traits and explained 63.0, 75.0, and 63.0% of the total variation in Gs, SPW, and HI, respectively. The MCARI explained 83.0 and 78.0% of the variation in GY and LAI, respectively, while when combined with NPQ, it explained 86.0% of the variation in TDW and 71.0% of the variation in WC. About 70.0% of the variation in GLA was explained by the combination of MCARI and NPQ with GNDVI (Table 6).
When the water-SRIs were incorporated into the SMLR model, eight (NDMI-2, NDMI-D, NDWI -1240 , NDWI-D, SWWI, WBI-D, DMCI-D, and CAI) were identified as the most influential SRIs explaining most of the variation in different destructively measured traits (from 41.0% of the variation in GLA to 78.0% of the variation in GY or TDW) ( Table 6). When all SRIs were incorporated in the SMLR model, NDVI, GNDVI, NPQ, and MCARI from vegetation-SRIs and NDWI -1240 , DMCI-D, and NDLI from water-SRIs were identified as the most influential SRIs explaining most of the variation in different destructively measured traits (63.0-86.0% of the total variation in the measured traits; Table 6). Table 5. Coefficient of determination (R2), root mean squared error (RMSE), and relative percent deviation (RPD) for calibration (R 2 cal , RMSE Cal , and RPD cal ) and validation (R 2 val , RMSE Val , and RPD val ) statistics of partial least squares regression models based on all spectral reflectance indices for the assessment of different destructively measured traits. Estimates were calculated across two seasons. variation in GLA to 78.0% of the variation in GY or TDW) ( Table 6). When all SRIs were incorporated in the SMLR model, NDVI, GNDVI, NPQ, and MCARI from vegetation-SRIs and NDWI-1240, DMCI-D, and NDLI from water-SRIs were identified as the most influential SRIs explaining most of the variation in different destructively measured traits (63.0-86.0% of the total variation in the measured traits; Table 6).     Table 1.

Discussion
Because the final grain yield (GY) always presents a high genotype × environment interaction and suffers from low heritability, especially under environmental stresses [73,74], plant breeders rely on several intermediate traits as indirect screening criteria to evaluate the drought tolerance of genotypes and increase the chances of capturing genotypes that are well adapted to water shortage under arid conditions [75,76]. Intermediate traits, such as different agronomical traits (biomass, yield components, leaf area, and harvest index) and physiological traits (traits related to plant water status and photosynthetic efficiency) reflect the integration of several plant processes either over time or in the whole plant canopy. For example, above-ground biomass (TDW), the leaf area index (LAI), and green leaf area (GLA) are important traits that reflect radiation use efficiency, photosynthesis and evapotranspiration rates, plant competition, and the status of crop growth under particular growing conditions [77][78][79][80]. The harvest index (HI) can support an understanding of the close relationships between sources and sink capacity [81,82]. Stomatal conductance (Gs) reflects the integration between photosynthesis and transpiration capacity and supports an understanding of the close relationships between external environmental conditions and internal plant characteristics induced by water stress [82,83]. Therefore, plant breeders have investigated the direct use of these intermediate traits as secondary traits for selection as an alternative to grain yield. However, selection accuracy based on these secondary traits will be efficient if these traits are well correlated with GY and are highly heritable [84,85].
In this study, there were significant differences (p < 0.001) in all destructively measured intermediate traits across genotypes (Table 2), with a wide range between the maximum and minimum values for each trait (Table 3), indicating the existence of genotypic variation for these traits. In addition, all traits were highly positively correlated with GY (r ranged from 0.78 to 0.94, Figure 1). Importantly, all traits had high heritability (h 2 > 60%) coupled with genetic gain (>20%), except for WC and HI, and provided almost equal values for both GCV and PCV ( Table 4). The results of these genetic analyses for destructively measured intermediate traits suggested they could be used as secondary traits to improve wheat production under water shortage conditions. High heritability along with high genetic advance indicated that these traits were tightly controlled by additive gene effects and were less affected by the environment; thus, these traits could be effectively estimated based on phenotypic selection [84].
However, phenotypic monitoring of these traits using traditional methods, which rely on destructive sampling, requires high labor intensity and time and poses challenges for the widespread incorporation of secondary traits in breeding programs, especially when plant breeders wish to evaluate a sufficient number of genotypes across different growth stages, years, and locations. In addition to the importance of heritability characteristics and genetic analyses for secondary traits, cost efficiency and measurement time must be considered when collecting phenotypic data for these traits. Therefore, in the following section, we will discuss the ability of different SRIs, as rapid, non-destructive, and cost-efficient tools, for the indirect phenotypic estimation of GY and different intermediate traits or for testing SRIs as indirect alternative breeding/selection tools for destructively measured intermediate traits and GY.

Ability of Different SRIs as Alternative Indirect Selection Tools for GY and Intermediate Traits
Two types of SRIs can be used as alternative indirect selection tools for GY and intermediate traits in breeding programs, including vegetation-SRIs and water-SRIs. The first type of SRIs is constructed from different wavelengths from the visible (VIS: 400-700 nm) and red-edge (NIR: 700-850 nm) regions. They are mainly sensitive to changes that take place in the leaf chlorophyll and associated pigment content and internal leaf structural factors and have been shown to be well associated with GY and different intermediate traits such as aboveground biomass, LAI, vegetative vigor, photosynthetic radiation use efficiency, stomatal conductance, different yield components, HI, and relative water content [5,9,[23][24][25][26]29,30,36,78,86]. Conversely, the second type of SRIs, which are formulated from different wavelengths from the near-infrared (NIR: 850-1300 nm) and shortwave-infrared (SWIR: 1300-2500 nm) regions, is sensitive to changes in plant water status and leaf biochemical characteristics such as lignin, cellulose, and sugar and have been reported to be well correlated with aboveground biomass and different traits related to plant water content [16,23,36,87,88]. These close relationships between plant characteristics and SRIs indicate that the use of different SRIs for the indirect assessment of GY and different intermediate traits, or their use as alternative indirect selection/breeding tools for these traits, is possible. To confirm this, heritability and genetic analyses were performed for different SRIs in this study. Significant genotypic variation was observed for most of the vegetation-SRIs (except REIP and SIPI) and seven out of 15 water-SRIs (Table 2), with a wide range between the maximum and minimum values (Table 3), confirming the existence of sufficient genotypic variation for these SRIs. The vegetation-SRIs showed high positive phenotypic correlations, while the water-SRIs showed negative moderate phenotypic correlations with GY and different intermediate traits (Figure 1). In addition, the vegetation-SRIs showed high heritability, while water-SRIs showed moderate to high heritability (Table 4). Most vegetation-SRIs had high genetic gain (>20%), while the genetic gain values for all water-SRIs were less than 20% (Table 4). Most vegetation-SRIs exhibited higher values of PCV and GCV than water-SRIs, and most importantly, the former SRIs provided almost equal values for both GCV and PCV. Conversely, PCV for the latter SRIs was higher than the corresponding GCV, with few exceptions (Table 4), which indicated the influence of environmental factors on the expression of water-SRIs. The results of these genetic analyses for two types of SRIs revealed that the vegetation-SRIs were more highly associated at the genetic level than the water-SRIs. In addition, high heritability along with high genetic advance for vegetation-SRIs and a few water-SRIs indicated that the variation in these SRIs among genotypes was associated with a high degree of additive effect; therefore, such SRIs could be effectively estimated based on phenotypic selection [84]. Thus, these SRIs could be used as alternative indirect selection/breeding tools for GY and several intermediate traits in different breeding lines under water shortage conditions. In addition, effective genetic gain in GY or intermediate traits could be obtained by selecting genotypes based on the SRIs, especially those showing high heritability and genetic gain. A similar result was found by Prasad et al. [24], Babar et al. [35], and Gutierrez et al. [36], who reported that it was possible to achieve significant genetic gains by including SRIs' measurements during selection in wheat breeding programs, practically if this selection is performed using mid-and late-breeding generations (F5 to F7).
To provide a comprehensive picture of the interrelationships between all destructively measured traits, SRIs, and tested genotypes, a PCA was used to visualize this interrelationship. The results of the PCA showed that the vegetation-SRIs (except REIP and SIPI), seven out of 15 water-SRIs, and all destructively measured traits were grouped together with an acute angle between vectors (the angles between the vectors of these traits were less than 90 • ). However, the remaining water-SRIs and REIP and SIPI were grouped together, with an angle higher than 90 • between the vectors of these SRIs and those of the destructively measured traits ( Figure 2). Generally, the acute angle between the vectors of SRIs and the destructively measured traits indicated that these SRIs were strongly and positively associated with the measured traits and were able to estimate these traits indirectly. Therefore, such SRIs can be used as alternative indirect selection/breeding tools in plant breeding trials. Conversely, the angle higher than 90 • between the SRIs vectors and the measured traits indicated that these SRIs exhibited negative associations with the measured traits and may be used as complementary selection tools with others SRIs in plant breeding trials. Interestingly, all destructively measured traits and most vegetation-SRIs identified the drought-tolerant parent Sakha 93 and the RILs with higher values for the measured traits. Conversely, the drought-sensitive parent Sakah 61 and RILs that exhibited low values for the measured traits were favored by most water-SRIs. Therefore, the PCA also confirmed that different SRIs could be incorporated as rapid, non-destructive, cost-efficient, and indirect tools in wheat breeding programs during the selection process.

Multivariate Statistical Analysis to Predict the Different Destructively Measured Traits
Each SRI focuses on only two-three sensitive wavelengths, which makes some SRIs highly sensitive to biomass saturation, soil background, and growth conditions effects. Therefore, different PLSR models based on full-spectrum wavelengths, multiple SRIs, or combined SRIs and destructively measured traits have been considered to increase the accuracy of estimating relevant measured traits [9,12,22,41,78,89]. For example, ElSayed et al. [22] reported that the accuracy of GY and canopy water content estimation under water stress conditions was increased when various SRIs were used simultaneously for PLSR models. Li et al. [12] also reported that incorporating the most significant SRIs (chosen via variable importance in projection) in the PLSR model increased the efficiency of LAI estimation. In addition, Yue et al. [27] found that a PLSR model based on 54 various SRIs provided the highest accuracy in leave-one-out cross-validation for estimation of the above-ground biomass of winter wheat. In the present study, when all SRIs were combined into a single dataset in PLSR to estimate different measured traits, the results showed that the PLSR models exhibited a high to moderate R 2 in both calibration (Cal.) and validation (Val.) datasets (R 2 ranged from 0.53 to 0.75 in Cal. and from 0.46 to 0.72 in Val.) ( Table 5). In addition, it provided a good relationship between observed and predicted values for the measured traits (Figure 3), whereas the range of R 2 values and the degree of accuracy between the observed and predicted values were strongly dependent on the measured traits and season. These findings indicated that PLSR models based on various SRIs may be used to estimate indirectly the GY and different intermediate traits used widely as direct screening criteria in breeding programs aiming to improve the drought tolerance of wheat genotypes in a fast and nondestructive manner. The ability of such PLSR models to estimate indirectly destructively measured traits may be related to their ability to minimize the effects of vegetation structure, the biophysical properties of plants, soil background, and atmospheric conditions because they include a wide range of sensitive wavelengths, which are directly related to the main changes in the target traits.
To identify the most influential SRIs that explained the most variability in each destructively measured trait across genotypes, all SRIs were further applied to SMLR analysis as independent variables. Interestingly, the results of this study showed that, in general, SMLR performed better than PLSR in the estimation of measured traits, especially GLA, WC, Gs, and HI (R 2 values are compared in Tables 5 and 6). In addition, the SMLR model identified NDVI, GNDVI, NPQ, and MCARI from vegetation-SRIs and NDWI 1240 , DMCI-D, and NDLI from water-SRIs as the most influential SRIs that explained much of the variation in measured traits among genotypes (explaining 0.63-0.86 of the total variability in the measured traits; Table 6). The high variability in GLA, TDW, LAI, and HI was explained by the indices belonging to vegetation-SRIs only, while the combined indices from the vegetation-SRIs and water-SRIs were the best to explain the variability in other traits (GY, SPW, WC, and Gs) ( Table 6). This indicated that the vegetation-SRIs, which combined information from blue, green, and red-edge regions and were related to the leaf chlorophyll and associated pigments contents, as well as photosynthetic efficiency, were suitable for monitoring the similar vegetation related traits within the crop canopy. Conversely, the indices from water-SRIs, which combined information from NIR and SWIR regions, could be used as complementary indices to vegetation-SRIs to assess the plant traits related to changes in plant water status, internal leaf structure, and leaf biochemical compounds successfully. Previous studies have also reported that SRIs combining information from blue, green, red, and red-edge wavebands are the most effective for phenotyping several vegetation parameters, such as above-ground biomass, LAI, GLA, and other morphological traits related to growth status [9,12,78,90,91]. Conversely, those combining information from NIR and SWIR wavebands are effective for indirectly estimating GY, WC, and different plant physiological characteristics [24,87,[92][93][94].

Conclusions
This study provides further support for the use of SRIs, for which data can be collected rapidly and cost effectively, as alternative selection tools to destructive GY and intermediate traits for enhancing the drought tolerance of spring wheat grown under limited irrigation conditions. The performance of SRIs for indirect assessment of the destructively measured traits was also evaluated using PLSR and SMLR models. The results allowed the following conclusions: 1.
The results of genetic analysis for destructively measured intermediate traits indicated that these traits could be used as direct traits to improve wheat production under water shortage conditions. 2.
Based on the genetic basis for the different SRIs, it seemed probable that rapid genetic progress could be obtained through the use of several SRIs as indirect alternative breeding/selection tools for destructively measured intermediate traits in plant breeding trials.

3.
The PLSR models based on all SRIs as a single index could be used for assessing the different intermediate traits that are used widely as direct screening criteria in plant breeding programs aiming to improve the drought tolerance of wheat genotypes. 4.
The SMLR model identified NDVI, GNDVI, NPQ, and MCARI from vegetation-SRIs and NDWI 1240 , DMCI-D, and NDLI from water-SRIs as the most influential SRIs that explained the most variability in the different destructively measured traits among genotypes.

5.
Based on all genetic parameters and multivariate analysis, the different SRIs could be incorporated into wheat breeding programs as effective indirect selection tools to enhance the genetic gains of drought tolerance in spring wheat.