Multivariate Analysis of Morpho-Physiological Traits Reveals Differential Drought Tolerance Potential of Bread Wheat Genotypes at the Seedling Stage

Drought is one of the foremost environmental stresses that can severely limit crop growth and productivity by disrupting various physiological processes. In this study, the drought tolerance potential of 127 diverse bread wheat genotypes was evaluated by imposing polyethylene glycol (PEG)-induced drought followed by multivariate analysis of several growth-related attributes. Results showed significant variations in the mean values of different morpho-physiological traits due to PEG-induced drought effects. Correlation analysis revealed that most of the studied traits were significantly correlated among them. The robust hierarchical co-clustering indicated that all the genotypes were clustered into four major groups, with cluster 4 (26 genotypes) being, in general, drought-tolerant followed by cluster 1 (19 genotypes) whereas, cluster 2 (55 genotypes) and 3 (27 genotypes) being drought-sensitive. Linear discriminant analysis (LDA) confirmed that around 90% of the genotypes were correctly assigned to clusters. Squared distance (D2) analysis indicated that the clusters differed significantly from each other. Principal component analysis (PCA) and genotype by trait biplot analysis showed that the first three components accounted for 71.6% of the total variation, with principal component (PC) 1 accounting for 35.4%, PC2 for 24.6% and PC3 for 11.6% of the total variation. Both PCA and LDA revealed that dry weights, tissue water content, cell membrane stability, leaf relative water content, root-shoot weight ratio and seedling vigor index played the most important discriminatory roles in explaining drought tolerance variations among 127 wheat genotypes. Our results conclude that the drought-tolerant and -sensitive wheat genotypes identified in this study would offer valuable genetic tools for further improvement of wheat productivity in arid and semi-arid regions during this time of unpredictable climate change.


Introduction
The bread wheat (Triticum aestivum L.) is an important cereal crop cultivated across diverse environments, ranging from warm lowlands to temperate highlands [1]. Among the cereals, wheat ranks first and second globally in terms of acreage (215.9 million ha) and more realistic, meaningful and accurate [39]. Among the multivariate techniques, the robust hierarchical co-cluster (RHCOC) approach produces a far lower clustering error rate than the conventional hierarchical clustering approaches in presence of outlying observations in the dataset [40]. PCA-biplot is one of the most effective multivariate analyses to evaluate the traits interaction and genotypic performance and extensively used to dissect the traits correlation in different crop plants [41]. Linear discriminant analysis (LDA), particularly useful in defining groups of the genotypes as prior classification criteria, identify misclassification error and measure the distance between groups, is effectively used for screening of flooding tolerant mungbean genotypes [42].
Therefore, the overall aim of this exploratory study is to evaluate a large number of wheat genotypes for drought tolerance potential based on the performance of seedling traits under PEG induced drought stress. Specific objectives are framed to-(i) determine the changes in morpho-physiological traits of the wheat genotypes exposed to PEG induced drought stress; (ii) classify wheat genotypes into different clusters using robust hierarchical co-cluster algorithm; and (iii) establish and verify the association between seedling traits and drought tolerance using different multivariate analysis tools.

Plant Materials and Stress Treatment
One hundred twenty-seven diverse wheat genotypes collected from different sources were used in this exploratory study. Among the genotypes, 14 mutant lines (material: BARI Gom 25; mutagen: 1% EMS) were collected from ACI Seed; 17 variety and 1 advanced line from Bangladesh Wheat and Maize Research Institute and Bangladesh Institute of Nuclear Agriculture; and 95 wheat accessions collected from Plant Genetic Resource Center of Bangladesh Agricultural Research Institute. The experiment was laid out in a two factor completely randomized design (CRD) with 3 replicates. The factors included 127 wheat genotypes and two levels of drought stresses simulated by adding polyethylene glycol (PEG)-6000 at two concentrations-0% (control), and 25% (w/v) having water potential, Ψ w ≈ −0.99 MPa; [43]). PEG is of a high molecular weight and cannot pass through the cell wall, therefore PEG is used to regulate the water potential in simulating the drought stress.
Uniform sized seeds of wheat genotypes were selected and surface-sterilized with 1% sodium hypochlorite for 10 min followed by washing several times with sterile distilled water. The seeds were then soaked in distilled water for 10 min and 30 seeds were sown in two sets of petriplates (11 cm diameter) filled with sterile sand moistened with distilled water for germination for five days in room temperature (28-33 • C). At the 6th day, seedlings were moved to a controlled environment chamber (Model: GC-560H, Firstek Scientific, Taiwan) maintaining 25 ± 1 • C temperature during day and night, relative humidity (RH) of 75-80%, 14 h of photoperiod with photosynthetic photon flux density (PPFD) of 200 µmol m −2 s −1 using cool-white florescent lamps. At the 9th day of seed sowing, the germination percent was calculated as n/N × 100 (where n was the number of total germinated seed; N was the number of total seeds) (Table S1) and the seedlings were then thinned to 15 in each petriplate. In one set of petriplate, 25% PEG-6000 solution was applied on days 11 and 13 (5 mL petriplate −1 ) of seeds sowing to induce drought stress. Only distilled water was applied to the other set of petri plates and considered as control. After PEG treatment is started, 7 mL of half-strength Hoagland's nutrient solution was added in petriplates every second day. Seedlings were then allowed to grow up to the 20th day of seed sowing to observe the visual effects of PEG induced drought stress. An example of the experimental setup and the differential effects of PEG-induced drought was illustrated in the Figure S1. Growth chamber temperature and RH were monitored by a digital humidity and temperature meter (Model: HD-306, HTC Instruments, Taiwan). Sources and types of 127 genotypes listed in Table S1.

Measurement of Seedling Traits
In Day 20 after seeds sowing and 9 days since first application of PEG, shoot and root length, fresh and dry weight was recorded from 10 seedlings. Shoot and root length was measured from the root-shoot junction to the tip of the longest leaf and root, respectively. Shoot and root dry weight was recorded after oven drying at 80 • C for 24 h. The root-shoot ratio was calculated as the root dry weight divided by the shoot dry weight. Tissue water content (TWC) in terms of the amount of water per unit shoot or root fresh weight, was calculated following the formula cited from Mickky [44]. TWC = (Fresh wt. − Dry wt.)/(Fresh wt.) Seedling vigor index (SVI) was calculated according to the formula proposed by Hellal et al. [22] as SVI = ((Shoot length + Root length) × Germination percent)/100 Leaf relative water content (LRWC) was estimated according to the procedure of Meher et al. [23]. Briefly, about 0.5 g of leaf sample was incubated in 100 mL of distilled water for 4 h. After that, the turgid weights of leaf samples were taken. The leaf samples were oven-dried at 80 • C for 24 h. The dry weights of the samples were taken till a constant weight was achieved.
Cell membrane stability (CMS) was determined following the procedure of Sairam et al. [45]. Briefly, leaf samples (0.1 g) were cut into uniformly sized squares and placed in test tubes containing 10 mL of deionized water in two sets. One set was kept at 40 • C for 30 min and another set at 100 • C in a boiling water bath for 15 min and their electric conductivities C 1 and C 2 , respectively, were measured by a conductivity meter (Model: EC-400L, HumanLab Instrument Co., Suwon, Korea).
Expanded leaves from 10 seedlings were collected to form a single replicate and the same repeated twice from an independent set of seedlings to obtain three biological replicates. All the measurements and assays were done in triplicate.
The stress tolerance index (STI) for all traits was calculated using the following formula used by Fernandez [46].
where X c and X s indicate the observed values of a trait in a given genotype under control and stress treatments, respectively, while X c is the average value of a particular trait examined in all genotypes under non-stress condition.

Statistical Procedure and Computer Analysis of Data Using Machine Learning Algorithms
All the statistical analyses were done using R-4.0.2 for win (http://CRAN.R-project. org/) (accessed on 23 February 2021) in Rstudio-1.3.1093 (https://rstudio.com/) (accessed on 23 February 2021). Data obtained were subjected to 2-factor (stress treatment × genotypes) analysis of variance (ANOVA) in the general linear model using the package lme4 [47] and the mean differences were compared by Tukey's HSD test using the library agricolae [48]. Differences at p < 0.05 were considered significant.
For the graphical presentation of the descriptive statistics of the traits, box and whisker plots were used. The relative magnitude of change in the trait values due to PEG-induced drought was displayed by radar plot. Boxplot and radar plots were prepared using the packages ggplot2 and fmsb, respectively, along with reshape2 [49,50].
The degree of association of the studied traits was determined by correlation coefficients among them. The correlation coefficient matrix and correlation heatmap were visualized using ggpair function of the packages GGally and ggplot2 [51].
To categorize the similar genotypes in clusters, the hierarchical co-cluster algorithm was used. The extracted clusters were distinct from each other, whereas the genotypes within each cluster were broadly similar to each other. The STIs of the studied traits were normalized and adapted to the library rhcoclust to generate robust hierarchical co-clusters and cluster heatmap [40]. Prior to cluster analysis, the number of clusters was determined using gap statistic algorithm in fviz_nbclust function of factoextra.
The principal component analysis (PCA) was used to reduce the dimensionality of the dataset without losing important information. The Eigen value, latent vectors and PCA-biplot extracted from the PCA. The PCA was carried out using the packages ggplot2, factoextra and FactoMineR [50,52].
Linear discriminant analysis (LDA) was used to determine the correctness of the prior classification of the genotypes by cluster analysis and to calculate the squared distances among the clusters. LDA was performed by using the packages MASS, tidyverse and caret [53].

Mean Variability in Seedling Traits
In the experimental setup including 127 wheat genotypes and two growing conditions, highly significant variation was unveiled among the wheat genotypes in both growing conditions for all seedling traits (Table S2). The descriptive statistics of the seedling traits were presented in the box plot ( Figure 1). The trait values were decreased significantly in PEG treated seedlings, except the root-shoot ratio (RSR), which showed a significant increase under PEG induced drought. The degree of association of the studied traits was determined by correlation coefficients among them. The correlation coefficient matrix and correlation heatmap were visualized using ggpair function of the packages GGally and ggplot2 [51].
To categorize the similar genotypes in clusters, the hierarchical co-cluster algorithm was used. The extracted clusters were distinct from each other, whereas the genotypes within each cluster were broadly similar to each other. The STIs of the studied traits were normalized and adapted to the library rhcoclust to generate robust hierarchical co-clusters and cluster heatmap [40]. Prior to cluster analysis, the number of clusters was determined using gap statistic algorithm in fviz_nbclust function of factoextra.
The principal component analysis (PCA) was used to reduce the dimensionality of the dataset without losing important information. The Eigen value, latent vectors and PCA-biplot extracted from the PCA. The PCA was carried out using the packages ggplot2, factoextra and FactoMineR [50,52].
Linear discriminant analysis (LDA) was used to determine the correctness of the prior classification of the genotypes by cluster analysis and to calculate the squared distances among the clusters. LDA was performed by using the packages MASS, tidyverse and caret [53].

Mean Variability in Seedling Traits:
In the experimental setup including 127 wheat genotypes and two growing conditions, highly significant variation was unveiled among the wheat genotypes in both growing conditions for all seedling traits (Table S2). The descriptive statistics of the seedling traits were presented in the box plot ( Figure 1). The trait values were decreased significantly in PEG treated seedlings, except the root-shoot ratio (RSR), which showed a significant increase under PEG induced drought.   PEG treated seedlings showed a substantial decrease in shoot and root traits. Due to PEG treatment, a comparatively higher decrease in shoot traits were recorded than the root traits ( Figure 1). Shoot length (SL), shoot fresh weight (SFW), shoot dry weight (SDW) and shoot tissue water content (STWC) were decreased by 29%, 36%, 20% and 4%, respectively, over control, whereas 26%, 29%, 14% and 5% decrease was recorded in root length (RL), root fresh weight (RFW), root dry weight (RDW) and root tissue water content (RTWC), respectively, by PEG induced drought treatment ( Figure 1).
Contrarily, a significant increase in root-shoot ratio (RSR) was observed in droughtstressed seedlings ( Figure 1). As a result of drought stress, 10% increase in the RSR was recorded while 25% increase was observed in the seedling vigor index (SVI) (Figure 1). Leaf relative water content (LRWC) and cell membrane stability (CMS) were decreased by 17% and 26%, respectively, due to PEG induced drought ( Figure 1).

Hierarchical Clustering and Co-Clustering of Genotypes and Traits
Prior to the cluster analysis, the number of projected clusters was determined using gap statistic algorithm ( Figure S2). Based on the variation in observed traits, 127 genotypes were grouped into four hierarchical row clusters and the traits were grouped into three by using a robust co-clustering algorithm and presented as a co-cluster heatmap ( Figure 3; list of the genotypes in each cluster is provided in Table S3). The highly similar genotypes were placed in a row cluster whereas the highly associated traits were placed in a column cluster. The SDW, LRWC and CMS formed column cluster 1. Column cluster 2 contained RDW, STWC, RTWC and RSR and the rest of the traits SL, RL, SFW, RFW and SVI were placed in cluster 3.
Among the row clusters, cluster 2 contained the highest number of wheat genotypes (55) followed by clusters 3 (27), 4 (26) and 1 (19) (Figure 3 and Table 1). In general, row cluster 1 was determined mostly by traits of column cluster 1 and then by cluster 2 traits, while row cluster 2 and 3 dominated by the traits of column cluster 2, though row cluster 4 was governed by the traits of almost all column clusters ( Figure 3). Mean stress tolerance index (STI) of SL, RL, SFW, SDW, SVI, LRWC and CMS were the highest in row cluster 4. Apart from these highest contributing traits, other seedling traits also substantially contributed to cluster 4. The genotypes in row cluster 1 maintained better STI of SDW, SVI, LRWC and CMS followed by cluster 4. The row of cluster 2 was chiefly characterized by the highest STI of RFW, RDW, STWC, RTWC and RSR and well maintained STI of SL, RL and SFW followed by cluster 4. Row cluster 3 had a better preserved STI of STWC, RTWC and RSR followed by row cluster 4 and a decent contribution of SDW.
Except BARI Gom 21 and BINA wheat-1, all the wheat varieties and advanced lines were distributed in the row cluster 1 (6 variety) and 4 (9 variety and 1 line), whereas, most of the mutant lines and wheat accessions were placed in the row cluster 2 and 3. (Figure 3 and Table S3). Plants 2021, 10, x FOR PEER REVIEW 8 of 20   Table S3. (SL-shoot length; RL-root length; SFW-shoot fresh weight; RFW-root fresh weight; SDW-shoot dry weight; RDW-root dry weight; STWC-shoot tissue water content; RTWC-root tissue water content; RSR-root-shoot weight ratio; SVI-seedling vigor index; LRWC-leaf relative water content; CMS-cell membrane stability).

Variability of the Genotypes in the Extracted Clusters
Results revealed that genotypes of cluster 4 performed significantly better under PEG treatment followed by cluster 1, while cluster 2 and 3 found to be affected substantially as a consequence of PEG-induced drought stress ( Figure 4 and Table 2). by row cluster 4 and a decent contribution of SDW.
Except BARI Gom 21 and BINA wheat-1, all the wheat varieties and advanced lines were distributed in the row cluster 1 (6 variety) and 4 (9 variety and 1 line), whereas, most of the mutant lines and wheat accessions were placed in the row cluster 2 and 3. (Figure 3 and Table S3).

Variability of the Genotypes in the Extracted Clusters
Results revealed that genotypes of cluster 4 performed significantly better under PEG treatment followed by cluster 1, while cluster 2 and 3 found to be affected substantially as a consequence of PEG-induced drought stress ( Figure 4 and Table 2). Shoot length (SL) found to be decreased by 21%, 36%, 33% and 14% in cluster 1, 2, 3 and 4, respectively, due to drought stress treatment ( Figure 4 and Table 2). The percent decrease in RL under PEG treatment compared to control was 21%, 33%, 31% and 9% in cluster 1, 2, 3 and 4, respectively. As a result of drought stress treatment, the least 25% decrease in SFW was recorded in cluster 4 followed by 28% in cluster 1, while on average 41% decrease was recorded in other two clusters, whereas RFW decreased by 23%, 33%, 33%, and 18% in cluster 1, 2, 3 and 4, respectively, compared to control.
Under PEG induced drought treatment, on average 24% decrease in SDW was recorded in cluster 2 and 3, while 17% and 12% decrease were recorded in cluster 1 and 4, respectively, compared to control (Figure 4 and Table 2). A similar pattern of the decrease was observed in RDW with 11%, 16%, 18% and 8% decrease in cluster 1, 2, 3 and 4, respectively, due to drought treatment. Shoot and root tissue water contents were also decreased in a similar trend.   Shoot length (SL) found to be decreased by 21%, 36%, 33% and 14% in cluster 1, 2, 3 and 4, respectively, due to drought stress treatment ( Figure 4 and Table 2). The percent decrease in RL under PEG treatment compared to control was 21%, 33%, 31% and 9% in cluster 1, 2, 3 and 4, respectively. As a result of drought stress treatment, the least 25% decrease in SFW was recorded in cluster 4 followed by 28% in cluster 1, while on average 41% decrease was recorded in other two clusters, whereas RFW decreased by 23%, 33%, 33%, and 18% in cluster 1, 2, 3 and 4, respectively, compared to control.
Under PEG induced drought treatment, on average 24% decrease in SDW was recorded in cluster 2 and 3, while 17% and 12% decrease were recorded in cluster 1 and 4, respectively, compared to control ( Figure 4 and Table 2). A similar pattern of the decrease was observed in RDW with 11%, 16%, 18% and 8% decrease in cluster 1, 2, 3 and 4, respectively, due to drought treatment. Shoot and root tissue water contents were also decreased in a similar trend.
Conversely, RSR was increased under PEG treatment in all clusters, recording an increase by 8%, 12%, 9%, and 5% in cluster 1, 2, 3 and 4, respectively, compared to control ( Figure 4 and Table 2). As a result of PEG treatment, the least 12% decrease in SVI was recorded in cluster 4 followed by 20% in cluster 1, while on average 33% decrease was recorded in other two clusters.
Due to PEG induced drought treatment, the least 14% decrease in LRWC was recorded in cluster 4 followed by 15% in cluster 1, while on average 19% decrease was recorded in other two clusters, whereas CMS decreased by 24%, 27%, 28%, and 23% in cluster 1, 2, 3 and 4, respectively, compared to control (Figure 4 and Table 2).

Principal Component Analysis
Principal component analysis (PCA) is a multivariate statistical analysis for examining and simplifying complex and large datasets. Based on the correlation among the traits and extracted clusters, the pattern of variation in wheat genotypes were also studied using principal component analysis (PCA) to evaluate the diversity of the genotypes and their association with the observed traits. The stress tolerance index (STI) of all studied traits were subjected to PCA. A total of 12 principal components (PCs) were obtained, but only three PCs that exhibited eigenvalues > 1 were measured as significant. The rest of the non-significant PCs (eigenvalue < 1) were not worthy of further interpretation. The values of the PCs explained all the characters influencing about 72% of the genotypic variability in PEG stress tolerance that accounted up to the first three components, while the first two PCs explained 60% of the variability (Table 3 and Figure 5). PCA-biplot showed the PC1 exhibited about 35% of the total variability and explained principally by SVI, SL, SFW, CMS, LRWC, RL, SDW and RFW (Table 3 and Figure 5). The second PC accounted for about 25% of the total variation and are mostly contributed by RSR, STWC, RFW, RTWC and SDW. The PC3 explained about 12% of total variability and are contributed by RDW, RSR and RTWC. PCA-biplot showed the PC1 exhibited about 35% of the total variability and explained principally by SVI, SL, SFW, CMS, LRWC, RL, SDW and RFW (Table 3 and Figure  5). The second PC accounted for about 25% of the total variation and are mostly contributed by RSR, STWC, RFW, RTWC and SDW. The PC3 explained about 12% of total variability and are contributed by RDW, RSR and RTWC.
A PCA biplot analysis can be utilized to select traits that can be categorized into main groups and subgroups based on homogeneity and dissimilarity. In our data set, three groups of traits were identified in the PCA biplot considering both PC1 and PC2 simultaneously ( Figure 5). The SDW, LRWC and CMS were clustered in group I, while SL, RL, RFW, and SVI were in group II; and RSR, STWC, RTWC, RFW and RDW with group III. A PCA biplot analysis can be utilized to select traits that can be categorized into main groups and subgroups based on homogeneity and dissimilarity. In our data set, three groups of traits were identified in the PCA biplot considering both PC1 and PC2 simultaneously ( Figure 5). The SDW, LRWC and CMS were clustered in group I, while SL, RL, RFW, and SVI were in group II; and RSR, STWC, RTWC, RFW and RDW with group III.
Interestingly, the PCA biplot revealed that group I traits, the major contributors in PC1, were strongly associated with genotypes of row cluster 1 and 4, while the traits of group II, also the contributors in PC1, were associated with the genotypes of row cluster 4 ( Figure 5 and Table 3). The traits of group III contributed to PC2 were found to be the most closely correlated with the genotypes of row cluster 2 and 3, however, some traits of group III (RDW, STWC and RTWC) closely linked with the genotypes of row cluster 4 ( Figure 5 and Table 3). PCA-biplot also indicated the cluster centroids (the multi-dimensional average of the cluster) and the approximation of distances among them ( Figure 5).

Linear Discriminant Analysis
LDA is used to reduce the number of dimensions (i.e., variables) in a dataset while retaining as much information as possible and to redefine groups of the genotypes as prior classification criteria. Seedling traits were ordered by the absolute size of the coefficients with the linear discriminant functions (LD) in Table 4. It was observed that the absolute coefficient for RDW (1.340) was ranked first of the discriminatory variables followed by SDW (0.729), STWC (0.640), CMS (0.538), RTWC (0.518) and LRWC (0.513) indicated the dominant role of the traits in explaining 72% variation under LD1 (Table 4). On the other hand, the largest absolute coefficients of RSR (1.238) followed by SVI (1.008), SDW (0.837) and RDW (0.589) mostly explained 18% variation in LD2. Variation in LD3 (10%) was dominantly played by RSR (1.631), STWC (1.058) and RDW (0.582). Taken together, the above seedling traits (RDW, SDW, STWC, CMS, RTWC, LRWC, RSR and SVI) played the most dominant discriminatory role in explaining the variation of the 127 wheat genotypes by stepwise linear discriminant analysis.

Verification of Cluster Grouping by LDA
The row clusters of wheat genotypes created using cluster analysis were verified with the predictive ability of linear discriminant analysis (LDA). Genotypes within the prior clusters were tested, compared and assigned in different groups based on LDA and then identified the misclassified genotypes that were re-assigned to the appropriate groups ( Table 5). Results of the LDA revealed that about 95%, 86%, 89% and 96% genotypes were correctly assigned to cluster 1, 2, 3 and 4, respectively, with an average of 90% correctness in assigning genotypes to different clusters.

Mahalanobis Distance Matrix
The Mahalanobis squared distance (D 2 ) among the clusters were calculated by LDA ( Table 6). The distance matrix revealed that cluster 2 and 4 were the most distant with 21.02 units followed by distance of cluster 1 and 2 (20.29 units). Cluster 4 assembled the genotypes that performed better in most of the characters under PEG stress followed by cluster 1. In contrast, the genotypes grouped in cluster 2 performed poorly ( Table 6). The most similar clusters in the present study were clusters 2 and 3 (distance 7.05 units).

Co-Cluster Based Selection of Genotypes
The robust hierarchical co-clustering is the robust approach for clustering as well as co-clustering row and column entities in the absence and presence of outlying observations in the dataset. The selection of the groups of genotypes depends mainly on the objectives of the breeders in the breeding program. To make the selection process more precise and convenient, we have prepared a co-cluster matrix from 4 row and 3 column clusters extracted from the robust hierarchical co-cluster algorithm (Table 7). Genotypes within the co-clusters were sorted in descending order by the transformed scores using a machine language algorithm adapted for the rhcoclust package. Table 7. Best performers of 127 wheat genotypes within different co-cluster combinations under PEG-induced drought stress.

Trait Variability under PEG-Induced Drought Stress
Due to drought stress, significant changes were observed among the wheat genotypes for all observed traits (Figure 1). Except for RSR, the average values of studied traits were decreased in PEG induced drought-stressed seedlings as compared to the control condition. Similar findings were reported by [54] in wheat plant against drought stress.
Length and fresh weight of shoot and root severely affected under PEG-induced drought stress in the present study ( Figure 1). Our results are in agreement with the findings of other researchers who have reported a significant decrease in shoot and root length, shoot and root fresh weight, and tissue water content in wheat [27][28][29]34] and in barley seedling [22] due to PEG-induced drought stress or soil drought. In this study, a higher relative decrease in the above traits of the genotypes of row cluster 2 and 3 suggested their relative sensitivity to PEG induced drought stress than the other clusters ( Figure 4 and Table 2).
The dry mass of wheat seedling is an important trait and is also affected by PEG and water-deficient stress. Dry matters accumulation a very strong parameter to realize how much biomass is gained by the seedling. We have recorded a significant but variable extent of decrease in shoot and root dry weight; with the least decrease in genotypes of row cluster 4 followed by cluster 1 (Figure 4 and Table 2). The decreasing trend in seedling dry weight was also reported by other researchers [28,29,34] who found that drought stress had a significant effect on dry matter production of seedlings. The maintenance in seedlings' dry weight under PEG-induced drought stress has been considered as a reliable drought-tolerant criterion for different plant species, including wheat [30]. A lower relative decrease in shoot and root dry weight of the genotypes of row cluster 4 and 1 indicated that these genotypes were able to endure drought stress better than the genotypes of row cluster 2 and 3 ( Figure 4 and Table 2).
The root-shoot ratio (RSR) imitates relative root and shoot growth patterns of a crop plant. In the present study, mean RSR was increased significantly due to drought stress and the increase was comparatively higher in the genotypes of row cluster 2 and 3 than the genotypes of row cluster 1 and 4 ( Figure 4 and Table 2). A higher root-shoot ratio indicates the root growth is less affected than shoot growth of the seedlings under PEG treatment. The less affected root growth under drought stress was considered as a drought escaping mechanism of wheat seedlings by other researchers [26,27,55]. Higher RSR of the genotypes of row cluster 2 and 3 in the current study was probably due to a higher reduction in SDW than RDW of the genotypes under PEG-induced drought. In accordance with this belief, the results obtained by Sani and Boureima [56] highlighted that physiological events of the root system are less sensitive to limited water, while the upward movement of the xylem sap is constrained in lower water potential media.
Seedling vigor index (SVI) is an approach to measure stress tolerance of crop at the seedling stage by considering germination percentage, shoot and root lengths together. In our study, the SVI significantly decreased due to PEG-induced drought stress. A lower relative decrease in SVI of row cluster 4 followed by row cluster 1 indicated that the genotypes in these two clusters performed better under PEG stress. Less affected shoot and root length of row cluster 4 and 1 coincided with lesser decrease of SVI ( Figure 4 and Table 2). A similar observation was reported by Duman [57] in lettuce, Radhouane [58] in pearl millet, Saha et al. [21] in wheat and Hellal et al. [22] in barley.
Leaf relative water content (LRWC) is a physiological trait that has a great importance on the screening of wheat genotypes for drought tolerance. Drought-induced reduction in the leaf relative water content has been reported in many crops including wheat [29,34,59]. In the present study, LRWC decreased significantly due to PEG stress and the decrease was more prominent in the genotypes of row cluster 2 and 3 than the row cluster 1 and 4 ( Figure 4 and Table 2). Relatively higher reduction in LRWC of row cluster 2 and 3 was an indication of the sensitivity of genotypes of those clusters to PEG-induced drought stress. Similar higher reduction in leaf relative water content in drought sensitive wheat genotypes as compared to tolerant ones has been observed earlier [29,60,61].
It is well established that water stress could greatly disturb the stability of plant cellular membranes [62]. Cell membrane stability (CMS) is an efficient physiological criterion while studying drought tolerance. The tolerant genotypes showed a minimal reduction in CMS values under drought stress. Research stated that higher RWC and CMS values were reflecting the higher capability to tolerate the drought stress [63]. The present study demonstrated that CMS decreased significantly under PEG induced drought stress with a relatively lower decrease in the genotypes of row cluster 4 than other row clusters ( Figure 4 and Table 2). A similar decrease in the CMS was also reported by Ahmed et al. [29] and Ahmed et al. [34] in wheat seedlings. It was evident that drought-tolerant genotypes maintained greater CMS, which ultimately increased the endurance of genotypes following drought treatments [36].
Taken together, genotypes of different hierarchical row clusters showed differential variation in the performance of studied seedling traits (Figure 4). Due to PEG-induced drought, genotypes of row cluster 4 exhibited a minimal change in the trait performance compared to control followed by cluster 1. Higher relative changes in the traits of row cluster 2 and 3 placed the genotypes to the drought-sensitive end. Those genotypes narrowed down the variation of performance of the studied traits in PEG that could be considered drought tolerant. Placement of wheat varieties in the row cluster 1 and 4 indicated that these varieties may have specific association with better tolerance to PEGinduced drought than most of the mutant lines and accessions ( Figure 3 and Table S3). Mutant line AS-10632 and 28 accessions were also distributed to the tolerant clusters.

Association between Traits, Genotypes and Drought Tolerance
Correlation studies illustrate the nature and degree of association between any pairs of parameters. It offers a core concept of the association among various traits, which is beneficial for plant breeders in choosing varieties having desired attributes [64]. With few exceptions, all the studied traits in the present study exhibited significant correlations with each other (Figure 2), which suggesting that change in any one of those traits correspondingly change the other traits. This indicated that these seedling traits of wheat play an important role under PEG induced drought stress conditions to determine the response of drought. It means that if one reliable trait is picked in drought stress and used as a selection criterion that will lead to affect other seedling traits for drought conditions [65].
Hierarchical cluster analysis clearly revealed that genotypes of row cluster 4 have the highest ability to tolerate PEG-induced drought stress as judged by almost all seedling traits ( Figure 3 and Table 1). The row of cluster 1 mostly contributed by dry weights and physiological traits could be placed in the second choice for further evaluation. Though row clusters 2 and 3 were less distant, but still have some potentials (STWC, RTWC and RSR) for the breeder's interest. Cluster analysis has been exploited to define the dissimilarity and grouping of the genotypes based on drought tolerance indices [66][67][68][69]. The successive linear discriminant analysis of the cluster groups exhibited that about 90% of the genotypes were correctly assigned to different clusters that means misclassification of the genotypes was lesser in our used clustering algorithm. We have used a comparatively newer R package rhcoclust for robust hierarchical co-clustering of our data. This algorithm showed far less error rate than other contemporary clustering algorithms when outlying observations present in the dataset [40]. It can also be used to create different trait-genotype clustering matrix that can allow researchers to select genotypes of their desired trait groups. We have also sorted the genotypes in a co-cluster matrix in descending order to select the best performers using this algorithm (Table 7).
PCA is a powerful statistical procedure to reduce the dimensions of the variables and to divulge constructive evidence-driven feedback from a highly correlated dataset [70]. In the PCA biplot, the cosine of the angle between the trait vectors approximates the correlation between them, where an acute angle (<90 • ) represents positive correlation, angle of >90 • indicates a negative correlation, while equivalent to 90 • angle denotes traits are independent of each other. It has been previously reported that the angles between the vectors of the traits in biplot analysis do not precisely translate into correlation coefficients [71]. However, our results clearly demonstrated that correlations of a trait pair were well coordinated with the approximation of the vector angles and contribution of the same trait pair in the PCA biplot (Figures 1 and 5). PCA biplot analysis has been used widely and effectively by other researchers for screening drought-tolerant cultivars of wheat [27,34,67,68].
Correlation study, hierarchical cluster analysis and PCA indicated that contrasting variations were present in 127 wheat genotypes due to differences in PEG induced stress tolerance and classified the genotypes into four distinct clusters. LDA then confirms the accuracy of the distribution of the genotypes into different clusters with overall correctness of about 90% (Table 5). Moreover, LDA measured the distance between the hierarchical clusters, an approximation of which was figured out as cluster centroids in the PCA-biplot ( Figure 5). LDA also depicted that RDW, SDW, STWC, CMS, RTWC, LRWC, RSR and SVI played the most discriminatory role in the classification of 127 wheat genotypes into four clusters. LDA was effectively used for the screening of flooding tolerant mungbean genotypes earlier [42]. Our results revealed that clustering of the wheat genotypes differing in PEG-induced drought stress tolerance based on the robust hierarchical co-clustering was well interpreted by the results obtained from PCA and LDA. Altogether, it is evident that 127 wheat genotypes exhibited significant variation in the PEG-induced changes in seedling traits and the multivariate analyses could be effective in the identification of the wheat genotypes of desirable traits for drought tolerance.

Conclusions
Based on seedling traits, some genotypes were found PEG-induced drought stresstolerant, and some were drought-sensitive. According to the magnitude of change in the seedling traits and the outcome of various multivariate analyses, genotypes of cluster 4 appeared as drought-tolerant trailed by genotypes of cluster 1 because these genotypes performed well in most of the seedling traits studied, whereas genotypes of cluster 2 and 3 were affirmed as sensitive to PEG-induced drought stress due to their poor growth and physiological capability under drought stress. Results of the present study will contribute to understanding the differential responses of bread wheat genotypes to PEG-induced drought stress based on the seedling traits. Furthermore, drought tolerance is not often discussed as an independent character by the plant breeders and thus by using co-cluster combinations, the breeders can effectively choose the genotypes of the trait groups they are interested in.
Supplementary Materials: The following are available online at https://www.mdpi.com/article/ 10.3390/plants10050879/s1, Figure S1: Sample of experimental setup in the sand filled petriplates. Panel A showing condition of seedlings of a tolerant genotype BD-9889 (belongs to cluster 4) in both control and 25% PEG, and panel B showing the same of a sensitive genotype BD-623 (belongs to cluster 2), Figure S2: Gap statistic showing an optimal number of clusters to be created in the hierarchical cluster analysis based on STIs of the seedling traits, Table S1: List of wheat genotypes used in the exploratory study, Table S2: Mean squares and their effects on seedling traits extracted from the ANOVA of the general linear model, and