Genetic Dissection of Drought Tolerance of Elite Bread Wheat (Triticum aestivum L.) Genotypes Using Genome Wide Association Study in Morocco

Drought is one of the most important yield-limiting factors in Morocco. Identification and deployment of drought-tolerant wheat varieties are important to cope with the challenge of terminal moisture stress and increase wheat productivity. A panel composed of 200 elite spring bread wheat genotypes was phenotyped for yield and agronomic traits for 2 years (2020 and 2021) in Morocco under rainfed and irrigated environments. The panel was genotyped using 20K SNPs and, after filtration, a total of 15,735 SNP markers were used for a genome-wide association study (GWAS) using a mixed linear model (MLM) to identify marker-trait associations (MTA) and putative genes associated with grain yield and yield-related traits under rainfed and irrigated conditions. Significant differences were observed among the elite genotypes for grain yield and yield-related traits. Grain yield performance ranged from 0.97 to 6.16 t/ha under rainfed conditions at Sidi Al-Aidi station and from 3.31 to 9.38 t/h under irrigated conditions at Sidi Al-Aidi station, while Grain yield at Merchouch station ranged from 2.32 to 6.16 t/h under rainfed condition. A total of 159 MTAs (p < 0.001) and 46 genes were discovered, with 67 MTAs recorded under rainfed conditions and 37 MTAs recorded under irrigated conditions at the Sidi Al-Aidi station, while 55 MTAs were recorded under rainfed conditions at Merchouch station. The marker ‘BobWhite_c2988_493’ on chromosome 2B was significantly correlated with grain yield under rainfed conditions. Under irrigated conditions, the marker ‘AX-94653560’ on chromosome 2D was significantly correlated with grain yield at Sidi Al-Aidi station. The maker ‘RAC875_c17918_321’ located on chromosome 4A, associated with grain yield was linked with the gene TraesCS4A02G322700, which encodes for F-box domain-containing protein. The markers and candidate genes discovered in this study should be further validated for their potential use in marker-assisted selection to generate high-yielding wheat genotypes with drought tolerance.


Introduction
Drought is the most significant abiotic factor impacting wheat yield and production, particularly in West and South Asia, North Africa, and Sub-Saharan Africa [1]. Climate change reduces precipitation, worsens drought stress, and has a significant influence on wheat production. Drought stress affects wheat grain production and quality at all phases of development and reduces wheat grain yield by around 30% [2,3]. Understanding drought-tolerance processes and identifying loci relevant to drought tolerance are essential for developing breeding tactics that promote drought tolerance [4]. In season 2022/2021, Morocco harvested 5.06 million tons of bread wheat on 3.2 million hectares, with an average of 1.5 tons per hectare.
Using molecular markers such as single nucleotide polymorphism (SNP) and diversity array technologies (DArT), marker-assisted selection (MAS) has become an essential  Table 1. Grain yield ranged from 3.31 to 9.38 t/h with a mean of 5.93 t/h and the heritability value of 0.45 for irrigated conditions at Sidi Al-Aidi station and grain yield for a rainfed experiment at Sidi Al-Aidi station ranged from 0.97 to 6.16 t/h with a mean of 4.19 t/h and the heritability value of 0.32. While grain yield at Merchouch station ranged from 2.32 to 6.85 t/h with a mean of 4.09 and a heritability value of 0.42. Irrigated conditions at Sidi Al-Aidi station recorded the highest value of grain yield followed by rainfed conditions at Merchouch and Sidi Al-Aidi stations, respectively. All traits showed significant results, wherein the irrigated conditions at Sidi Al-Aidi the NTP value had the highest coefficient of variance (CV) with a value of 23.28% followed by biomass and NSS with a value of 22.68 and 20.59%, respectively. While days to maturity (DMA) recorded the lowest value of CV (2.41%). In the rainfed experiment at Sidi Al-Aidi station, the number of seeds per spike (NSS) had the highest CV with a value of 22.85% followed by biomass and grain yield (GY) with a value of 20.39 and 15.04%, respectively. While canopy temperature (CT) recorded the lowest value of CV (3.57%). In the rainfed experiment at Merchouch station, the GY had the highest CV with a value of 18.29% followed by biomass and number of plants in m 2 (NP/m 2 ) with a value of 16.85 and 14.39%, respectively. While DMA recorded the lowest value of CV (2.97%). Traits in irrigated conditions recorded a higher heritability than in rainfed conditions. Under irrigated conditions at Sidi Al-Aidi station, NSS recorded the highest value of heritability with 0.80 followed by biomass and plant height (PLH) with 0.77 and 0.72, respectively. While GY traits showed the lowest value of heritability. Under rainfed conditions at Sidi Al-Aidi station, PLH recorded the highest value of heritability with 0.65 followed by NSS and days to heading (DHE) with 0.63 and 0.62, respectively. While the GY trait showed the lowest value of heritability. Under rainfed conditions at Merchouch station, Biomass recorded the highest value of heritability with 0.70 followed by PLH and DHE with 0.67 and 0.64, respectively. While the GY trait showed the lowest value of heritability. The top 10 elite genotypes for grain yield are indicated in Supplementary Table S1. The genotypes G-146 (SUDAN#3/SHUHA-6/4/BOW/PRL//BUC/3/WH576) and G-145 (STAR *3/LO-TUS 5//TNMU/MILAN/3/QAMAR-2) were recorded the highest grain yield under the rainfed condition with a value of 6.16 and 5.75t/h, respectively.  Supplementary Table S1. The genotypes G-146 (SUDAN#3/SHUHA-6/4/BOW/PRL//BUC/3/WH576) and G-145 (STAR *3/LOTUS 5//TNMU/MILAN/3/QAMAR-2) were recorded the highest grain yield under the rainfed condition with a value of 6.16 and 5.75t/h, respectively.

Principal Components Analysis (PCA) and Correlation
Principal component bi-plot analyses show the association between grain yield and yield-related traits ( Figure 2). Under the rainfed regime at Sidi Al-Aidi station, the two principal components explained 20.5% and 13% of the total variation for grain yield and yield-related traits. Under the irrigated regime at Sidi Al-Aidi station, the two principal components explained 17% and 15.5% of the total variation for grain yield and yield-related traits. Under the rainfed regime at Merchouch station, the two principal components explained 22.3% and 15.8% of the total variation for grain yield and yield-related traits. Figure 3 shows four groups under rainfed and irrigated regimes at both stations. Figure 4 shows the correlation between agronomic traits under rainfed and irrigated regimes at both stations. Under the rainfed regime at Sidi Al-Aidi station, DHE was positively correlated with days to maturity with a value of 0.53 and biomass was positive with PLH and biomass with a value of 0.31. Whereas, canopy temperature showed a significant negative correlation with grain yield and biomass with values of −0.21 and −0.20, respectively. Under the irrigated regime at Sidi Al-Aidi station, days to heading and days to maturity were correlated with a value of 0.65, grain yield is correlated with biomass with a value of 0.49; also, grain yield was positively correlated with plant height with a value of 0.32. Whereas canopy temperature was negatively correlated with NPM and NSS with a value of −0.18 and −0.15, respectively, and also, grain yield was negatively correlated with days to heading with a value of −0.18.     Figure 5 shows the LD decay of A, B, and D genomes. Genome B had the highest linkage disequilibrium, with an average r 2 value of 0.18, followed by genomes A and D with r 2 = 0.17 and r 2 = 0.15, respectively. Figure 6 shows the density of markers over 21 chromosomes of 200 bread wheat genotypes. For all traits under the rainfed regime at Sidi Al-Aidi station, a total of 67 significant marker-trait associations (MTAs) with p < 0.001 were found (Table 2), the A sub-genome has the most MTAs (28), followed by the B and D sub-genomes, which had 26 and 13 MTAs, respectively. For all traits under the irrigated regime at Sidi Al-Aidi station, a total of 37 significant marker-trait associations (MTAs) with p < 0.001 were found (Table 3), the A sub-genome has the most MTAs (16), followed by the D and B sub-genomes, which had 12 and 9 MTAs. For all traits under the rainfed conditions at Merchouch station, a total of 55 significant marker-trait associations (MTAs) with p < 0.001 were found (Table 3), the B sub-genome has the most MTAs (21), followed by the A and B sub-genomes, which had 20 and 13 MTAs. The GWAS results is the mean of two years (2020/2021).  Figure 5 shows the LD decay of A, B, and D genomes. Genome B had the highest linkage disequilibrium, with an average r 2 value of 0.18, followed by genomes A and D with r 2 = 0.17 and r 2 = 0.15, respectively. Figure 6 shows the density of markers over 21 chromosomes of 200 bread wheat genotypes. For all traits under the rainfed regime at Sidi Al-Aidi station, a total of 67 significant marker-trait associations (MTAs) with p < 0.001 were found (Table 2), the A sub-genome has the most MTAs (28), followed by the B and D sub-genomes, which had 26 and 13 MTAs, respectively. For all traits under the irrigated regime at Sidi Al-Aidi station, a total of 37 significant marker-trait associations (MTAs) with p < 0.001 were found (Table 3), the A sub-genome has the most MTAs (16), followed by the D and B sub-genomes, which had 12 and 9 MTAs. For all traits under the rainfed conditions at Merchouch station, a total of 55 significant marker-trait associations (MTAs) with p < 0.001 were found (Table 3), the B sub-genome has the most MTAs (21), followed by the A and B sub-genomes, which had 20 and 13 MTAs. The GWAS results is the mean of two years (2020/2021).        The results of significant marker-trait associations for grain yield and yield-related traits under rainfed and irrigated conditions at Sidi Al-Aidi and Merchouch stations are shown in Tables 2-4, and Figure 7 as well as Supplementary Figures S1-S3. For rainfed conditions at Sidi Al-Aidi station, seventy-four MTAs were recorded (Table 2), where 11 MTAs were recorded for biomass, these markers are located on chromosomes 2B, 2D, 3A, 3B, 5D, and 6D, on the same trait the markers "BobWhite_c2988_493" and "Ex_c14755_1362" located on chromosomes 2B and 2D, respectively, with a value of −Log10(p) = 4.17. Twelve markers were significantly associated with canopy temperature located on chromosomes 1B, 2B, 3D, 5B, 6A, 6B, 6D, and 7A, where the marker "AX-94463982" recorded the highest −Log10(p) = 7.03 located on chromosome 7A, whereas the marker "Bob-White_rep_c64068_241" recorded the lowest −Log10(p) = 3.09 located on chromosome 2B. Chlorophyll content recorded three significant markers located on chromosomes 3A, 5A, and 1B. Days to heading recorded two significant markers located on chromosomes 1B and 2B. Days to maturity recorded two significant markers located on chromosomes 5B and 7D. Grain yield recorded three significant markers located on chromosomes 4A, 2B, and 2D. Fifteen MTAs were recorded for NPM, which are located on chromosomes 1A, 1D, 3A, 3D, and 4A, where the marker "Ra_c29200_300" recorded the highest −Log10(p) = 3.63 located on chromosome 1A. For the NSS trait, 11 MTAs were recorded located on chromosomes 1A, 1B, 1D, 5B, 6B, and 7D. The marker "Excalibur_c5329_1335" had the highest −Log10(p) = 4.85. The number of tillers per plant (NTP) trait recorded one significant marker located on chromosome 3B. Seven MTAs were recorded for plant height located on chromosomes 1B, 3B, 4A, 5A, 5B, 5D, and 6A. The marker "RAC875_c41731_321" located on chromosome 2A was recorded for canopy temperature and chlorophyll content.  For irrigated conditions at Sidi Al-Aidi station, forty-eight MTAs (Table 3) were recorded, where 5 MTAs were recorded for biomass, these markers are located on chromosomes 5B and 7A. Chlorophyll content recorded seven MTAs located on chromosomes 1D, 2A, 3A, 3B, 5B, and 7A, where the marker "IAAV5819" recorded the highest −Log10(p) = 3.95 located on chromosome 3B. Seven significant markers were associated with canopy temperature, where the markers were located on chromosomes 1D, 2A, 3A, 3B, 5B, and 7A. DHE recorded one significant marker located on chromosome 1D with −Log10(p) = 3.81. DMA recorded 3 MTAs located on chromosomes 2A and 2D. Six significant MTAs were recorded for NPM located on chromosomes 2A, 2D, 5A, 5B, and 6A, where the marker "Excalibur_rep_c69730_391" had the highest −Log10(p) = 4.87 located on chromosome 5B. NSS trait recorded 5 MTAs located on chromosomes 4A and 5D. NTP trait recorded two significant markers located on chromosome 1B. Plant height recorded one MTA located on chromosome 7A with −Log10(p) = 3.07. The marker "AX-94513795" located on chromosome 3A was recorded for chlorophyll content and biomass.
Under rainfed conditions at Sidi Al-Aidi station, out of the two highest significant MTAs found for grain yield, BobWhite_c2988_493 and Ex_c14755_1362 located on chromosomes 2B and 2D, respectively, (Figure 8) expressed the highest frequency of favorable alleles for grain yield in the population. The results showed that wheat genotypes carrying the cytosine and guanine bases at the BobWhite_c2988_493 and Ex_c14755_1362 markers combination, respectively, significantly out-yielded the other base combinations. The combination of the cytosine and thymine bases is shown to have low-yield performance and that is linked to the low cumulative number of favorable alleles of the two significant markers. Under irrigated conditions at Sidi Al-Aidi station, two of the highest significant MTAs found for grain yield, AX-94653560 and wsnp_Ex_c6748_11659366 located on chromosomes 2D and 2B, respectively, (Figure 8) expressed the highest frequency of favorable alleles for grain yield in the population. The results showed that wheat genotypes carrying the cytosine and thymine bases at the AX-94653560 and wsnp_Ex_c6748_11659366 markers combination, respectively, significantly out-yielded the other base combinations. The combination of the thymine and cytosine bases has low-yield performance, which is linked to the low-cumulative number of favorable alleles of the two significant markers. Under rainfed conditions at Sidi Al-Aidi station, out of the two highest significant MTAs found for grain yield, BobWhite_c2988_493 and Ex_c14755_1362 located on chromosomes 2B and 2D, respectively, (Figure 8) expressed the highest frequency of favorable alleles for grain yield in the population. The results showed that wheat genotypes carrying the cytosine and guanine bases at the BobWhite_c2988_493 and Ex_c14755_1362 markers combination, respectively, significantly out-yielded the other base combinations. The combination of the cytosine and thymine bases is shown to have low-yield performance and that is linked to the low cumulative number of favorable alleles of the two significant markers. Under irrigated conditions at Sidi Al-Aidi station, two of the highest significant MTAs found for grain yield, AX-94653560 and wsnp_Ex_c6748_11659366 located on chromosomes 2D and 2B, respectively, (Figure 8) expressed the highest frequency of favorable alleles for grain yield in the population. The results showed that wheat genotypes carrying the cytosine and thymine bases at the AX-94653560 and wsnp_Ex_c6748_11659366 markers combination, respectively, significantly out-yielded the other base combinations. The combination of the thymine and cytosine bases has low-yield performance, which is linked to the low-cumulative number of favorable alleles of the two significant markers.

Gene Annotation
Forty-six putative candidate genes were discovered for all significant markers related to the grain yield and yield-related traits (Supplementary Table S2), where 34 candidate genes were recorded for rainfed conditions and 12 genes were recorded for irrigated conditions. In rainfed conditions, for biomass, the highest significant marker on chromosome 2D, 'Ex_c14755_1362' was linked to the TraesCS2D02G053300 gene, which encodes the Receptorlike serine/threonine-protein kinase protein involved in serine/threonine/tyrosine kinase activity. Seven genes were identified for canopy temperature, where the highest significant marker 'AX-94463982' on chromosome 7A was linked to the TraesCS7A02G088800 gene, which encodes Pre-mRNA-splicing factor 38 involved in the mRNA splicing, via spliceosome. In irrigated conditions, the gene TraesCS5B02G032700 was linked to the marker 'AX-94826800' located on chromosome 5B which is associated with chlorophyll content, the gene encodes for Kinesin-like protein involved in microtubule binding. The gene TraesCS5B02G300100 was linked to the marker 'wsnp_Ex_c6748_11659366' located on chromosome 5B which is associated with biomass the gene encodes for Trigger_N domain-containing protein involved in peptidyl-prolyl cis-trans isomerase activity.

Phenotypic Variability for Grain Yield and Yield-Related Traits
Bread wheat is a major crop in many parts of the world, particularly in Central and West Asia and North Africa (CWANA) as well as in Sub-Saharan Africa (SSA). In the CWANA region, water shortage is becoming a severe challenge for wheat production. The creation of novel drought-tolerant bread wheat genotypes might be a long-term solution for increasing wheat yield in challenged regions. ICARDA's bread wheat breeding program targets generating genotypes that are abiotic and biotic stress-resistant and have acceptable end-use quality. Grain yield is a complicated property governed by numerous genes, the expression of which is controlled by environmental conditions. Therefore, this study was carried out to evaluate the drought-tolerance capacity of 200 bread wheat genotypes and to discover markers that are strongly related to grain yield and yield-related traits. Phenotype observations showed that wheat in rainfed conditions had lower GY and fewer DHE than wheat in irrigated conditions. The days to heading and maturity exhibited greater heritability than grain yield at all conditions, as previously reported [35,36]. As a result, grain yield heritability is often lower than that of other characteristics [37]. The correlation between days to heading and maturity was strongly positively correlated in both conditions which were consistent with previous research [38]. Biomass was positively correlated with grain yield in both conditions, but the correlation in irrigated conditions was upper than in rainfed conditions, this result is due to the ample irrigation and optimum conditions.

Marker Trait Association
GWAS has been used in several studies on bread wheat to uncover novel MTAs and QTLs linked with distinct traits in different conditions. We found a total of 104 significant markers under rainfed and irrigated conditions associated with agronomic traits across chromosomes in this study, where 44 MTAs found on the A genome, followed by the B and D genomes with 35 and 25, respectively, similar to another study that found the majority of MTAs for grain yield and yield-related traits on the A and B genomes, followed by the D genome [39]. Drought-stressed environments often have more major MTAs than non-stressed environments, which explains why drought tolerance is highly influenced by genotype by environmental factors [38]. Sixty-seven and thirty-seven MTAs at p < 0.001 were found associated with grain yield and agronomic traits under rainfed and irrigated conditions, respectively.
The highly significant marker for grain yield in rainfed conditions was identified on chromosome 2B with −Log10(p) = 3.81 at position 338 Mbp. We identified two other significant markers on chromosomes 2D and 4A. Several studies detected many MTAs/QTLs for grain yield in the same chromosomes [5,35]. In this study, we identified five MTAs for biomass in irrigated conditions, which are located on chromosomes 7A and 5B where the highest marker was 'BS00022169_51' located on chromosome 7A with −Log10(p) = 4.14 at position 691 Mbp, a previous study found the same marker on chromosome 7A [40]. On chromosome 5B, a highly significant marker associated with a number of seeds per spike under rainfed conditions was discovered with −Log10(p) = 4.85. Several MTAs/QTLs on chro-mosome 5B have been linked to NSS in dry conditions, according to a previous study [41]. Five MTAs were found under irrigated conditions located on chromosomes 4A and 5D for grain number per spike. Using 192 bread wheat genotypes planted at the same study site in Morocco, Alemu et al. [42] discovered that the marker 'wsnp_Ex_c20386_29451037' was linked to grain number per spike on chromosome 4A. Twelve MTAs regulating canopy temperature traits under rainfed conditions were found on chromosomes 1B, 2B, 3D, 5B, 6A, 6B, 6D, and 7A. MTAs for canopy temperature on chromosomes 7A and 2B were previously reported using different wheat panels cultivated in various settings [43]. Utilizing a population of 287 advanced elite lines of spring wheat, Sukumaran et al. [44] discovered the QTL RAC875_rep_c114561_587 associated with canopy temperature on chromosome 6A. The investigation was conducted using 90,000 wheat SNP created by Infinium iSelect SNP assays. The marker 'BS00037225_51' was associated with plant height on chromosome 3B with −Log10(p) = 7.1 at position 590 Mbp. Previously, numerous MTAs/QTLs associated with plant height were discovered on the same chromosome that we discovered [5]. The important markers were mapped to the bread wheat genome reference database at Ensemble Plant and UniProt in order to discover the probable genes linked with the many assessed features previously published. The findings indicated 46 genes involved in ATP binding, mRNA binding, protein serine kinase activity, protein heterodimerization activity, defense response, and oxidoreductase activity, among other biosynthetic pathways [45][46][47].

Plant Material and Experimental Conditions
This study employed a population of 200 bread wheat genotypes from the International Center for Agricultural Research in Dry Areas (ICARDA) (Supplementary Table S3). The field experiment was conducted at ICARDA's experimental stations in Sidi El-Aidi, Morocco (33 • 07 27.6 N 7 • 37 43.8 W, 406 m a.s.l.) and Marchouch, Morocco (33 • 36 26.8 N 6 • 42 43.9 W) during the 2019-2020 and 2020-2021 cropping seasons. The genotypes were grown in a 3 m 2 plot following an alpha lattice design with two replications. In Sidi EL-Aidi station the experiments were carried out in two different water regimes, rainfed and irrigated. For Marchouch station, the experiment was carried out under rainfed conditions. During the growth cycle in the irrigated experiment, the plots were watered by drip irrigation with three irrigations each week. The Sidi El-Aidi station is characterized by vertisol soil type with annual precipitation of about 300 mm over the two years. The Marchouch station is characterized by cambisol soil type with annual precipitation of about 400 mm over the two years. The stations are characterized by moderate humidity with an annual temperature ranging between 10 and 40 • C. The panel was seeded at a rate of 100 kg ha −1 during the first week of December.

Phenotyping and Statistical Analyses
The GWAS panel was phenotyped for 11 traits: grain yield (GY), days to heading (DHE), days to maturity (DMA), plant height (PLH), number of tillers per plant (NTP), number of plants in m 2 (NPM), number of seeds per spike (NSS), biomass, chlorophyll content (CC), canopy temperature (CT), and thousand kernel weight (TKW). After threshing each plot, GY was measured in kg/plot and then converted to t/ha. When 50% and 90% of the plants reached the heading and maturity phases, respectively, DHE and DMA were recorded for each plot. When the plants reached maturity, the PLH in centimeters was recorded in each plot by measuring the distance from the soil surface to the top of the spike, without measuring the length of the awns. The biomass was measured for all plots in kilograms. The leaf chlorophyll content was measured on the flag leaf using a SPAD-502Plus chlorophyll meter at the heading stage on a clear day between 9:00 and 12:00. Canopy temperature was measured at the booting stage using a thermometer. TKW was calculated in grams using a grain counter. The "RcmdrMisc" package in R was used to determine the means, maximum, minimum, and standard deviation, as well as the coefficient of variance.

Genotyping
Fresh leaf samples were used to obtain genomic DNA from 2-week-old seedlings. Before extraction, the samples were frozen in glace liquid and maintained at −80 • C according to the methodology outlined in [48]. A total of 15,735 SNP markers were retained after quality filtering (the markers were used to generate the genotypic data for the panel used).

Linkage Disequilibrium Analyses and Population Structure
Linkage disequilibrium (LD) was calculated using Tassel software, where LD values of the 200 bread wheat genotypes were investigated using the SNP markers distributed across wheat genomes. Remington and his collaborator's approach was used to conduct the LD decay study [49]. The population structure was estimated as a principal component (PC) analysis.

Genome-Wide Association Mapping and Genes Annotation
Genome-wide association analysis for grain yield and yield-related traits was carried out using Genomic Association and Prediction Integrated Tool version 3 (GAPIT 3) in the R environment [50]. A kinship K matrix and PCs (MLM: PCs+K) were used in the analysis. The important markers were shown on Manhattan plots, and markers with a −Log10(p) > 3.0 were designated as significant MTAs across the 21 chromosomes. Only significant markers (−log10p > 3.0) were traced using the CMplot program in the R environment [51]. Functional genes associated with grain yield and yield-related traits were identified using the Triticum aestivum genome annotation files downloaded from the Ensemble Plants database (http://plants.ensembl.org/Triticumaestivum/Info/Index, accessed on 5 April 2022 and the National Center for Biotechnology Information (NCBI) database. Uniport (https://www.uniprot.org, accessed on 6 April 2022) was used to identify protein function.

Conclusions
Some genotypes demonstrated good grain yield and acceptable agronomic performance in this study, and they were included in the worldwide nurseries given by ICARDA to national partners in the CWANA and SSA areas for release and future usage as crossing block parents. Using 200 spring bread wheat genotypes, we discovered 104 MTAs and 46 genes linked with various agronomic parameters under rain-fed and irrigated conditions. The markers discovered in this work will be utilized for marker-assisted selection after validation using a new set of elite genotypes.