Genome-Wide Association Mapping for Stomata and Yield Indices in Bread Wheat under Water Limited Conditions

: Genome-wide association study (GWAS) was performed for stomata- and yield-related attributes with high-density Illumina 90 K Inﬁnium SNP (single nucleotide polymorphism) array in bread wheat to determine genetic potential of germplasm for scarce water resources with sustainable yield potential. Major yield and stomata attributes were phenotyped on a panel of Pakistani and foreign accessions grown in non-stressed and water shortage environments during two seasons. Highly signiﬁcant variations were shown among accessions in both conditions for examined characteristics. Water shortage conditions reduced the overall wheat yield and strong positive correlation existed among stomatal frequency, leaf venation and grain yield per plant. Population structure analyses based on 90,000 SNP data classiﬁed the accessions into four sub-populations which indicated the presence of genetic variability. Marker-trait association (MTA) analyses revealed that 422 signiﬁcant SNPs at p ≤ 10 − 3 , after crossing the false discovery rate (FDR) <0.05 threshold, were linked with examined attributes. Pleiotropic loci ( wsnp_Ex_c8913_14881924 and Tdurum_contig10598_304 ) were associated with ﬂag leaf area (FLA), stomata size (SS), stomata frequency (SF), leaf venation (LV), number of grain per spike (NGS) and grain yield per plant (GYP), which were located on chromosome 4B and 6B at the positions 173.63cM and 229.64cM, respectively, under water shortage conditions. Pleotropic loci wsnp_Ex_c24167_33416760 , wsnp_Ex_c5412_9564046 and Tdurum_contig81797_369 on chromosomes 7A, 2A and 4B at the positions 148.26cM, 261.05cM and 173.63cM, respectively, were signiﬁcantly linked with stomata and yield indices such as FLA, SS, SF, LV, NGS and GYP under normal and water shortage conditions. The current experiment not only validated several MTAs for studied indices reported in other studies but also discovered novel MTAs signiﬁcant under water shortage environments. Associated and signiﬁcant SNPs will be useful in discovering novel genes underpinning water shortage tolerance in bread wheat for producing high-yielding and drought tolerant wheat varieties to fulﬁll the wheat demand for growing populations.


Introduction
Bread wheat (Triticum aestivum L.) is one of the earliest cereals ever domesticated and is currently one of the major sources of food and feed in the world. Wheat is adapted to diverse climatic zones including drought prone areas [1]. Changes in global climate and expansion of wheat production to less optimum production zones are causing severe crop losses annually. Water stress is one of the grand challenges limiting crop growth and productivity in various parts of the world [2]. The characteristics of the plants that result in the increase of their abilities to tolerate the physiological processes can result in the increase in their development and yields. Confirming sustainable wheat yield to accomplish the requirements of a growing population under continuous variation in climate is a remarkable problem for wheat scientists and growers [3]. Therefore, wheat needs more advancement with biotic and abiotic stress tolerance. The demand of wheat is rising with the constantly growing human population and it is predicted to increase up to 40% by 2030. Therefore, there is an urgent requirement to escalate wheat yield to settle sustainable food security [4,5]. Wheat crop ranks first among other cereal crops due to its nutritional importance and more uses. The prompt growth in population and healthier ways of life have made new tasks for wheat breeders to produce new wheat genotypes with higher yield, good quality and resistance to biotic and abiotic stresses [6].
Wheat actual yield is not in accordance with its potential due to certain restriction factors, one of which is water shortage, primarily due to irregular patterns of rainfall [7]. It was reported that yield losses in wheat are 17-70% due to a shortage of water. Water shortage is the main problem of wheat yield loss. Due to this problem, the efficiency of translocation of photosynthetic activity is affected, and ultimately wheat grain weight and total yield is reduced [8]. Another study reported a 50-90% yield reduction under water shortage environments as compared to its well-watered potential in several regions of the world. The wheat plant undergoes a severe reaction to water shortage stress at the tillering, jointing, booting, anthesis and grain filling stages. The tillering phase in wheat is a crucial phase in which the plant produces tillers, spike primordia, spikelet, and floret. Water shortage environment in this stage can decrease total wheat yield by 46% [6].
Water shortage stress is a complex mechanism because it depends on several factors, for example, species of crop, intensity of water shortage, water shortage duration, and plant growth stages [9]. To create water shortage tolerant varieties, it is imperative to know the mechanisms and behavior of plants under water shortage conditions. To survive, plants in a water shortage environment can adopt more than one tolerance mechanism. In that regard there are three basic phenomena in which a plant can compete with water shortage stress: (1) escape mechanisms, (2) avoidance mechanisms, and (3) tolerance mechanisms [10,11]. In the first phenomena, a plant completes its life cycle earlier in response to water shortage stress. In the second phenomena, plants compete parallel with water shortage stress, e.g., by the closing of stomata and/or decreasing transpiration rate. In the third phenomena, a plant proceeds phases against water shortage stress by increasing the pigmentation of photosynthates and retaining the ratio of root/shoot to efficiently partition the overall assimilate [12,13]. Although breeding progress for improvement in wheat has been attained for normal water environments, much less success has been had for water shortage environments.
Due to climate change, the frequency and severity of drought stress will significantly increase in the future and pose a threat to the food security of the rapidly increasing world population [14]. Crop productivity in dry areas can be improved through appropriate exploitation of available genetic variability of crop plants to better adapt to climate change [15]. Drought-related plant traits including leaf venation, stomata size and stomatal frequency are important to select for developing drought resistant varieties/lines [16]. To identify the ploidy level of different species of plants, the frequency and the size of the stomata have been utilized as the morphological markers. Scientists observed that there was enough variation in frequency of the stomata among the ploidy levels for the leaves in Aegilops neglecta and Triticum, respectively, but it was reported that the frequency and the size of the stomata are associated with each other negatively [17]. The problem is more acute under water shortage conditions, and the production differences are large among the maximum production regions and dryland farming regions. Scientists proposed that the larger frequency of the stomata, along with the pathways of photosynthesis, are interconnected with the larger proficiency of use of water in plants that have C 4 pathway than the plants that have C 3 pathway. The mechanisms of drought tolerance are very complicated, therefore enhancing drought tolerance selection for productivity must be utilized together with all the unidentified factors that are necessary for the increasing drought tolerance [13].
Approaches to decrease this gap include the advancement of genetics for water shortage environments by recognizing causes of water shortage tolerance and subsequent introgression of genes for associated specific characteristics to domesticated wheat genotypes. The challenges to the application of such approaches in breeding schemes are the knowledge of the desired characteristics in a time-efficient and cost-effective manner in water shortage environments [18]. A higher number of veins per unit leaf area will increase leaf area, therefore greater flag leaf area played an important role in yield increase under normal as well as water stress conditions. Stomata play an important role in regulating drought condition, as gaseous exchange and transpiration rate depends upon the size of stomata aperture and its number per unit leaf area. A lower number and smaller sized stomata can help the plant to use available water more efficiently [17,19]. Similarly, stomata size and frequency also help the plant to stand well, even under shortage of water. Denser leaf venation helps plants to survive under a shortage of water by adequately consuming available moisture, thus maintaining yield [16]. Similarly, stomatal size and number of stomata also contribute to adequate water consumption, thus, these traits can be selected for when developing normal and drought resistant wheat varieties Producing maximum yield and water shortage tolerant genotypes has been a slow process, since yield and water shortage tolerance-related genes are complex and polygenic in nature. Several methods have been considered for selecting bread wheat varieties for water shortage conditions. The choice of wheat genotypes on the basis of these drought-related traits is informal, inexpensive, and less protracted [5,20].
Understanding drought tolerance mechanisms and identifying loci responsible for mediating drought tolerance are key steps for any breeding approach aimed at increasing stress tolerance induced by water deficiency in bread wheat [13,21]. Several efforts have been made by applying conventional breeding approaches to enhance wheat production under water shortage stress, but these strategies have contributed no more than a 1% increase in annual production [22][23][24]. As time goes on, wheat breeding systems will depend upon discovering genetic and molecular potential of heat and water shortage tolerance using QTLs (quantitative trait loci) and Association Mapping (AM) studies. Genome-wide association mapping studies (GWAS) with genotypic and phenotypic data of association panels has been proved to be a robust method to detect QTLs linked with target attributes [25]. They permit the use of a different set of germplasms and make available wider genomic regions/allelic frequency with a maximum resolution, without any biparental mapping population [26]. GWAS explore the genetic basis of desired characteristics and their related genes. Genome-wide association mapping is a valuable technique with the best results due to contained maximum genetic diversity (GD) and favorable recombinant alleles among the relevant panel [27]. It is useful to identify the genomic parts linked with water shortage tolerance and yield-related indices in diverse association mapping panels. In the current experiment, a diverse panel of bread wheat accessions was genotyped with a 90,000 SNPs array, as well as phenotyped under normal and water shortage conditions for two cropping cycles in 2017-2018 and 2018-2019. Novel SNPs were found to be associated with key drought tolerance-related traits. Favorable alleles controlling grain yield and drought-related traits were of vital importance, and incorporation of these alleles after validation through marker assisted selection and fine mapping could be helpful in wheat yield improvement under stress and non-stress conditions.

Plant Material
A panel of 96 Pakistani and foreign bread wheat accessions preserved at the department of Plant Breeding and Genetics (PBG), Faculty of Agriculture and Environment, The Islamia University of Bahawalpur, Pakistan was used in this study. The detailed information about genotype code, origin, name, and pedigree record can be found in the Supplementary Materials Table S1. The association panel was sown in non-stressed and water shortage conditions in the research area of PBG in triplicates with a randomized complete block design in crop seasons of two consecutive years, 2017-2018 and 2018-2019. In the non-stressed experiment, irrigation was given at three critical stages; the first irrigation was applied at tillering stage (35 DAS (Days after sowing)), second irrigation was given at booting stage (85 DAS) and third irrigation was given at milking stage (112 DAS) [28]. In the present study, water stress was induced at tillering stage by skipping the irrigation. Each genotype was sown in a 1 m long experimental unit with 3 replicates maintaining a P × P distance of 15 cm and R × R distance of 30 cm. Only 2 seeds of each genotype were dibbled/hole and a healthy seedling was kept after germination by thinning. All other agronomic practices such as fertilizer applications, hoeing, weeding, etc., were applied homogeneously to reduce the experimental error in both conditions during both seasons.

Phenotyping
At maturity, data were recorded from ten plants from each replication for flag leaf area (FLA), stomata size (SS), stomata frequency (SF), leaf venation (LV), number of grain per spike (NGS) and grain yield per plant (GYP) under non-stressed and water shortage environments. Flag leaves of fully matured selected plants of each genotype were used to estimate the FLA. Maximum length and width of these leaves were measured in cm 2 and FLA was calculated by using the formula: flag leaf area = flag leaf length × flag leaf width × 0.74 [29]. The number of stomata per unit area was measured from the top (upper) surface of the third leaf of each randomly selected plant. The leaf strips taken from the center of the leaf were immersed in Carnoy's solution (10% acetic acid, 30% chloroform and 60% ethanol) to stop stomatal activities and remove chlorophyll from the leaf tissue [30]. After 48 h, the leaf strips were removed from the solution, peeled off with a razor, and inspected with a 40× microscope to determine the frequency of stoma in the tissue. Each leaf strip was observed and the average value was taken for analysis [16,19]. The same leaf strip was also used to measure the size of the stomata under the objective lens of a 10× microscope. For three stomata of each sample, the length and width of the stomata were measured randomly using a microscope eyepiece, and the length and width were multiplied by the microscopic normalized value of 3.33 µ by 10×, then the average mean values were scored. Using a micrometer, the length and width of the epidermal cells on the surface of the same flag leaf strip at 40× were measured [8,19]. Leaf veins were observed as the number of longitudinal veins entering the microscope field of view at 10×, which is the low magnification [8,16].

Statistical Analysis
Scored data were used for pooled analysis of variance (ANOVA) procedure applying the GenStat ® Ver.17 [31] for studied attributes in non-stressed and water shortage environments. Broad-sense heritability was also calculated from recorded data average over years. Pearson's Correlation coefficients (r) were calculated to accomplish the correlation between yield-and stomata-related characteristics using the SPSS ver.23 [32] in both environmental conditions. For ANOVA and correlation analysis significance levels, α = 0.01 was used for highly significant effects and α = 0.05 was used for significant effects. The broad-sense heritability for each trait under individual stress treatment was calculated using the following equation [33]  where H 2 is the broad-sense heritability, σ 2 g is genotypic variances, σ 2 g×e is variance for genotype and environment interaction, ny is the number of years and nr is the number of replications.

Genotyping of Bread Wheat Genotypes
Each genotype was sown in polyethene bags. The fresh leaf of 15-day wheat seedlings was used for DNA extraction according to the CIMMYT Molecular Genetics Manual procedure [34]. The DNA of each genotype (70-100 ng/µL) was preserved in 96-well plates and sent for genotyping to CapitalBio ® at Beijing in China with high-density Illumina 90 K Infinium SNP array [35]. The genome-wide positions of SNPs in terms of genetic distance (cM) situated on chromosomes founded on a consensus genetic map of bread wheat in 2015 were used in the current experiment [35]. During data analysis, monomorphic SNPs, more than 20% missing SNPs, minor alleles and allelic frequency <5%, were excluded in this study.

Population Structure
Bayesian clustering procedure was used with SNPs to categorize clusters of genetically similar genotypes through the statistical analysis through STRUCTURE v.2.3. Burn-in iterations of 10 4 cycles and admixture model selection was applied. An ad-hoc method based on the online tool "Structure Harvester v0.6.93" was practiced [36] to attain high value or peak of "K" for authentication to know the STRUCTURE outcomes. We choose the K value range from one to ten and performed six independent runs to achieve reliable effects.

GWAS Analysis
GAPIT (Genome Association and Prediction Integrated Tool) was practiced with modal selections preference to check the dependability of the outcomes [37]. It is an R package that performs a genome-wide association study and genome prediction. It implements unconventional statistical methods including compressed mixed linear modal (CMLM) and CMLM-based genomics prediction selections and determines the false discovery rate (FDR) [37,38]. The threshold for describing a marker to be significant (p-values) was measured at 10 −3 or above [25,33]. The significance level for p-values were measured after Bonferroni adjustment (p = 1/n, n = total number of SNPs) after crossing the FDR <0.05 threshold [39]. To determine the relevance of the applied model for GWAS, quantile-quantile (QQ) plot was derived among the observed and expected log10(p) value. To describe the unclear correlation obtained from population structures, covariates from STRUCTURE [40] were measured as fixed effects and using the principal components through GAPIT [37,41]. The mysterious associations among genotypes were estimated using a kinship matrix in the incorporated MLM [42]. Overall, 33,210 of the functional iSelects beads chip analyses visually exhibited polymorphisms and were detected on the available genetic map [35] in the genotypes being studied.

Phenotypic Evaluation
In the present study, analysis of variances (ANOVA) for 96 bread wheat accessions showed significant influence on phenotypic variation (p < 0.01). All traits showed significant effects among studied accessions as exhibited in Table 1. Descriptive statistics data of observed characteristics in non-stressed and water shortage environments based on average data over years are presented in Table 2. Broad-sense heritability of the studied indices was calculated and is given in Supplementary Materials Stomata frequency (SF) and leaf venation (LV) exhibited mean performance of 13.25 and 7.81 under normal environments, and 12.7 and 7.54 under water shortage environment, respectively. The average value of grain yield per plant was 22.50 g and 15.4 g under non-stressed and stressed environments, respectively. Pearson's Correlation coefficient of the studied attributes based on data averaged over years under normal and water shortage environments are mentioned in Table 3. Stomata size (SS) was negatively associated with all studied traits in non-stressed and water shortage conditions. The SF was positively correlated with LV, NGS and GYP under all studied conditions. Water shortage tolerancerelated characteristics such as FLA and NGS were strongly and positively linked among themselves, while a negative correlation with stomata size under both environments was observed Table 1. The yield-related traits such as FLA, LV, SF, NGS and GYP exhibited positive correlation with each other in non-stressed and water shortage conditions. Table 1. Mean sum of squares of 96 spring wheat genotypes for studied attributes under both seasons, data averaged over environments.

Population Structure
The results from the structure harvester displayed that the uppermost peak at K = 4 depended on the rates of changes in the log probabilities of the data among successive K-values (Supplementary Figure S1) that were found on the second order derivation on the variance of the maximum probability of the model to give a specific K. Delta K exhibits only the highest clustering level and number of sub-populations in main populations. Results from STRUCTURE analysis (Supplementary Figure S2) exhibited that a total of 12 accessions (genotype 1 to genotype 10, genotype 27 to genotype 28) were included in the first group. In the second group, a total of 14 genotypes (genotype 11 to genotype 22, genotype 29 to genotype 30) were recorded. The third group included a total of 39 genotypes from genotype 34 to genotype 72. Three genotypes from genotype 31 to genotype 33 exhibited mixed genetic material from above-mentioned (2nd and 3rd) groups. The 4th group contained a total of 17 genotypes from genotype 73 to genotype 89. A total of 6 accessions from genotype 90 to genotype 96 revealed the collective genetic material from the third and fourth group [3].

Genome-Wide Marker-Trait Associations (MTA)
In this study, 33,210 high-density SNP markers from the 90 K Illumina iSelects SNP array were evaluated to perceive SNPs associated with water shortage tolerance and yieldrelated indices. Before analyses of GWAS and genomic prediction, scientists should validate and maintain genotype quality. The GAPIT provides a series of diagnostic tools to help users perform quality control on genotypes. Marker-trait associations for these indices in non-stressed and water shortage conditions were examined. A total of 422 significant SNPs were correlated with the observed characteristics, and out of those 178 and 244 MTAs were scored in non-stressed and water shortage conditions, respectively, at or above −log 10 (p < 0.0001) threshold level using MLM (mixed linear model) for yield-and stomatarelated traits (Supplementary Tables S2 and S3). Manhattan plots Figure 1A,B, Figure 2A,B, Figure 3A,B, Figure 4A,B, Figure 5A,B and Figure 6A,B present the site of significant SNPs at −log10(p) which significantly linked with the desired characteristics under studied conditions. The quantile-quantile (QQ) plot of p-values was created (Supplementary Figures S3-S14) to confirm the results of Manhattan plots as mentioned in the figures. The Y-axis is the observed negative base 10 logarithm of the p-values, and the X-axis is the expected observed negative base 10 logarithm of the p-values under the assumption that the p-values follow a uniform (0,1) distribution. The dotted lines show the 95% confidence interval for the QQ plot under the null hypothesis of no association between the SNP and the trait.

Genome-Wide Multiple Traits Loci Associations
The highest phenotypic variability of 29.30% and 28.10% occurred in A-genome represented by the significant SNPs wsnp_RFL_Contig2699_2402527 and BS00038787_51 from chromosomes 3A and 7A, having the position of 321.73cM and 369.49cM, and were associated with SF and LV in non-stressed and water shortage conditions, respectively. The significant markers, namely wsnp_Ra_c3176_5975986 and Excalibur_c35339_106 were correlated with LV in B-genome from chromosome 7B (246.93cM) and 4B (220.48cM), explaining 29.50% and 27.97% variation in non-stress and water shortage conditions, respectively. The significant SNPs, namely Excalibur_c17237_688 and Excalibur_c17237_688 located on chromosome 2D (222.41cM) were associated with LV in D-genome showing more phenotypic variability of 27.55% and 25.60%, respectively, in the examined conditions. The lowest phenotypic variation of 6.53% and 13.74% existed in D-genome as depicted by the markers Excalibur_c53680_124 and tplb0025h02_1383 on the same chromosome, 6D, with the positions 50.3cM and 293.15cM linked with NGS in normal and water shortage conditions, respectively.
In normal conditions, pleiotropic locus (BS00038787_51) at chromosome 7A on the positions 369.49cM was also linked with SS, LV, NGS and GYP. The examined characteristics such as SF, LV, NGS and GYP were also influenced by a pleiotropic locus (wsnp_Ex_c24167_33416760) on chromosome 7A at 148.26cM under normal environments. Under water shortage conditions, a pleiotropic locus (BobWhite_c44691_648) for FLA, SS, SF, LV, NGS and GYP was identified on chromosome 4B at the position 173.63cM. The studied characteristics governed by a pleiotropic locus (wsnp_Ex_c8913_14881924) were associated with all studied yield-and stomata-related traits which were located on chromosome 4B at 173.63cM under water shortage conditions as shown in Supplementary Table S3.

Phenotypic Evaluation
In this study, highly significant differences were observed among genotypes and environmental conditions, which showed the occurrence of variation and environmental effects on genotype performance in all studied traits as exhibited in Supplementary Table S2 [5,18]. An overall decline in stomatal frequency was observed under water shortage conditions in respect to non-stressed conditions. Enhanced drought stress declined the stomatal frequency which caused a reduction in transpiration rate and exchange of gases [19]. Heritability estimates provided the information about the extent of which a genetic character can be transmitted to the successive generations. In current experiments, high heritability was reported in the studied traits under normal environments, such as FLA (73.52), followed by SS (69.21), LV (68.29) and NGS (64.33) which indicated (Supplementary Table S2) that these were simply inherited traits and most likely the heritability was due to additive gene effects, and selection may be effective in early generations for these traits. Previous studies have also reported the high heritability in NGS and GYP as complex traits. The average value of grain yield per plant was 22.50 g and 15.4 g under nonstressed and stressed environments, respectively. Adverse possessions of drought stress on wheat performance in addition to genotypic variation in response to water shortage have already been observed [13,18]. It is also reported that in wheat crop the maximum loss of grain development due to water shortage environments at anthesis and post-anthesis stage ultimately reduce overall grain yield per plant [43]. The occurrence of water deficiency at booting stage directly reduced the number of spikelets per spike, which ultimately reduced the number of grains per spike and grain yield per plant. The findings agreed with our results regarding water deficit stress at the booting stage. Few motives are associated with the limit of physiological and biochemical pathways due to water shortage, which cause yield losses in wheat [21]. Drought-related plant traits including leaf venation, stomata size and stomatal frequency are important to select for developing drought resistant varieties/lines [16]. Stomata size (SS) was negatively associated with all studied traits in non-stressed and water shortage conditions in a previous study, which supports our findings [8], and showed that no association was established with other observed parameters. Leaf venation had positive association with frequency of stomata. A higher number of veins contained maximum stomatal frequency; some wheat scientists have also observed that higher stomatal frequency may lead to thick veins. As a result, the selection for stomata size is not a favorable criterion for these genotypes. The yield-related traits such as FLA, NGS and GYP showed positive association with each other in non-stressed and drought conditions and the selection of mentioned traits will be fruitful for this germplasm [28,44]. Stomata frequency (SF) and leaf venation have already been used to select tolerant varieties [45]. In wheat, grain filling is closely related to flag leaf characteristics and function. Stomata are specialized leaf epidermal cells which regulate photosynthetic CO 2 uptake and water loss by transpiration. Understanding the mechanisms controlling stomatal size, and their opening under drought, is critical to reduce plant water loss and maintain a high photosynthetic rate which ultimately leads to elevated yield [17]. The mean variability of drought-and yield-related traits was observed. Among the genotypes with variable performance in water shortage conditions, the best performers were categorized as drought tolerant genotypes. Based on these evidences, the tolerant genotypes were G6, followed by G1 and G21, as shown in Table 4. Table 4. Best performance wheat genotypes under both environments on data averaged over years.

Population Structure
The Bayesian approach in statistical software package STRUCTURE version 2.3.3 was used to assess the genetic structures of 96 bread wheat genotypes. The observed accessions were classified into four sub-groups based on molecular data. Various kinds of colors in Supplementary Figure S2 show the discrete groups and overall studied genotypes assigned into the four sub-groups. This practice was also previously used in a wheat breeding program by several experts and they attained descriptive results [20]. STRUCTURE analysis recommended that 96 genotypes originated from diverse progenies. However, the known origin indication according to preserver of 96 wheat genotypes and pedigree records, exhibited three types of populations, but genetically these genotypes were divided into four sub-populations. According to the maintained sources (Supplementary Table S1), the first group consisted of the accessions G1 to G22 from advance-breeding lines, established in PBG-UAF, while the second group had accessions G23 to G46 from a foreign source and the third group had the accessions G47 to G96 from locally cultivated Pakistani wheat varieties.
In wheat breeding programs these practices have also been used by several wheat breeders [21,46]. STRUCTURE analysis recommended the divergences among 96 bread wheat genotypes which represented the genetic resemblance within groups and genetic differences between the groups. Mainly, outcomes were practically deliberating to the already identified pedigree record and origin of wheat accessions. Determination of genetic diversity would be helpful to recognize the diverse accessions for the improvement and progress of future wheat breeding programs [47]. Those accessions having diverse genetic makeup can be designated for required combinations to produce multiple and significant traits to gain a better yield [3]. Discernment of wheat accessions based on their genetic basis will be beneficial for effective and early screening of anticipated accessions in wheat breeding programs for producing high-yielding wheat varieties.

Genome-Wide Marker-Traits Associations
Genes and QTLs related to drought tolerance and yield characteristics in whole wheat genomes were distributed across 21 chromosomes described by several wheat breeders [33]. A marker-trait association (MTA) study recognized the connection between a particular morphological and genetic variation within a genome, which ultimately perceived locus underpinning-related characteristics at the end [48]. In this study, 33,210 high-density, polymorphic SNP markers from 90 K Illumina iSelects SNPs array [35] were examined to study SNPs associated with stomata-and yield-associated indices. Markertrait associations (MTA) for these characteristics in non-stressed and water shortage conditions were detected. The red horizontal line on Manhattan plot entitles the threshold level (p < 0.0001) of significance for SNPs with specific traits ( Figure 1A,B, Figure 2A,B, Figure 3A,B, Figure 4A,B, Figure 5A,B and Figure 6A,B). The quantile-quantile (QQ) plot (Supplementary Figures S3-S14) is a useful tool for assessing how well the model used in GWAS accounts for population structure and familial relatedness. In this plot, the negative logarithms of the p-values from the models fitted in GWAS are plotted against their expected value under the null hypothesis of no association with the trait. Because most of the SNPs tested are probably not associated with the trait, the majority of the points in the QQ plot should lie on the diagonal line [49,50]. Deviations from this line suggest the presence of spurious associations due to population structure and familial relatedness, and that the GWAS model does not sufficiently account for these spurious associations. It is expected that the SNPs on the upper right section of the graph deviate from the diagonal. These SNPs are most likely associated with the trait under study. By default, the QQ plots in GAPIT show only a subset of the larger p-values (i.e., less significant p-values) to reduce the file size of the graph [37,50].
MTAs for stomata size were dispersed overall on seven chromosomes, eleven MTAs at A-genomes, seventeen at B-genomes and two at D-genomes under water shortage stress conditions. Earlier researches had reported MTAs for and stated that, two SNPs situated at chromosomes 1A and 2D are linked with stomata-related attributes under water shortage environments [51]. Thirteen MTAs/QTLs controlling stomatal traits such as stomata size were located on chromosomes 3A, 4A, 6A, 6B, 2D, and 3D, accounting for 7.69 to 22.83% of the phenotypic variation under different water regimes [17,51]. In this experiment, MTAs for SS were dispersed across the three genomes under water shortage environments. Wheat scientists [33] found the significant SNPs on 5A for physiological and yield-related trait under water shortage conditions. Nine significant MTAs/QTLs for stomata frequency were detected on chromosomes 4A, 5A, 6A, 1B, and 2B, explaining from 7.65 to 30.93% of the phenotypic variation under the normal and drought stress conditions [51,52]. MTAs for SF were dispersed across 8 chromosomes with 20 SNPs at A-genomes, 19 at B-genomes and 2 at D-genomes under drought conditions (Supplementary Table S5). The significant MTAs for leaf venation were distributed across A-and B genomes, including 12 SNPs at A-genome, 24 at B-genome, and 2 at D-genome. Thomelin et al. [53] evaluated a water shortage and high temperature tolerant QTL qDHY.3BL in~1Mbp on chromosome 3B, having twenty-two responsible genes for physiological and yield indices. In this trait, 20 SNPs were reported, which showed a significant association [45]. These were situated at chromosomes 1A, 1B, 2A, 3A, 4A, 4B, 6B and 7B, and similar to our findings. Wheat scientists [54] reported that in water shortage environments, grain per spike was significantly linked with eight SNPs situated at chromosomes 1B, 2B, 4B, 6B, 2D and 5D; where in similar conditions, seven significant marker-trait associations were reported which were situated at chromosomes 5A, 1B, 4B, 5B, 6B and 2D and which supports our results for this study.
Highly phenotypic attribute associations can be described in relation to directly or indirectly influencing one attribute to other attributes. Within the genome, these attributes could be controlled by pleiotropic loci. It is proved by the presence of numerous MTAs in which one gene will show the pleiotropic influence on more associated attributes. MTA for NGS were dispersed across three genomes; nineteen at A-genomes, twenty-seven at B-genomes and four at D-genomes under water shortage conditions. Remarkably, wheat scientists identified [33,46] that chromosome 5B harbors a region regulating numerous yield-related attributes which have genomic parts linked with grain per spike. The record noteworthy and steady MTAs/QTL responsible for NGS was observed on chromosomes 1A and 2A commonly found across different environmental conditions [55]. Earlier detected MTAs for controlling NGS and GYP in wheat crop were at chromosomes 1B, 2B, 3A, 7A, 7B and 3D and 5B 1A, 4B, 6B, and 7D, respectively [18].
The marker locus on 5A was correlated with GYP in normal conditions described previously by wheat scientists [33]. Similarly, Edae et al. [46] described significant MTAs for GYP at chromosomes 4A, 1B, 5B, and 2B. Moreover, Lozada et al. [56] also specified MTAs for this trait on chromosomes 5A, 1B, 2B and 4B. MTAs for GYP on 2D, 4A and 1B have been informed by several scientists which described 27% of variations in water shortage conditions [26]. MTAs for GYP recognized in the present study were precise to different environmental conditions suggesting the dynamic nature of genetics corresponding to the grain yield per plant. The highest number of marker-trait associated genome regions and chromosomes are mentioned trait-wise in Table 5.   (Table 5). QTLs/MTAs were identified by several wheat scientists for stomatal traits on chromosomes 1A, 1B, 2B, and 7A in bread wheat genotypes under different environmental conditions [17]. Therefore, one of the important aspects in wheat breeding for increasing drought tolerance lies in a better understanding of the molecular mechanisms and genetic control of stomatal distribution and opening associated with growth rate and grain yield under abiotic stress [52]. The distal areas of chromosomes 7A and 7B are described to comprise QTL for yellow pigment content of grain, which is directed by Phytoene synthase1 (PSY-1) genes, and the occurrence of these genes can affect the pleiotropic link at 7B [60]. The TaSnRK2 gene that determined sucrose non-fermenting 1-related protein kinase was situated at chromosomes 4A and 4B. Its role is important for responding against several environmental factors, and it depicted significant association with stomata-and yield-related characteristics [25]. Wheat yield-related genes, such as TaGW2-A1 at 6A, TaTGW6-A1 at 3A, TaCwi at 2A, TaGS5-A1/TaGS-D1 at 4A, TaSus1 at 7A, and TaSus2 at 7B have been reported by wheat scientists [60].
A pleiotropic locus is correlated and alters the appearance of several phenotypic attributes. In the current experiment, many pleiotropic loci were recorded, under normal and water shortage environments. Pleiotropic locus (wsnp_Ex_c5412_9564046) on chromosome 2A at 261.05cM was significantly linked with yield and stomata-related characteristics such as FLA, SF, SS, NGS, GYP, NGS and GYP. The markers BS00038787_51 and wsnp_Ex_c24167_33416760 on the same chromosome (7A) at 369.49cM and 148.26cM, respectively, showed pleiotropic effects for observed stomata and yield indices under both environments as shown in Supplementary Table S4. In this experiment, multi-character loci for yield-and stomata-related characteristics were recognized on chromosomes 2A, 5A, 4B, 7A, 6B and 7B in both conditions as mentioned in Table 3. These regions would be a useful target for selection in breeding programs. As the interval is still large, fine genetic mapping will be necessary to demonstrate that these MTAs are a unique locus with pleiotropic effects.

Conclusions
Drought stress reduced grain yield and strong positive correlation among yield and drought-related attributes such as stomata frequency and leaf venation with grain yield per plant exists. Therefore, selection based on these characteristics would be fruitful, such as the parameter stomata size, which was negatively associated and can affect the performances of other attributes during selection procedures. In this experiment, numerous pleiotropic loci were known, and the examined characteristics such as SF, LV, NGS and GYP were also influenced by a pleiotropic locus (wsnp_Ex_c24167_33416760) on chromosome 7A at 148.26cM under normal environments. The studied characteristics governed by a pleiotropic locus (wsnp_Ex_c8913_14881924) were associated with all of the studied yieldand stomata-related traits which were located on chromosome 4B at 173.63cM under water shortage conditions. The markers BS00038787_51 and wsnp_Ex_c24167_33416760 on the same chromosome (7A) at 369.49cM and 148.26cM, respectively, showed pleiotropic effects for observed stomata and yield indices under both environments. In this experiment, multi-trait loci for yield and water shortage tolerance linked attributes was recognized on chromosomes 2A, 5A, 7A, 4B, 6B and 7B under both conditions. Supplementary Materials: The following are available online at https://www.mdpi.com/article/10 .3390/agronomy11081646/s1, Figure S1: The result obtained of 96 spring wheat genotypes using studied SNPs from Structure Harvester analysis. Figure S2: Population structure of 96 bread wheat accessions. Figures S3-S14: Quantile-quantile (QQ)-plot of p-values.