Canopy spectral reflectance indices correlate with yield traits variability in bread wheat genotypes under drought stress

Location, soil, and climatic conditions

The experiments were carried out at the experimental field of Bangabandhu Sheikh Mujibur Rahman Agricultural University (24.038°N latitude, 90.397°E longitude), Gazipur, Bangladesh, in two consecutive wheat growing seasons (November to April; 2017–18 and 2018–19) (Fig. S1). This location has a typical subtropical monsoon climate. In the study area, bread wheat is typically sown between mid-November and mid-December. The texture of the field soil was silt loam with sand 26%, silt 50%, and clay 24%; and 30.6% volumetric soil water content was required to reach full soil field capacity. The measured soil water content (%) of control and irrigation-limited trials and mean air temperature throughout the reproductive stages were presented in Fig. 1. Monthly weather data of the two wheat growing seasons as well as of the past 10 years (2010–2019) averaged data at the experimental site were summarized in Table S2, and the daily rainfall pattern was presented in Fig. S2.

Experimental design and treatments

All wheat genotypes were cultivated in two trials maintaining two different moisture regimes: normal irrigated (referred to as “control”) and restricted irrigated (discontinued irrigation after 45 days of seed sowing and termed as “drought”). The experiments were repeated for two consecutive years (2017–’18 and 2018–’19). A randomized complete block design with three replications was employed in all trials.

Seed sowing and agronomic management

Seeds of all wheat genotypes were sown on 5 December and 18 November during the first (2017–’18) and second (2018–’19) growing years, respectively. Before sowing, seeds were treated with a commercial fungicide containing carboxin and thiram to reduce seedling infection and improve germination and then sown in 20-cm-apart rows at the rate of 12g m⁻² in the 1.5 m × 1 m unit plots, maintaining a target plant density of approximately 325 seedlings m⁻². The N, P, K, S, Zn, and B fertilizer dosage of 120−25−90−15−2.5−1 kg ha⁻¹ was applied in the commercial form of urea, triple superphosphate (TSP), muriate of potash (MoP), gypsum, zinc sulphate, and boric acid, respectively (Ahmmed et al., 2018). Two-thirds of the nitrogen fertilizer was given as a basal dose in conjunction with the rest of the fertilizers at final land preparation. After the first irrigation (20 days after sowing (DAS)), the remaining ⅓ urea was used as a top dress. During the whole cropping season, five flood irrigations were applied at 20, 45, 60, 70, and 80 DAS in control plots, while only the first two irrigations (20 and 45 DAS) were applied to drought treated plots to allow vegetative growth and then the irrigation was stopped. Approximately 4.35 cm of water was applied to each plot at every irrigation with a mean flow rate of 0.0652 m³ min⁻¹. Plots were kept weed-free to avoid interference during reflectance measurement. A fungicide containing Azoxystrobin and Difenoconazole was applied at a 10-day interval after 25 days after seedling emergence following label instructions to reduce the perplexing effects of stripe and leaf rust.

Crop phenology and soil water content

Canopy reflectance, canopy temperature, and leaf greenness (SPAD readings) were collected at heading (GS55), anthesis (GS65), and after 7 (GS71), 14 (GS75), and 21 (GS85) days of anthesis (DAA) (Zadoks, Chang & Konzak, 1974). Heading and anthesis of the genotypes were determined as the visibility, and flowering of 50% of the spikes in a plot, respectively. Since the phenological events differed across genotypes, a schedule was created at the time of anthesis for each genotype for the succeeding measurements of 7, 14, and 21 DAA. Soil water status (%) was recorded daily in randomly 15 selected plots at a depth of 15 cm in each of the control and drought-treated blocks, using a microprocessor-equipped soil water status monitoring device (The PMS-714 from Lutron Electronic Enterprise Co., Ltd. in Taiwan). Averaged daily soil water data were presented in Fig. 1. Using the method of Ogbaga, Stepien & Johnson (2014), the field capacity of the experimental soil was measured at the beginning of the experiment.

Canopy temperature depression measurement

The canopy temperature was measured according to the procedure of Mohi-Ud-Din et al. (2022). Briefly, a portable infrared thermometer (Model: IR-818, URCERI, USA; the distance-spot ratio of 13:1) was used to record the canopy and ambient air temperature between 11.30 a.m. and 12.30 p.m. (Fig. S3). It took about 30° angle to get to the spotted canopy, which was about 1 m away from the horizontal line. From each plot, a total of 10 records were captured. Canopy temperature depression (CTD) was calculated using the method of Fischer et al. (1998) using the following formula:

$CTD=AT-CT$where, AT and CT are the ambient air and canopy temperatures, respectively.

Flag leaf greenness and estimation of stay green

A handheld leaf greenness reader (SPAD-502, from Konica Minolta Inc. of Osaka in Japan) was used to test flag leaf greenness as the Soil Plant Analysis Development (SPAD) index, from the middle portion of 10 preselected well-developed flag leaves of each plot (Fig. S3). The stay-green (SG) feature of genotypes was then derived from SPAD values of the flag-leaf using a slightly modified method of Rosyara et al. (2007)’s approach for estimating the area under the SPAD decline curve (AUSDC):

$SG=\sum _{i=1}^{n-1}\left[{\displaystyle \frac{S_{\left(i\right)}+S_{\left(i+1\right)}}{2}}\right]\ast \left[D_{\left(i+1\right)}-D_{\left(i\right)}\right]$where, S_(i) and S_(i+1) are the successive flag leaf SPAD values, and D_(i) and D_(i+1) are days after sowing for S_(i) and S_(i+1) measurements, respectively.

Canopy reflectance measurement

From 40 cm above the plant canopy, a handheld radiometer called an MS-720 by EKO Instruments, Japan, was used to measure the plant canopy’s reflectance (Fig. S3). The radiometer was used to look at the plant canopy from a 25° angle. The radiometer can collect reading from 350 to 1050 nm of the light spectrum at a 1 nm interval, and it can do this every time. The reflectance measurements were taken between 11 a.m. and 2 p.m., when there were no shadows, clouds, or strong winds. The mean of three readings from three random places in each plot was used to calculate reflectance indices. A white reflectance panel was used to figure out how much light was coming in for every canopy reflectance measurement. The calibrated reflectance values were used to make spectral reflectance indices (SRIs) under control and drought stress at five growth stages over two growing seasons. The formulae used to calculate the SRIs are shown in Table 1.

Measurement of phenotypic and yield traits

Days to heading (DTH) were calculated as the time necessary from seed sowing to the visibility of 50% of a plot’s spikes. The plant height (PH) was established by measuring the height of the completely emerging spike from base to tip, excluding the awns. At the physiological maturity, four linear meters of plants were picked at the ground level in the center of the plot. Spikes were extracted and counted from collected samples to determine the spike density per square meter of land area (NSM). Spikes were gathered in a cotton bag, and both the straw and the spikes were sun-dried. Spikes were manually cleaned after being threshed into the bag. The grain yield (GY) was calculated by weighing the grain and adjusting it to 12% moisture content. Biological yield (BY) was defined as aboveground biomass (straw + grain + chaff) that had been sun-dried, and t ha⁻¹ was used to denote both BY and GY. The number of kernels per spike (NKS), the weight of kernels per spike (WKS), and the hundred kernel weight (HKW) of 10 main-shoot spikes collected separately from each plot and were determined by sun-drying, threshing, weighing, and counting grain numbers.

Analysis of data using a general linear model and machine learning algorithms

Statistical analyses were done with META-R (multi-environment trial analysis in R) software version 6.0 (Alvarado et al., 2020) and R-4.1.0 for the win (R Core Team, 2022).The growth phase SRIs, SG, and CTD values were averaged and along with the yield traits, fitted to a combined analysis in the linear model with META-R for the randomized complete block design. The linear model used in META-R implemented in the R package lme4 (Bates et al., 2015) that uses the restricted maximum likelihood (REML) to estimate variance components (Alvarado et al., 2020). The analysis of variance (ANOVA) was carried out using the following equation to ascertain if there was a significant effect of genotype, environment (growing season, growing condition), and genotype × environment interactions on the traits:

$${\mathrm{Y}}_{ijk}=\mu +\mathrm{G}\mathrm{e}{\mathrm{n}}_i+\mathrm{E}\mathrm{n}{\mathrm{v}}_j+\mathrm{R}\mathrm{e}{\mathrm{p}}_k\left(\mathrm{E}\mathrm{n}{\mathrm{v}}_j\right)+(\mathrm{G}\mathrm{e}{\mathrm{n}}_i\times \mathrm{E}\mathrm{n}{\mathrm{v}}_j)+{\epsilon }_{ijk}$$where Y_ijk is observation for the i-th genotype in the j-th environment (trial) and k-th replication; µ is the grand mean effect; Gen_i is the effect of the i-th genotype; Env_j is the effect of the j-th environment; Rep_k(Env_j) is the effect of the k-th replication within the j-th environment; (Gen_i × Env_j) is the interaction effect of the i-th genotype with the j-th environment; and ε_ijk is the effect of the error associated with the i-th genotype, j-th environment and k-th replication, which is supposed to be independently and identically distributed (iid) normal with mean zero and variance σ²_ε.

META-R also estimated the best linear and unbiased estimators (BLUEs) for the model simultaneously (Alvarado et al., 2020), as the estimate (adjusted arithmetic mean, ȳ_ij) of the trait of the i-th genotype in j-th environment (Searle, 1987).

$${\overline{y}}_{ij.}=\mu +\mathrm{G}\mathrm{e}{\mathrm{n}}_i+\mathrm{E}\mathrm{n}{\mathrm{v}}_j+(\mathrm{G}\mathrm{e}{\mathrm{n}}_i\times \mathrm{E}\mathrm{n}{\mathrm{v}}_j)$$

BLUEs were calculated for the control and drought environments across growing seasons and used in the subsequent statistical analyses. The differences between the mean BLUEs of different traits under control and drought environments were compared using the paired t-test and visualized in boxplots with the R packages ggplot2 and ggpubr (Wickham, 2016). The relative BLUEs (ratio of BLUE_drought to BLUE_control) were used in cluster analysis, principal component analysis (PCA), and linear discriminant analysis (LDA).

The unbiased REML estimator maximizes the likelihood functions for variance components taking into account the degrees of freedom involved in calculating the fixed effects $[{\sigma }_{REML}^2=\sum {\left(x_i-\overline{x}\right)}^2/\left(n-1\right)]$ (Chytrasari & Kartiko, 2019). Using the REML-adjusted variance components from the above model, the broad-sense heritability across environments was calculated for SRIs and yield traits, adopting the method proposed by Alvarado et al. (2020):

$$H^2={\displaystyle \frac{{\sigma }_g^2}{{\sigma }_g^2+{\sigma }_{ge}^2/n\mathrm{E}\mathrm{n}\mathrm{v}+{\sigma }_{\epsilon }^2/\left(n\mathrm{E}\mathrm{n}\mathrm{v}\times n\mathrm{R}\mathrm{e}\mathrm{p}\right)}}$$where σ²_g, σ²_ge, and σ²_ε indicate the variances of genotype, genotype × environment interaction, and error, respectively. nEnv and nRep designate number of environments and replications, respectively.

Lollipop and radar plots were created using R packages ggplot2 and fmsb, respectively, along with reshape2 (Wickham, 2016). Correlation among BLUEs of the studied parameters was visualized by the web-based MVApp application (Julkowska et al., 2019). Using the relative BLUEs of the traits (ratio of BLUE_drought to BLUE_control), hierarchical clustering, principal component analysis (PCA), and linear discriminant analysis (LDA) were performed. Heatmaps and hierarchical clusters (the distance is Euclidean and the method is ward.D2) were generated by the package pheatmap (Kolde, 2019). The number of clusters was obtained prior to cluster analysis using the gap statistic algorithm in the factoextra’s fviz_nbclust function (Fig. S4). PCA was carried out using the R packages ggplot2, factoextra, and FactoMineR (Lê, Josse & Husson, 2008). LDA was performed using Minitab 19 statistical software (Minitab, 2022). ANOVA was performed among the clusters for biological and grain yields, means were compared using the estimated marginal means (EMM) test, and boxplots were constructed along with the statistics using the R packages tidyverse (Wickham et al., 2019), ggplot2-based ggpubr (Wickham, 2016), rstatix (Kassambara, 2020), and emmeans (Lenth, 2022).

Dynamic changes in canopy reflectance

Figure S5 depicts the dynamic variations in canopy reflectance during the reproductive growth stages under drought stress and control. Both spectra have a similar typical pattern and display the usual reflectance properties of green vegetation. Canopy reflectance throughout the VIS region of the spectrum (400–700 nm) was normally low due to light absorption by leaf pigments and moderately increased by dryness, but reflectance in the NIR portion of the spectrum (700–1,050 nm) was higher and dropped dramatically under drought conditions (Fig. S5). With the progression of reproductive growth, the canopy reflectance decreased in both control and drought conditions, especially in the NIR zone of the spectra, and the decrease was substantially higher in the drought-affected plants compared to the control (Fig. S5).

Combined analysis of variance

It was demonstrated that the genotypic main effect for all SRIs and yield parameters was extremely significant (p < 0.001) (Table 2). The genotype-environment interaction impact was incredibly significant (p < 0.001) for all SRIs, although for other yield characteristics, the interaction effect was significant at least at p < 0.01, with the exception of SG, PH, WKS, and HKW, for which the interaction effect was not significant (Table 2).

The effect of growing environments on the SRIs and yield traits

The drought treatment showed significant differences in the best linear unbiased estimators (BLUE) of SRIs and yield traits (Fig. 2). Among the SRIs, NDVI, SR, PRI, GNDVI, NWI, PNSI, EVI, MSAVI, and OSAVI were the greatest under control and decreased critically (p < 0.001) due to drought stress. Conversely, NCPI, ARI, mCRI, XES, SIPI, and PSRI were the lowest in control and increased significantly (p < 0.001) under drought stress conditions. The SG, CTD, DTH, PH, yield attributes, BY and GY decreased significantly (p < 0.01) under drought conditions. Heritability in the broad sense (H²) was calculated over growing seasons for SRIs, yield, and yield contributing traits (Table 2), and all the traits showed a higher level of heritability (>60%) across environments, ranging from 62% to 98%. Yield traits had greater heritability (>87%), while most SRIs had a heritability of less than 90%.

Functional relationships between mean SRIs and yield

Both BY and GY had insignificant and irregular correlations with SRIs under control conditions; however, these correlations became significant and stronger under drought stress (Fig. 3). Among the yields, the BY exhibited a comparatively stronger association with mean SRIs under drought condition.

The BY was positively and significantly correlated (p < 0.05) with GNDVI, EVI, NWI, SR, NDVI, PRI, PNSI, MSAVI, OSAVI, CTD, and SG, and negatively correlated with NCPI, ARI, XES, SIPI, and PSRI under drought stress. The SR, NDVI, GNDVI, EVI, PRI, CTD, and SG were positively (p < 0.05) correlated with GY, while NCPI, ARI, and XES were negatively correlated under drought. The mCRI did not show any substantial correlation with either BY or GY across the growing conditions.

Association of growth stage-specific SRIs to biological and grain yield

The growth stage-specific relationship of SRIs to BY and GY revealed that the SRIs were strongly correlated with BY and GY under drought compared to control; however, the SRIs exhibited higher correlation coefficients with BY than GY for drought treatments. Generally, SRIs taken during the mid-grain filling phases (anthesis to 14DAA) offered a stronger correlation with BY and GY related to early and late reproductive development under drought stress (Table 3). During the mid-time grain filling phases, the NWI, mCRI, SIPI, PSRI, PNSI, and MSAVI measured under drought were the SRIs having insignificant correlations with GY, whereas the insignificant relationships for BY was recorded only with mCRI, and PSRI.

Association between yield and yield contributing traits

The association between SRIs and yield traits was determined after integrating the data as BLUEs from the two seasons. Stronger significant (p < 0.05) correlations among studied parameters were recorded under drought compared to the control conditions (Fig. 3). The vegetation- and water-based indices, SG, and CTD were correlated positively amongst them, whereas they correlated negatively with pigment-specific spectral indices. Both BY and GY showed significant positive correlations with vegetation- and water-based indices, SG, and CTD and negative correlations with pigment-specific indices under drought stress (Fig. 3).

BY exhibited a significant positive correlation with WKS and HKW across growing conditions, and with NKS and the number of kernels per square meter land area (NSM) under drought, whereas GY was positively correlated with WKS, HKW, and BY under both growing conditions and with NKS under drought stress (Fig. 3). Both BY and GY exhibited negative correlations with DTH under control and drought conditions. The DTH also negatively correlated with WKS and HKW, but the correlations were significant only in the drought condition (Fig. 3). A negative relationship between GY and PH was found in both growing conditions, but the association was significant only in the control condition, while PH exhibited a significant positive correlation with NSM under both growth conditions. WKS was positively correlated with NKS and HKW across the growing conditions (Fig. 3).

Grouping of the genotypes using hierarchical cluster analysis (HCA)

Using the relative BLUEs of yield characteristics and SRIs, wheat genotypes were classified into separate clusters of comparable genotypes based on HCA (method = ward.D2 and distance = Euclidean) (Fig. 4). Using the machine-learning algorithm of gap statistics, the number of predicted clusters was computed (Fig. S4). Based on the dissimilarity in studied traits and SRIs, three clusters were generated at the genotype level, and two groups were separated at the variable level and presented as a two-way cluster heatmap (Fig. 4). All the yield traits and SR, NDVI, GNDVI, NWI, PRI, PNSI, EVI, MSAVI, OSAVI, SG, and CTD were placed in the variable group 1, whereas NCPI, ARI, mCRI, XES, SIPI, and PSRI confined to group 2. Cluster II has the most wheat genotypes (26) among the row clusters, followed by clusters I (15), and III (15) (Fig. 4). Generally, cluster I was determined primarily by the variables of group 1. In contrast, cluster III was determined mainly by the parameters of group 2. However, the genotypes of cluster II exhibited a diverse pattern of variations among the variables of the two studied groups (Fig. 4).

Principal component analysis (PCA)

The PCA based on the relative BLUE values of yield traits and SRIs unveiled a notable diversity of the bread wheat genotypes (Fig. 5). There were discovered 15 principal components (PCs), although five PCs only appeared with eigenvalues greater than one were considered significant. The first five PCs explained 73% of the variation in traits across trials, the first two of which were responsible for 56% (Table S3). In the PCA-biplot, the PC1 explained 48.2% of the total variability. All yield traits (except DTH) and SRIs contributed (positively or negatively) to the yield (excluding DTH) (Fig. 5, Table S4). The second PC, on the other hand, responded for about 7.8% variation of the total, and it was mostly contributed by ARI, NSM, BY and GY (Fig. 5, Table S3). The SRIs and yield traits that were discriminated across the first PC are assumed to be more explanatory for the evaluation of the studied genotypes.

Among the SRIs, NCPI, XES, SIPI, mCRI, and PSRI were clustered together and separated from the other attributes including BY and GY (Fig. 5), which is consistent with the negative association revealed using correlation analysis (Fig. 3). These SRIs also greatly contributed to the genotypes of cluster III (Figs. 4 and 5). Similarly, except for DTH, the positions of the rest of the traits and SRIs were placed within a 90° radius comparative to BY and GY signifying the positive phenotypic association among the traits and these traits substantially contributed to the genotypes of cluster I (Figs. 4 and 5).

Verification of hierarchical clustering by linear discriminant analysis (LDA)

The prediction capacity of LDA was used to validate the hierarchical cluster analysis. The distribution of the genotypes within the clusters was verified and compared to identify the misclassified genotypes and to re-assigned them to suitable groups based on LDA. The LDA result indicated that genotypes were assigned to distinct clusters with 98% overall correctness (Table 4).

The Mahalanobis squared distance (D²) was calculated between the clusters by LDA (Table 5). The values unveiled that clusters I and III were most distantly positioned, separated by 85.21 units, followed by clusters II and III (38.77 units). Cluster I comprised the genotypes that outperformed the majority of the phenotypes and SRIs under drought stress, while Cluster II was in the second position. On the other hand, the genotypes in cluster III appeared with poor performance (Fig. 4). The shortest distance measured in this study was between clusters I and II (36.29 units).

Cluster-wise response in SRIs, yield traits, and grain yield to drought

Spectral reflectance indices which belong to group 1 in cluster analysis (SR, NDVI, GNDVI, NWI, PRI, PNSI, EVI, MSAVI, and OSAVI) declined sharply under drought, and the magnitude of decrease was lower in genotypes’ in cluster I than in clusters II and III (Figs. 4 and 6). The trend of changes over the growth stages due to drought for the above SRIs was highest in the genotypes of cluster III, intermediate in cluster II, and lowest in cluster I (Fig. S6). On the contrary, group 2 SRIs (NCPI, ARI, mCRI, XES, SIPI, and PSRI) enhanced significantly under drought-stress with a lower increment in the genotypes of cluster I compared to the genotypes of the other two clusters (Fig. 6). The increase in NCPI and XES over the growth stages was the lowest in the genotypes of cluster I, while intermediate and the highest in clusters II and III genotypes, respectively (Fig. S6).

Under drought conditions, DTH decreased by 4%, 4%, and 5% in clusters I, II, and III, respectively, whereas, 8%, 10%, and 11% decrease was recorded in PH (Fig. 7). The SG decreased by 8%, 13%, and 15% in the genotypes of clusters I, II, and III, respectively under drought condition (Fig. 7). Growth stage-specific change in the SG revealed comparatively a lower decline in the genotypes belonging to cluster I than the genotypes of clusters II and III under drought conditions (Fig. S6). The lowest decrease (16%) in CTD was recorded in the genotypes of cluster I, while far more decrease was observed in the genotypes of clusters II and III (31% and 39%) (Fig. 7). A consistent and lower decline in CTD was observed in cluster I compared to clusters II and III in all growth stages under drought stress (Fig. S6).

Genotypes of cluster I exhibited a comparatively lower decrease in NSM, NKS, WKS, and HKW than the genotypes of clusters II and III under drought conditions (Fig. 7). NSM decreased by 9%, 24%, and 27% in clusters I, II, and III, respectively, while 8%, 16%, and 25% decreases were recorded in NKS under drought conditions (Fig. 7). About 10% and 8% decrease in WKS and HKW, respectively, was exhibited by the genotypes of cluster I under drought conditions, whereas the decrease was 22% and 14% in cluster II, and 25% and 15% in the genotypes of cluster III (Fig. 7).

The BY decreased by 10%, 29%, and 33% in the genotypes of clusters I, II, and III, respectively, due to drought (Figs. 7 and 8). The GY for clusters I, II, and III ranged from 2.77 to 6.15 tons ha⁻¹, from 2.43 to 5.22 tons ha⁻¹, and from 2.08 to 5.93 tons ha⁻¹ under control, and from 2.69 to 6.06 tons ha⁻¹, from 1.68 to 3.95 tons ha⁻¹, and from 1.55 to 4.57 tons ha⁻¹ under drought, respectively (Fig. 8). The GY decreased by 7, 32, and 35% in clusters I, II, and III, respectively, due to drought (Figs. 7 and 8). The clusters and growing conditions showed significant interactions (p < 0.05) for BY and GY. The decreases in the BY and GY of clusters II and III were significant (p < 0.001) according to the EMM (estimated marginal means) test, while the decreases were not significant in the case of cluster I (Fig. 8).

Relation between SRIs and crop yield

Under drought, the significant phenotypic association between SRIs and grain yield and biological yield were detected across the growing years (Fig. 3). The vegetation- and water-based indices always gave positive correlation coefficients, while negative correlations occurred in the case of the pigment-specific spectral indices (NCPI, ARI, and XES). Similar relations were also reported for wheat in rainfed environments (Aparicio et al., 2000; Gizaw, Garland-Campbell & Carter, 2016; Royo et al., 2003). The indices based on vegetation and water assess the canopy’s health and water content, and the greater correlations between the SRIs with biological and grain yield suggest that canopy health and water status are important factors in a genotype’s productivity under drought conditions. Besides, the significant associations between SRIs and plant height, along with biological and grain yields also indicated that plant height and biomass development during stem elongation are the major yield determinants under limiting environment. Under poor available water, the greater correlations of vegetation- and water-based indices, CTD, and SG indicate that these indices are associated with canopy coolness and greenness. In turn, the positive correlation of CTD, SG, and grain yield imitates the cooler and greener canopies that contribute to the wheat productivity in water-limiting environments. In fact, time-resolved canopy temperature and leaf stay-green features were considered simultaneously under the short-term drought, and variations among the tested wheat genotypes were also revealed (Anderegg et al., 2021).

SRIs calculated at the mid-grain-fill stage (anthesis to 7 DAA) were, on average, strongly related to yields than SRIs measured at other growth stages (Table 3). Earlier research studies have indicated that the mid-grain-fill stage is more explanatory since genotypic differences in leaf area index, and as a result, SRIs appeared as modest in the heading and also in the early vegetative growth phases (Babar et al., 2006; Gizaw, Garland-Campbell & Carter, 2016). However, better associations between SRIs and yields were achieved when five growth stages were combined, which provided a more accurate reflection of total plant growth conditions in the growing cycle as determined by the SRIs. Nonetheless, several assessments throughout growth stages are preferable to exploit the response of numerous SRIs that collectively affect crop yield. Canopy reflectance information over different growth stages was used to forecast the yield of grains previously for wheat cultivated in spring (Babar et al., 2006) and winter wheat (Gizaw, Garland-Campbell & Carter, 2016; Sun et al., 2019).

Relationship among CTD, SG, PH, biomass and SRIs

Canopy temperature depression (CTD) showed a significant positive correlation with all vegetation- and water-based indices (Fig. 3), imparting better plant health and water content were associated with cooler canopies resultant of higher rates of C-fixation linked with greater stomatal conductance (Bita & Gerats, 2013; Bonari et al., 2020). Apart from the SRIs, CTD also showed strong correlations with SG, biomass, plant height, yield, and yield attributes under drought in this present study. Various studies have stipulated that the CTD under drought is strongly associated with delayed leaf senescence, leaf and stem wax content, root system’s depth and spreading, the biomass of shoot, phenology, morphological attributes, spike sterility, and canopy attributes like leaf area index, plant height, coverage to ground, and orientation of leaves and infloresence, etc. (Lepekhov, 2022). In previous investigations, a strong association between CTD and the wheat attributes that reflect plant water status was observed (Lepekhov, 2022; Thakur, Rane & Nankar, 2022). Earlier, the water-based index NWI was suggested as an alternative to CTD by Gutiérrez-Rodríguez et al. (2004). Our findings have been reaffirmed across bread wheat genotypes, and insight on these two features combined can assist plant breeders in selecting the genotypes having higher water content and cooler canopies for water-limiting environments. SG was positively linked with vegetation- and water-based indices in the current study, except PNSI, MSAVI, and OSAVI (Fig. 3). It has been proposed that this relationship could well be utilised to create a phenotyping approach for the SG traits in wheat (Gizaw, Garland-Campbell & Carter, 2016). Staying green in wheat entails not only delayed leaf senescence but also a gradual and progressive greenness declination of the inflorescence axis and spikelets (Reynolds et al., 2009). CTD, together with canopy greenness (Kumar et al., 2017), and their temporal trends (Anderegg et al., 2021) can be exploited as a crucial feature of leaves in genotype selection for drought avoidance and tolerance to low soil moisture stress situations. Introgressing vegetation- and water-based indices, low canopy temperature, and SG traits into the wheat genotypes can produce cumulative responses that promote stress tolerance of wheat (Lopes & Reynolds, 2012), which also reflected in the results of the present investigation which revealed significantly correlated responses of the vegetation- and water-based indices with CTD and SG.

Our data also demonstrated that plant height and biomass declined substantially as a result of the drought and that they are strongly linked with the vegetation- and water-based indices. The correlated response of different vegetation indices with plant height (Ryu, Jeong & Cho, 2020) and crop biomass (Prabhakara, Hively & McCarty, 2015) under different abiotic stresses was reported in the previous studies. The significant decline in the plant height and biomass under drought could be due to stem elongation being retarded before flowering, as also indicated by the major decline in the vegetation- and water-based indices before anthesis as a result of drought stress.

Phenotypic association of SG, CTD, and yield traits

Stay green is a key yield characteristic that permits plants to keep their leaves photosynthetically active and, as a result, improve grain filling even under stressful conditions (Borrell et al., 2014; Kamal et al., 2018; Zhang et al., 2019). In this study, SG maintained significantly positive correlations with parameters like NKS, WKS, HKW, BY, and GY (Fig. 3), implying its contribution to grain yields when water is scarced. Additionally tolerance against high temperatures and water limitence (Christopher et al., 2016; Kamal et al., 2018), the association between SG and NKS has also been reported by Luche et al. (2015). The decline in SG under drought indicates restrained photosynthetic activity, early onset of senescence, and shortened assimilate translocation to developing kernels, which ultimately contributed to yield loss. It has been observed that SG is a genotypic trait that permits agricultural plants to maintain their water status., cellular integrity, sustained photosynthetic capacity, and assimilate reallocation longer after anthesis, particularly under drought and high temperature (Christopher et al., 2008; Kamal et al., 2018). CTD is considered a trustworthy indicator of plant hydration and drought tolerance, and a positive value for CTD indicates that the cooler plant canopies than the atomosphere around (Balota et al., 2008). In the present investigation, CTD showed a positive correlation with PH, NSM, NKS, BY, and GY under drought condition (Fig. 3), suggesting the cooler canopy can be attributed to the morphological and yield traits, biomass, and grain yield during drought stress. Previous research indicated that under water limiting environment, the cooler canopies at anthesis are strongly correlated with greater spikelets spike⁻¹ and a greater number of grains per spike in rain fed environment (Ozturk, 2020; Thakur, Rane & Nankar, 2022).

Drought stress generally decreases the time required to onset the heading or blooming owing to the earliler reproductive development. In this study, DTH has reduced by 3 days among the wheat genotypes due to drought. A reduced number of DTH is crucial in the breeding programs for drought stress tolerance because it allows for drought escape (Lopes & Reynolds, 2012). However, in the present study, the longer vegetative period (late heading) under control and drought did not increase grain yield, which implies that the assimilates produced over the extended pre-blooming time are not fully used during the short grain filling stage (Gizaw, Garland-Campbell & Carter, 2016; Wheeler et al., 1996). Drought-susceptible genotypes onset heading earlier during drought, resulting in a shortned life span. In contrast, drought-resistant genotypes maintained heading time which is comparable to the normal condition (Majer et al., 2008).

In the present study, drought stress negatively affected PH, NSM, NKS, WKS, HKW, BY, and GY. These results are consistent with those of Pour-Aboughadareh et al. (2020) and Wasaya et al. (2021). Reduced PH may be associated with the inhibition of cell enlargement and division, as well as faster leaf senescence, which are linked to the drought-induced reduction in turgor potential and protoplasm dehydration (Anjum et al., 2011; Khakwani et al., 2012).

In this study, significant positive correlations among NKS and WKS, HKW, BY, and GY were observed (Fig. 3), suggesting that NKS could be a reliable index of drought stress tolerance. These findings are in line with Liu et al. (2015) and Bennani et al. (2022). Earlier, NKS has been identified as the key contributor to enhanced yield in bread wheat, as well as the most affected yield traits under water-scarced environments (Chen et al., 2012; Dolferus et al., 2019). Under drought stress, reduction in the NKS along with other yield traits contributed to decreasing BY and GY (Fig. 3). The decrease in NKS is the result of floral deformities, lower viability of pollens, and fertilization inhibition induced by water scarcity throughout booting, initiation of spike, and flowering phages, which eventually accelerate yield loss (Dolferus, Ji & Richards, 2011; Dolferus et al., 2019).

Classification and selection of genotypes

Instead of looking at each trait separately, we adopted multivariate techniques to look at all of them together. Cluster analysis classifies genotypes that are similar in a cluster and the extracted clusters are distinct from each other. In the present study, the hierarchical cluster analysis yielded three different clusters from 56 wheat genotypes by evaluating the contribution of the studied traits (Fig. 4). Cluster I has the best ability to withstand drought, as determined by the lower magnitude of changes in SRIs and yield traits, and hence is designated as drought-tolerant. In contrast, clusters II and III are the moderately tolerant and susceptible clursters, respectively (Figs. 5 and 6). Cluster analysis was used to establish genetic divergence and classification of drought tolerance in wheat by Grzesiak et al. (2019), Islam et al. (2021) and Mohi-Ud-Din et al. (2021). The subsequent linear discriminant analysis exhibited that genotypes were assigned to different clusters with 98% correctness, which verified there was no significant misclassification error while clustering the genotypes (Table 4). Mahalanobis distance matrix (D²) unveiled that the clusters were significantly different from one another (Table 5).

The PCA-biplots clearly visualized the relationship between the SRIs and yield traits, as well as the coordinates of the wheat genotypes in the components (Fig. 5). All of the vegetation and water-based indices were within a 90° radius of BY and GY, indicating a positive phenotypic association. Because of their physiological relationship to vegetative growth, biomass, water content, and grain yield, vegetation and water-based indices, yield attributes, and yields were clustered into one dimension (Fig. 5). The close proximity among the vegetation and water-based indices indicates strong co-linearity and potential information redundancy. Similar vicinity was also found in the case of pigment-specific indices. This analysis may aid in the development of a dimension-reduced prediction model for grain yield under drought conditions. However, further evaluation is warranted to decrease dimensionality and information redundancy by defining unique indices that reflect different physiological features. According to the previous research, many of these indices are linked to various physiological aspects of yields and yield stability, such as NDVI to biomass, pigments, and SG; NCPI to photosynthetic efficacy under stress; NWI to crop water status (Babar et al., 2006; Reynolds, Pask & Mullan, 2012).

Thus, it is evident that 56 genotypes of wheat showed substantial differencein the SRIs and yield traits under drought and the multivariate analyses effectively classified the wheat genotypes into three clusters differing in drought tolerance. Among the genotypes, those of cluster I exhibited minimal changes in SRIs and yield traits, imparting better yield stability under drought compared to the genotypes of the other two clusters (Figs. 4 and 8), which perpetuated the relative tolerance of the genotypes of cluster I to drought.