Association Analysis for Some Biochemical Traits in Wild Relatives of Wheat under Drought Stress Conditions

In the present study, we estimated genetic diversity and population structure in 186 accessions of Triticum and Aegilops species using 24 simple sequence repeat markers (SSR). Furthermore, an association analysis was performed for antioxidant activities, including guaiacol peroxidase (GPX), ascorbate peroxidase (APX), peroxidase (POX), catalase (CAT), and dry matter (DM) under two control and drought stress conditions. Our findings showed that drought treatment significantly decreased DM, whereas activities of all antioxidant enzymes were increased compared to the control conditions. The results of correlation analysis indicated that, under drought stress conditions, all biochemical traits had a positive and significant association with each other and with dry matter. In the molecular section, the results of the analysis of molecular variance (AMOVA) indicated that the molecular variation within species is more than within them. The dendrogram obtained by cluster analysis showed that grouping the investigated accessions was in accordance with their genomic constitutions. The results of association analysis revealed 8 and 9 significant marker–trait associations (MTA) under control and drought stress conditions, respectively. Among identified MTAs, two associations were simultaneously found in both growing conditions. Moreover, several SSR markers were associated with multiple traits across both conditions. In conclusion, our results could provide worthwhile information regarding marker-assisted selection for the activity of antioxidant enzymes in future breeding programs.


Introduction
Water availability is one of the most important parameters in determining global agricultural crop productivity [1]. Among abiotic stress, drought is known as one of the most important effects-reducing crop productivity worldwide. It has been reported that inaccessibility to water is responsible for more than 50% of the variance in final plants' performance. Hence it is the most common source of performance losses throughout the world, and developing more drought-tolerant varieties is one of the key strategies for food security [2]. Plants will be faced with a water deficit at least of some duration in their life cycle; hence, they use several defense mechanisms to survive and adapt to water scarcity conditions [3]. When drought stress is prolonged, it will result in the over-generation of reactive oxygen species (ROS), which in turn cause oxidative stress [4]. ROSs Table S1). A factorial pot experiment with three replications was carried out in a research glasshouse at the Department of Agronomy and Plant Breeding, Tehran University, Karaj, Iran, in 2018-2019. After breaking seed dormancy, all samples were sown in plastic pots (40 cm height and 20 cm diameter), filled with equal parts soil, humus, and sand. After complete germination, all seedlings were grown under optimal conditions of temperature (25/15 ± 2 • C during the day and night) and photoperiod (16/8 h light/dark). After the appearance of three true leaves, drought stress was initiated based on the field capacity (FC) of each pot. The FC was determined according to the method suggested by Souza et al. [32]. Half of the pots were kept under full field capacity (FC = 95 ± 5%) as the control conditions, and the other half were subjected to drought stress (FC = 30 ± 5%) as the drought stress conditions, for 20 days. After 30 days, the shoot biomass, along with the activities of several antioxidant enzymes, were determined.

Phenotyping Assay
Thirty days after sowing and applying drought stress, leaves of tested samples were subjected to the extraction of crude enzymes based on a method proposed by Pagariya et al. [33]. For this, 100 mg of fresh leaves was homogenized in 1 mL extraction buffer. The extraction buffer contained 0.1 mM EDTA, phosphate buffer (pH 7.4), 0.1% x-Triton, and 1% polyvinyl pyrrolidone (PVP). The extracts were filtered and immediately centrifuged at 15,000× g for 20 min at 4 • C. Antioxidant enzyme activities, such as peroxidase (POD), catalase (CAT), guaiacol peroxidase (GPX), and ascorbate peroxidase (APX), were determined based on the separated supernatants as described by Manoranjan and Dinabandhu [34], Hadwan [35], Chance and Maehly [36], and Nakano and Asada [37], respectively. Finally, the shoots of each sample were harvested and oven-dried at 70 • C for 72 h. After that, the dry tissues were weighed and recorded as dry matter (DM).

Genotyping Assay
The total DNA for all investigated samples was isolated based on the CTAB protocol [38]. The quality of isolated DNAs was confirmed using 0.8% agarose gel electrophoresis. To achieve molecular data, a set of microsatellite primers was selected based on the D genome of bread wheat, which was previously developed by Roder et al. [39] (Table 1). All PCR reactions were run out in a 20 µL reaction mixture comprising 10 µL master mix 2XPCR (ready-to-use PCR master mix 2X; Ampliqon), 6 µL double-distilled water, 2 µL of template DNA, and 2 µL of each primer. Amplification reactions were carried out as follows: an initial denaturation step at 95 • C for 5 min, followed by 35 cycles of denaturation at 95 • C for 45 s, a primer-annealing temperature that varied between 51.3 and 69.3 • C for 45 s, and primer elongation at 72 • C for 60 s; the final extension at 72 • C was held for 5 min. The amplified products were visualized on a 2% agarose gel stained with safe view II using a UV-based imaging system.

Data Analysis
A general linear model (GLM) was used to assess the main [drought stress treatments (D) and accessions (A)] and interaction effects (D × A) using R software [40]. The percentage changes in antioxidant profiles and dry biomass were estimated as used by Pour-Aboughadareh et al. [9]. To identify the most tolerant individuals from each species, the stress tolerance index (STI) was calculated using the iPASTIC toolkit [41]. To group the measured traits and to study the relationships among biochemical features with dry biomass, principal component analysis (PCA) was performed using the 'factoextra', 'ggdendro', and 'ggplot2' packages of R software [40].
In the molecular section, all banding patterns obtained from PCR products were used to make the binary matrices. To test the efficiency of used primers, several informativeness parameters such as average heterozygosity (He), polymorphism information content (PIC), and gene diversity (H) were calculated. An analysis of molecular variance (AMOVA) was done with the aim of the dissection of genetic variation between and within species using the GenAlEx package ver. 6.5 [42]. To compare the level of genetic diversity in the investigated species, some genetic indices, including the number of observed alleles (Na), the effective number of alleles (Ne), Nei's gene diversity (H), Shannon's information index (I), and the percentage of polymorphic loci (PPL), were calculated using GenAlEx package [42]. The phylogenetic tree was drawn based on Jaccard's genetic similarities coefficients using MEGA ver. 5.1 software [43]. Furthermore, principal coordinate analysis (PCoA) was computed using GenAlEx package [42]. A Bayesian clustering algorithm was used to analyze population structure using STRUCTURE 2.3.4 software [44]. Using a web-based STRUCTURE HARVESTER v2.3.4 [45], the optimum number of subpopulations (∆K) was obtained. A marker-trait association analysis was performed using SSR markers and measured biochemical traits through a mixed linear model (MLM) using TASSEL ver2.1 software [46].

Phenotypic Variation
The results of ANOVA indicated significant variations in growth conditions [control and drought stress environments] and investigated samples for all measured biochemical traits and dry matter (DM). Also, the 'drought × accession' interaction effect was significant for all traits except DM. The mean values for CAT activity across all accessions increased under drought stress conditions by 79% compared with the control conditions ( Figure 1A). CAT ranged from 0.01 to 0.15 U mL −1 in the control conditions and varied between 0.03 and 0.25 U mL −1 in stress conditions. Under drought stress conditions, Ae. cylindrica showed the highest activity of CAT compared with other species ( Figure 1A). Among the tested species, Ae. crassa and Ae. cylindrica showed the maximum values for this antioxidant enzyme under drought stress conditions ( Figure 1A). POX activity ranged from 0.01 to 0.14 with an average of 0.06 U mL −1 in the control condition and from 0.05 to 0.32 with an average of 0.13 U mL −1 in the stress conditions. Indeed, drought stress increased the POX activity by 109% compared with control conditions. Among species, Ae. crassa showed the maximum value for POX under drought stress conditions ( Figure 1B).
Like other antioxidant enzymes, the activity of GXP increased from drought stress by 108%. Under control conditions, GXP varied between 4.31 and 31.95, with a mean of 10.35 U mL −1 , while under stress conditions, it ranged from 9.04 to 36.40, with an average of 21.59 U mL −1 . Furthermore, the highest activity of GPX was observed for Ae. crassa species under drought stress conditions ( Figure 1C). The APX activity was significantly affected by drought stress. This enzyme ranged from 0.03 to 0.26, with an average of 0.09 U mL −1 under the control conditions, whereas under the drought stress conditions, it varied between 0.03 and 0.29, with an average of 0.14. As a result, drought stress increased the activity of this antioxidant enzyme across all samples by 44.38% when compared with the corresponding values under control conditions. From a species viewpoint, Ae. cylindrica showed the highest APX activity under the stress conditions ( Figure 1D and Ae. tauschii species, respectively. Furthermore, only one sample from Ae. cylindrica was placed in group A. Pearson's correlation coefficients showed that DM was significantly correlated with all antioxidant activities under drought stress conditions. Hence, these biochemical traits can be used as biomarkers to identify drought-tolerant accessions ( Figure 3A). A principal component analysis (PCA) was used to assay the existence of biochemical variation among the instigated wild wheat accessions under drought stress conditions. Accordingly, the first two principal components (PCs) accounted for 64.52% (PC1: 47.86% and PC2: 16.66) of the total biochemical variation. According to the PCA-based biplot, all accessions were widely separated from each other, suggesting a high level of variability among the investigated germplasms ( Figure 3B). Like other antioxidant enzymes, the activity of GXP increased from drought stress by 108%. Under control conditions, GXP varied between 4.31 and 31.95, with a mean of 10.35 U mL −1 , while under stress conditions, it ranged from 9.04 to 36.40, with an average of 21.59 U mL −1 . Furthermore, the highest activity of GPX was observed for Ae. crassa species under drought stress conditions ( Figure 1C). The APX activity was significantly affected by drought stress. This enzyme ranged from 0.03 to 0.26, with an average of 0.09 U mL −1 under the control conditions, whereas under the drought stress conditions, it varied between 0.03 and 0.29, with an average of 0.14. As a result, drought stress increased the activity of this antioxidant enzyme across all samples by 44.38% when compared with the corresponding values under control conditions. From a species viewpoint, Ae. cylindrica showed the highest APX activity under the stress conditions ( Figure 1D). As a result, drought stress decreases the DM parameter. Under stress conditions, DM de-   Pearson's correlation coefficients showed that DM was significantly correlated with all antioxidant activities under drought stress conditions. Hence, these biochemical traits can be used as biomarkers to identify drought-tolerant accessions ( Figure 3A). A principal component analysis (PCA) was used to assay the existence of biochemical variation among the instigated wild wheat accessions under drought stress conditions. Accordingly, the first two principal components (PCs) accounted for 64.52% (PC1: 47.86% and PC2: 16.66) of the total biochemical variation. According to the PCA-based biplot, all accessions were widely separated from each other, suggesting a high level of variability among the investigated germplasms ( Figure 3B).  Pearson's correlation coefficients showed that DM was significantly correlated with all antioxidant activities under drought stress conditions. Hence, these biochemical traits can be used as biomarkers to identify drought-tolerant accessions ( Figure 3A). A principal component analysis (PCA) was used to assay the existence of biochemical variation among the instigated wild wheat accessions under drought stress conditions. Accordingly, the first two principal components (PCs) accounted for 64.52% (PC1: 47.86% and PC2: 16.66) of the total biochemical variation. According to the PCA-based biplot, all accessions were widely separated from each other, suggesting a high level of variability among the investigated germplasms ( Figure 3B).

Molecular Variation
In this section, 48 polymorphic fragments were generated using the 24 polymorphism SSR primers across all Triticum and Aegilops accessions.
The PIC values varied between 0.15 (Xgwm-272) and 0.38 (Xgwm-157, Xgwm-469, and Xgwm-325), with an average of 0.32 (Table 1). The average heterozygosity (He) was 0.69, and it ranged from 0.18 (Xgwm-272) to 1 (Xgwm-157, Xgwm-469, and Xgwm-325). The gene diversity (H) ranged from 0.16 to 0.50, with an average of 0.41. Among the utilized primers, Xgwm-16, Xgwm-296, Xgwm-301, Xgwm-325, Xgwm-469, Xgwm-157, and Xgwm-484 showed the highest gene diversity ( Table 1). The results of AMOVA revealed that the percentage of variance was higher within species than between them (64% vs. 36%). In terms of genetic variation parameters, Ae. tauschii showed the maximum values of Na, Ne, I, He, and PPL compared with other species. To investigate phylogenetic genetic relationships among all investigated accessions, cluster analysis was computed ( Table 2). Based on the rendered dendrogram from this analysis, all samples were grouped into four main clear clusters. The two first clusters separately consisted of all Ae. cylindrica and Ae. crassa accessions. All Ae. tauschii accessions (except seven samples) created the third cluster. The fourth cluster included all T. aestivum accessions along with the remaining samples of Ae. tauschii ( Figure 4A).

Molecular Variation
In this section, 48 polymorphic fragments were generated using the 24 polymorphism SSR primers across all Triticum and Aegilops accessions.
The PIC values varied between 0.15 (Xgwm-272) and 0.38 (Xgwm-157, Xgwm-469, and Xgwm-325), with an average of 0.32 ( Table 1). The average heterozygosity (He) was 0.69, and it ranged from 0.18 (Xgwm-272) to 1 (Xgwm-157, Xgwm-469, and Xgwm-325). The gene diversity (H) ranged from 0.16 to 0.50, with an average of 0.41. Among the utilized primers, Xgwm-16, Xgwm-296, Xgwm-301, Xgwm-325, Xgwm-469, Xgwm-157, and Xgwm-484 showed the highest gene diversity ( Table 1). The results of AMOVA revealed that the percentage of variance was higher within species than between them (64% vs. 36%). In terms of genetic variation parameters, Ae. tauschii showed the maximum values of Na, Ne, I, He, and PPL compared with other species. To investigate phylogenetic genetic relationships among all investigated accessions, cluster analysis was computed ( Table 2). Based on the rendered dendrogram from this analysis, all samples were grouped into four main clear clusters. The two first clusters separately consisted of all Ae. cylindrica and Ae. crassa accessions. All Ae. tauschii accessions (except seven samples) created the third cluster. The fourth cluster included all T. aestivum accessions along with the remaining samples of Ae. tauschii ( Figure 4A).  The SSR data matrix was used for estimating structure analysis. The results of structure analysis showed that the optimum ∆K was 5 ( Figure 4B), and the Q matrix was extracted for marker-trait association analysis. The association analysis between biochemical traits and SSR data was computed based on the MLM method. Under control and drought stress conditions, eight and nine significant associations were found, respectively. The coefficient Genes 2022, 13, 1491 9 of 14 of determination (R 2 ) ranged from 2.30 to 13.93% under control conditions, while it varied between 5.86 and 13.70% under stress conditions (Table 3). Under control conditions, the highest R 2 values were observed for markers Xgwm-455, Xgwm-484, and Xgwm-296, while markers Xgwm-455, Xgwm-515, and Xgwm-608 showed the highest R 2 values under drought stress conditions (Table 3). Under control conditions, some SSR markers were significantly associated with DM and enzymatic activities as follows: Xgwm-212 and Xgwm-484 with APX; Xgwm-272, Xgwm-296, and Xgwm-232 with CAT; Xgwm-349, Xgwm-455, and Xgwm-515 with DM. Under drought stress conditions, markers Xgwm-271 and Xgwm-608 were significantly associated with APX. The activity of GPX showed a significant association with Xgwm-455 and Xgwm-232 markers. The association between Xgwm-515 and POX was also significant. DM and STI were simultaneously associated with Xgwm-455 and Xgwm-515 markers ( Table 3). As a part of the results, the Xgwm-455 marker was simultaneously associated with GPX, DM, and STI under drought stress conditions. Furthermore, Xgwm-515 was significantly associated with POX, DM, and STI. Under both control and drought stress conditions, markers Xgwm-455 and Xgwm-515 showed a significant association with DM (Table 3).

Discussion
Aegilops and Triticum species are the most important relatives of bread wheat. Knowledge of genetic diversity in these natural resources can provide an opportunity for future wheat-breeding programs [47]. Despite various studies that have exhibited the excellent breeding potential of these species, their genetic potential is still mainly underutilized [48]. Hence, estimating the population structure in wheat germplasm is one of the main steps of the utilization of breeding potential.
Drought or water-deficit stress is one of the main consequences arising from changes in climate. Thus, screening wheat germplasm in different growth and development stages is an important task for improving its drought tolerance in future breeding programs [49]. Therefore, we assessed the responses of a core collection of bread wheat genotypes and wild relatives possessing the D genome to drought stress at the early stage of growth in terms of dry matter and four antioxidant enzymes. Our results showed a high level of physiological and biochemical variation among the investigated accessions (Figures 1 and 2). This result was in accordance with the previous studies on the response of the wild relatives of wheat to drought stress [11,13,18,50].
As part of the phenotypic assay, highly positive and significant associations were observed among dry matter and activities of CAT, APX, GPX, and POX enzymes under drought stress conditions ( Figure 3A). Indeed, this finding confirms the value of the worthiness of the data in the marker-trait association analysis [15]. Multivariate methods provide powerful classifiers to capture phenotypic diversity and can identify separate groups related to functional plant adaptation and regional origin [51]. In this work, principal component analysis (PCA) showed that the biochemical traits and dry matter captured 64.52% of the total variation, suggesting that the antioxidant profiles were efficient in showing the genetic variation among the investigated wheat accessions ( Figure 3B). The biplot-based PCA simultaneously indicated the interrelationships among between biochemical profile and dry matter (DM), as well as the distribution of samples based on the two first PCs. The first PC showed a positive and significant association with DM and other biochemical traits; hence, it is named the effective component under drought stress. As shown by the angles among traits' vectors, all antioxidant enzymes were positively correlated with DM. Formerly, these associations were confirmed by the Pearson's coefficients ( Figure 3B). Likewise, several reports indicated that PCA is an efficient multivariate method to display the level of phenotypic variation and intercept associations between measured traits [9,15,[52][53][54][55].
To identify the most drought-tolerant accessions, we used a 3D plot rendered using STI index and dry matter values under both control and drought stress conditions. Among drought-tolerance indices, STI has a good capability of identifying the most tolerant genotypes, and this fact has been confirmed in numerous studies [56][57][58][59]. Accordingly, out of 186 accessions, 41 samples were identified as the most tolerant accessions. All recognized samples (except one accession) belonged to Ae. crassa and Ae. cylindrica species ( Figure 2B).
Furthermore, a high level of molecular variation was also revealed among the investigated wheat accessions through SSR markers. The used primers generated two polymorphic bands across all investigated accessions (Table 1). Based on three markerinformativeness parameters (PIC, He, and H), we found that the used primers had a good efficiency for assessing the genetic diversity and population structure analyses. In other words, this finding reveals an acceptable competence of the used primers in marker-trait association analysis in the investigated germplasm. According to the result of ANOVA, a high degree of molecular diversity was found within species (Table 2). Previously, several researchers estimated the genetic diversity in wheat germplasm through various molecular marker techniques. For instance, Moradkhani et al. [60] using microsatellite markers reported a vast genetic diversity in some Aegilops populations. Pour-Aboughadareh et al. [61] using the ISSR technique revealed the high genetic diversity among Triticum boeoticum populations collected from different regions of Iran. Furthermore, Naghavi et al. [62], Pour-Aboughadareh et al. [63], Arabbeigi et al. [28], and Etminan et al. [64] reported a high level of genetic diversity among Iranian wild relatives of wheat using various markers such as ISSR, SSR, SCoT, and CBDP markers. As shown in Table 2, the genetic diversity in the Ae. tauschii population was greater than what was found among other species. This result suggests a good breeding potential of this species for exploring new genes and even alleles with the aim of utilization in breeding programs. In this regard, numerous studies mentioned the fact that Ae. tauschii can serve as an ideal breeding resource especially for improving tolerance to abiotic stresses. Ahmadi et al. [54] examined a set of the wild relatives of wheat and reported that Ae. tauschii responded well to drought stress due to its expanded root system. This species has also shown good potential to cope with salt stress when it was investigated under severe salinity conditions at the seedling stage [10].
After detecting a high level of biochemical and molecular variability in the studied germplasm, we used a MLM model-based marker-trait association analysis (MTA) to incorporate phenotypic and genotypic data [65]. The results of association analysis for the biochemical profile and other parameters indicated that the MLM model was effective in detecting significant MTAs. This finding is in agreement with previous reports that showed a high efficiency of the MLM method to perfectly present the significant association [15,16,[66][67][68]. In this study, a total of eight and nine significant SSR markers associated with DM and activities of antioxidant enzymes in the control and drought stress conditions, respectively ( Table 3). The explained coefficient of determinations (R 2 ) for the identified MTAs were high (varied between 2.30 and 13.93% in the control and between 5.86 and 13.70% in the drought stress conditions), suggesting many genes contribute to the main part of the measured quantitative traits [15].
Under control conditions, significant MTAs were as follows; Xgwm-212 and Xgwm-484 with APX, Xgwm-272, Xgwm-296, and Xgwm-232 with CAT, and Xgwm-349, Xgwm-455, and Xgwm-515 with DM. Under drought stress conditions, markers Xgwm-271 and Xgwm-608 were significantly associated with APX. Two markers Xgwm-455 and Xgwm-232 showed a significant association with GPX. The association between Xgwm-515 and POX was also significant. The DM and STI were simultaneously associated with Xgwm-455 and Xgwm-515 markers. Furthermore, the Xgwm-455 marker was simultaneously associated with GPX, DM, and STI under drought stress conditions (Table 3); this result is more desirable in breeding programs [16]. Among the significant associations, Xgwm-455 and Xgwm-515 markers revealed a significant association with DM under both growing conditions. In other words, this finding reveals that the growing conditions are not effective in these MTAs. Indeed, this result confirms that these genomic regions are associated with different growth and development features across two growing conditions. Hence, these markers are important to wheat breeding for biochemical adaptation in drought-prone environments [69].

Conclusions
Antioxidant enzymes play a key role in the defense of plants and help them to cope with the oxidative stress induced by various types of environmental stresses. The identification of molecular markers associated with these biochemical parameters can accelerate the screening of tolerant genetic materials in breeding programs. The results of this study illustrated that wild relatives of wheat have a high level of biochemical and molecular variability. Furthermore, the SSR marker system was a powerful technique for identifying significant marker-trait associations in this germplasm. Our results indicated that the interaction between growth conditions and accessions is not influenced by some MTAs. Therefore, the results obtained from the present study contribute to the knowledge of genetics and breeding for key biochemical activities in wheat.

Supplementary Materials:
The following supporting information can be downloaded at: https: //www.mdpi.com/article/10.3390/genes13081491/s1, Table S1: Passport of the genetic materials used in this study.

Conflicts of Interest:
The authors declare no conflict of interest.