Genome-Wide Association Study Reveals Novel Powdery Mildew Resistance Loci in Bread Wheat

Powdery mildew (PM), caused by the fungal pathogen Blumeria graminis f. sp. tritici (Bgt), significantly threatens global bread wheat production. Although the use of resistant cultivars is an effective strategy for managing PM, currently available wheat cultivars lack sufficient levels of resistance. To tackle this challenge, we conducted a comprehensive genome-wide association study (GWAS) using a diverse panel of 286 bread wheat genotypes. Over three consecutive years (2020–2021, 2021–2022, and 2022–2023), these genotypes were extensively evaluated for PM severity under field conditions following inoculation with virulent Bgt isolates. The panel was previously genotyped using the Illumina 90K Infinium iSelect assay to obtain genome-wide single-nucleotide polymorphism (SNP) marker coverage. By applying FarmCPU, a multilocus mixed model, we identified a total of 113 marker–trait associations (MTAs) located on chromosomes 1A, 1B, 2B, 3A, 3B, 4A, 4B, 5A, 5B, 6B, 7A, and 7B at a significance level of p ≤ 0.001. Notably, four novel MTAs on chromosome 6B were consistently detected in 2020–2021 and 2021–2022. Furthermore, within the confidence intervals of the identified SNPs, we identified 96 candidate genes belonging to different proteins including 12 disease resistance/host–pathogen interaction-related protein families. Among these, protein kinases, leucine-rich repeats, and zinc finger proteins were of particular interest due to their potential roles in PM resistance. These identified loci can serve as targets for breeding programs aimed at developing disease-resistant wheat cultivars.


Introduction
Bread wheat (Triticum aestivum L.) is a vital source of nutrients and serves as the primary staple food crop for approximately 35% of the global population [1].In 2022, worldwide wheat consumption increased from 783 million metric tons to 792.69 million metric tons [2].This rise in wheat consumption is necessary for global food security, but has led to extra burdens in the wheat production system, which is already challenged by the prevalence of pests and diseases such as powdery mildew (PM), a devastating wheat disease.
PM is caused by the fungal pathogen Blumeria graminis (DC.)E.O.Speer f. sp.tritici Em.Marchal (syn.Erysiphe graminis f. sp.tritici) (Bgt).It is the third most destructive disease of wheat, causing yield losses ranging from 13% to 34% under low infestation and 50% Plants 2023, 12, 3864 2 of 18 to 100% under severe infestation in various wheat-growing regions [3,4].Environmental factors like temperature and relative humidity significantly influence the development and spread of this disease, contributing to epidemics [5].In the Himalayan regions of India, PM is a significant factor contributing to reduced grain yield [6]. Breeding approaches have proven effective in developing PM-resistant wheat varieties to overcome yield losses and increase economic production.Using resistant varieties is currently considered an efficient, environmentally safe, and successful strategy to control PM and mitigate production losses.
PM resistance genes have primarily been transferred from both domesticated and wild relatives and from other species such as rye (Secale cereale) to wheat [29].However, their utilization in wheat breeding is challenged by linkage drags [15].For instance, the Pm8 gene, which derived from S. cereale and also significantly contributed to PM resistance in wheat breeding during the 1990s, is linked with secalin glycopeptide locus on 1RS chromosome segment which causes a decline in flour quality [30].Furthermore, many of the identified resistance genes have been overcome by the pathogen in different countries, leaving only a few that are still effective against current Bgt isolates in the field [31,32].Thus, the continuous identification, characterization and deployment of new resistance genes is essential to minimize economic losses.
The declining cost of genotyping and the increasing availability of diverse singlenucleotide polymorphism (SNP) fingerprinting platforms has furthered the identification and deployment of new resistance genes.Many genetic mapping methods are now available to accommodate different types of traits, mapping populations, and study objectives [33].These methods can be categorized into three main groups based on the principle and genetic material used: (i) QTL mapping using a biparental population, (ii) association mapping or genome-wide association studies (GWAS) in a diverse set of germplasm, and (iii) genetic mapping utilizing the multiparent populations.The first two approaches, QTL mapping and GWAS, have successfully identified genes and their markers associated with disease resistance [34][35][36][37].However, the QTL mapping approach, which relies on a segregating population generated from parents exhibiting contrasting trait performances [12], faces limitations in QTL detection power.The scarcity of recombination events and genetic variation in such populations often leads to identifying only a few QTLs.To overcome this limitation, advanced wheat populations such as 'multiparent advanced generation intercross' (MAGIC) have been introduced to enhance the resolution of traditional mapping [38].Nevertheless, challenges arise in validating map order and analyzing recombination due to using bi-allelic markers in these populations [39].Conversely, GWAS circumvents the need for developing populations by conducting analyses on diverse germplasm [40] and offers several advantages such as the utilization of larger sample sizes, which facilitates the exploration of a broader range of genetic variations including those with minor effects and allelic diversity [36,41].Additionally, GWAS enables higher-resolution mapping, identifying candidate genes associated with the trait of interest [42].
Thus, the objectives of the present study were to (i) identify potentially novel loci that confer effective PM resistance at the adult plant stage utilizing a panel of 286 wheat accessions known as the wheat association mapping initiative (WAMI) panel created by CIMMYT, Mexico, (ii) compare identified loci with the previously reported genes/QTLs, and (iii) identify potential candidate genes associated with PM resistance using compar-Plants 2023, 12, 3864 3 of 18 ative analyses.Leveraging landmark resources such as the high-density Illumina 90K Infinium iSelect assay [43] and the Ensemble plants database [44], we aim to address the aforementioned objectives and shed light on the genetic basis of PM resistance in wheat using a GWAS study.

Phenotypic Analysis
The phenotypic analysis conducted in specific environments showed significant variation for disease severity in the WAMI panel.Table 1 provides a summary of the phenotypic variation for PM disease severity.The broad-sense heritability for PM disease severity ranged from 65% to 89% (Table 1).The analysis of variance indicated that both genotype and environment as well as genotype-environment interactions had highly significant (p ≤ 0.001) estimated variance components for disease response to PM (Table 2).However, despite environmental variation, significant positive correlations (>0.6 for each environment) were observed for the disease severity among the three environments (Figure 1A).CIMMYT, Mexico, (ii) compare identified loci with the previously reported genes/QTLs, and (iii) identify potential candidate genes associated with PM resistance using comparative analyses.Leveraging landmark resources such as the high-density Illumina 90K Infinium iSelect assay [43] and the Ensemble plants database [44], we aim to address the aforementioned objectives and shed light on the genetic basis of PM resistance in wheat using a GWAS study.

Phenotypic Analysis
The phenotypic analysis conducted in specific environments showed significant variation for disease severity in the WAMI panel.Table 1 provides a summary of the phenotypic variation for PM disease severity.The broad-sense heritability for PM disease severity ranged from 65% to 89% (Table 1).The analysis of variance indicated that both genotype and environment as well as genotype-environment interactions had highly significant (p ≤ 0.001) estimated variance components for disease response to PM (Table 2).However, despite environmental variation, significant positive correlations (>0.6 for each environment) were observed for the disease severity among the three environments (Figure 1A).[45].The red curve line is the smoothing spline regression model fitted to LD decay.The horizontal yellow line is the standard critical r 2 value of the genome (r 2 = 0.3), and the vertical blue line is the genetic distance (1.5 centimorgan, cM) at the intersect between the standard critical and the LD decay curve.The vertical green line is the genetic distance (3.0 cM) at which the LD half-decay (r 2 = 0.23, the horizontal purple line) intersects with the LD decay curve.

Population Structure Analysis
Principal components analysis revealed largely even distribution for the allele frequencies within the WAMI panel (Figure 1B).PC1 accounted for 43.3% of the total variance, while PC2 explained 18.8%.These results suggest that the combination of PC1 and PC2 captures a significant portion of the underlying genetic variation within the panel, allowing for identifying and differentiating the subgroups present in the WAMI population.

Genome-Wide Marker-Trait Association (MTA) Analysis
A scatter plot for r 2 values of pairwise markers showing genome-wide linkage disequilibrium (LD) decay for 286 genotypes of the WAMI panel assessed using the Hill and Weir formula [45].LD analysis was performed on data generated from TASSEL using a 100-SNP sliding window.The average genome-wide r 2 was 0.11, and LD decay began at r 2 = 0.46 and reached half-decay at r 2 = 0.23 (Figure 1C).The LD decay curve intersected with the half-decay and standard critical (r 2 = 0.3) lines at 3.0 and 1.5 centimorgans (cM), respectively.This defines 1.5 cM as the genome-wide critical distance to detect linkage as described earlier [46].
The comprehensive genome-wide MTA analysis encompassed three distinct models: the general linear model (GLM), mixed linear model (MLM), and the fixed and random model circulating probability unification (FarmCPU); however, the final GWAS analysis exclusively incorporated the best-fitting FarmCPU model.Over the course of the threeyear investigation, a total of 19, 85, and 9 MTAs were successfully identified in each respective year (Figure 2 and Table 3).Figure 2 displays the frequency distribution of the best linear unbiased prediction (BLUP) values for PM disease severity for three consecutive years namely, 2020-2021, 2021-2022, and 2022-2023, the Manhattan plot illustrating the results of genome-wide association scans, and the quantile-quantile (Q-Q) plot of p-values, comparing the observed −lg (P) for PM resistance to Bgt isolates in FarmCPU method. 1 = Marker position (genetic) on consensus map in cM; 2 = marker position (physical) on chromosome in bp; 3 = MTA reported previously or found associated with a known powdery mildew resistance gene or region; $ = marker with candidate genes that encode a protein; PNM = potentially novel marker-trait association.Associated SNP markers (above threshold) that were common over the years are underlined.In the initial year (2020-2021), the distribution of MTAs spanned multiple chromosomes, with notable occurrences on the chromosome 1A (2 MTAs), 2B (4 MTAs), 6B (9 MTAs), and 7A (4 MTAs).Notably, a striking concentration of 9 MTAs emerged at a specific genetic position (59.0-60.0cM) on chromosome 6B.This intriguing observation suggests the presence of a potential genomic "hotspot" linked to PM resistance, signifying the collective impact of multiple genetic variants.

Candidate Genes
In our study, we investigated genomic regions to find candidate genes within a specific window of 200 kilobases (kb) surrounding each MTA for their association with PM resistance.This window consisted of 100 kb on each side of the MTA.We discovered a total of 94 unique candidate genes (89 encoding known protein domains) in 51 genomic intervals for 65 MTAs (Table S1), while no candidate gene were spotted in defined genomic regions of the remaining 48 MTAs.
Further analysis of these 94 candidate genes involved screening them for the presence of genes that are already known to play a role in various pathways related to the interactions between pathogens and their host organisms.The results of this screening, as presented in Table S2 revealed that 30 candidate genes encoded a diverse range of protein domains relevant to plant defense and the interactions between pathogens and hosts.The identified proteins encompassed a wide array of functional domains and families.They included the following: (i) ABC transporter-like, ATP-binding domain; (ii) Ankyrin repeat; (iii) aspartic peptidase A1 family; (iv) peroxidase; (v) cytochrome P450; (vi) disease resistance protein (NB-LRR); (vii) glutathione S-transferase; (viii) Ploop-containing nucleoside triphosphate hydrolase; (ix) kinase-like domain superfamily (including serine-threonine/tyrosine/cysteine-protein kinase); (x) wall-associated receptor kinase, galacturonan-binding domain; (xi) zinc finger proteins; (xii) F-box superfamily.These candidate genes and their corresponding proteins are directly or indirectly involved in the response of the host organism to pathogen attacks.

Discussion
To uncover novel PM resistance genes/MTAs, a series of field experiments were conducted over three years, 2020-2021, 2021-2022, and 2022-2023, utilizing the diverse WAMI panel of common wheat.The analysis of PM disease severity unveiled a wide spectrum of phenotypic variations, displaying normal distributions in two of the three environments (2020-2021 and 2022-2023), and a near binomial distribution in the remaining environment (2021-2022) (Figure 2).This distribution pattern suggests the involvement of major and multiple quantitative trait loci (QTLs)/genes that contribute to PM resistance.As indicated by the phenotypic analysis, GWAS analysis also identified a total of 113 MTAs associated with adult plant resistance against PM (Table 3).This can be explained by the extensive diversity present for PM resistance within the WAMI panel [63].The presence of multiple PM resistance QTLs/MTAs in wheat germplasm sourced from diverse countries has also been reported previously [16], underlining the intricate and diverse genetic foundation of PM resistance.These findings highlight the complexity and diversity of the genetic basis underlying PM resistance and reinforce the importance of the WAMI panel in dissecting PM resistance and emphasizing the need for further investigations to unravel the specific loci and PM resistance mechanisms involved in wheat.

Population Structure Analysis
Population structure is a crucial consideration in GWAS to mitigate false-positive as-sociations [64].Previously, the WAMI genotype panel has been analyzed using the same marker set employing STRUCTURE software v2.3.4 [65].Thus, in this study, to avoid redundancy, we opted for PCA analysis, which revealed that the allele frequencies were largely evenly distributed, indicating a lack of distinct subpopulations due to the absence of strong correlations among the accessions within the panel.In contrast, Singh et al. [65], employed a model-based analysis using STRUCTURE software v2.3.4,identifying six distinct subpopulations with varying genetic purity and admixture.These results could be attributed to the composition of the panel, which includes many synthetically derived lines [64], possibly corresponding to elite parents.Conversely, the divergence between these two sets of results underscores the importance of the methodologies employed.Our PCA analysis provides an overarching view of genetic diversity without highlighting subpopulations, while previous model-based analysis dissects the panel into distinct subgroups, revealing the intricate population structure.While this panel has been successfully used for GWAS for spot blotch resistance previously, we found it suitable for this study on PM resistance.

Investigating the Concordance among Significant MTAs and Previously Mapped QTL/Pm Resistance Genes
This study identified a total of 113 MTAs (19, 85 and 9 MTAs in 2020-2021, 2021-2022 and 2022-2023, respectively) with some notable findings such as common associated markers over the years (Table 3).We observed that the highest number (85) of MTAs were detected in the second year (2021-2022), possibly because the presence of multiple markers at the same position in linkage disequilibrium [66].As indicated by the slightly higher disease pressure in the second year (Figure 2), the more favorable environmental conditions might have also contributed to this outcome.
Upon comparing the position of markers associated with significant MTAs identified in this study with the previously reported MTAs/QTLs and Pm resistance genes, it was observed that several MTAs were located to the positions where Pm resistance genes/QTLs had been reported previously while others represent novel loci.For instance, on chromosome 1A, we identified an MTA at the 27.0 cM position during the 2020-2021 season.This MTA shares a location with previously known Pm resistance loci, including resistance genes like Pm3 [47] and Pm223389 [48], as well as QTLs such as QPm.mgb-1AS [20], QPm.osu-1A [49], QPm.caas-1A [50], and a meta-QTL MQTL1 [29].Additionally, on the same chromosome, during the same 2020-2021 season, we found another MTA at the 74.0 cM position that is co-located with five MTAs previously identified by Liu et al. [12], as well as SNP markers like tplb0041a22_935, Excalibur_c15098_591 [17], and BS00021714_51 [46].
The presence of two MTAs at 68.0 cM on chromosome 5B (detected during 2021-2022) may be associated with a previously reported broad-spectrum PM resistance gene Pm16 [55], highlighting its potential for marker-assisted selection in breeding programs targeting PM resistance.
Among the chromosomes examined, chromosome 6B stood out with the largest number of MTAs, totaling 57, one at 0.0 cM while remaining within the 59.0 to 65.0 cM region.The PM resistance MTA at 0.0 cM appears to overlap with a previously identified QTL on chromosome 6B from the Swedish winter wheat cultivar Folke [53].However, no previously identified loci were found in second region, suggesting the presence of novel genes or regulatory elements associated with PM resistance in this region.
Chromosome 7B exhibited four MTAs, with one MTA detected at 73.0 cM during 2022-2023 and the other three detected at 136.0 cM during 2021-2022.A QTL QPm.mgb-7BL [20] has been previously reported in vicinity of the MTA detected at 73.0 cM position on chromosome 7B.Previously, Simeone et al. [20] revealed that QPm.mgb-7BL co-localizes with Pm resistance genes PmE [60], Pm40 [61] and Pm47 [62].It appears that markers detected at 136.0 cM position represent a novel PM resistance locus since no resistance loci has been observed in this region previously.
In addition to the above-discussed MTAs, a number of other novel MTAs were detected on chromosomes 3A (at 158.0 cM, detected during 2022-2023), 4B (at 66.0 cM, detected during 2022-2023), and 5A (at 144.0 cM, detected during 2021-2022).No QTLs/MTAs or genes were previously reported by other studies in their vicinity.These findings emphasize the importance of further investigations focused on these and other novel MTA regions such as those on chromosome 6B (between 59.0 and 65.0 cM, detected during 2020-2021 and 2021-2022) for a deeper understanding of the genetic mechanisms underlying PM resistance.
The co-localization of MTAs with previously identified genetic markers, MTAs, QTLs and genes provides further evidence for their involvement in the resistance against PM.Identified novel MTAs could be valuable for diversifying the sources of genetic resistance in breeding programs and provide valuable targets to enhance PM resistance by combining with repetitively detected loci using the associated SNP markers identified in this study.Moreover, these findings provide valuable information for breeders to prioritize genomic regions and genetic markers associated with PM resistance and contribute to a better understanding of the genetic basis of PM resistance in common wheat.

Candidate Genes for PM Resistance
To identify potential target genes for breeding, GWAS is often used as a starting point [69].However, a challenge arises when significant markers associated with the trait of interest encompass a wide range of genes within their confidence intervals, making it difficult to pinpoint the exact causal genes [36].Thus, in this study, we addressed this issue by focusing on MTAs found associated with genomic intervals carrying disease resistance-associated candidate genes.A total of 30 out of the 96 candidate genes broadly represent twelve plant protein families that exhibited clear connections to disease resistance and host-pathogen interactions.These proteins have been previously reported to impart disease resistance to various pathogens such as Blumeria graminis f. sp.tritici, Fusarium graminearum, Puccinia striiformis, Bipolaris sorokiniana, Puccinia triticina, Parastagonospora nodorum and Pyrenophora tritici-repentis in plants [70][71][72][73][74][75][76][77][78][79].
It may be recalled that 30 candidate genes encode proteins with distinctive domains, including the ABC transporter-like domain, ATP-binding domain, ankyrin repeat, aspartic peptidase A1 family, peroxidase, cytochrome P450, disease resistance protein (NB-LRR), glutathione S-transferase, NB-ARC P-loop-containing nucleoside triphosphate hydrolase, kinase-like domain superfamily, wall-associated receptor kinase, galacturonan-binding domain, zinc finger proteins, and F-box superfamily domain (Table S2), which are widely recognized as essential components of disease resistance mechanisms.These results reinforce findings from previous research.For instance, Peng and Yang's [80] comprehensive analysis delved into ABC, NLR, and START genes in hexaploid wheat, revealing their colocalization with disease resistance quantitative trait loci (QTLs) associated with leaf rust resistance.Molecular characterization of the leaf rust-resistance gene Lr34, which encodes an ABC transporter with transmembrane (TM) and nucleotide binding site (NB) domains, was elucidated in studies by Krattinger et al. [81,82].Similarly, Lr14a, a race-specific leaf rust resistance gene, was found to encode an ANKTM protein [83].The YrU1 protein, conferring wheat stripe rust resistance, also featured an integrated ANK domain from ANKTM proteins, suggesting their role in disease resistance [83,84].The up-regulation of A1 aspartic peptidase and G1 families during Zymoseptoria tritici infection [72] underscores their defensive roles.Similarly, peroxidases' positive contribution against Puccinia striiformis infection [85] highlights their significance in immunity.Disease sensitivity genes like Tsn1 which encode a protein with S/TPK and NBS-LRR domains [71,74] reveal their significance in recognizing necrotrophic effectors.The critical role of the NB-ARC domain in regulating R protein activity [70] deepens our understanding of pathogen recognition.Wall-associated receptor kinases (WAKs), exemplified by Stb6 [73] and TaWAK6 [76], emerge as potent defenders against Septoria tritici blotch and leaf rust.During defense, identified proteins exhibit versatile functions [75,79], encompassing stress resilience, apoptosis, transcription, and interactions.The F-box family protein, known for its role in numerous biological processes, including biotic stress resistance [86], highlights the multifaceted engagement of these proteins in wheat's resistance responses.However, despite the identification of several strong SNPs associated with PM resistance, no candidate genes were found in genomic intervals of 48 MTAs which may potentially carry other novel structural variants responsible for regulating plant defense against Bgt and wheat-Bgt interaction; however, further investigations are required to determine the significance and potential functional relevance of these specific SNPs in Bgt resistance.
In conclusion, our findings unveiled the intricate genetic landscape of PM resistance within the WAMI panel.This collection showcases a valuable pool of significant PM resistance MTAs/genes embedded in elite genetic backgrounds.Finally, through a comparison with previous studies on disease resistance in wheat, we have investigated the functions of key MTAs.

Plant Material
The WAMI panel, consisting of 286 genetically diverse and elite advanced wheat lines, was utilized for GWAS of PM resistance.The panel was assembled and distributed through the International Wheat Improvement Network (IWIN) by the International Maize and Wheat Improvement Center (CIMMYT) and possesses a narrow range of variation for days to heading and plant height which is appropriate for gene discovery without the confounding effects of phenology and plant height [63].The seed of the WAMI panel was obtained from CIMMYT.A highly PM-susceptible common wheat cultivar WL711 [87] was used as a susceptible check and disease spreader.

Experimental Design and Trait Evaluation
The field evaluations of genotypes were conducted at the research farm (30.75 • N, 77.30 • E) of Eternal University, Baru Sahib, Himachal Pradesh, India, spanning three consecutive growing seasons (2020-2021, 2021-2022, and 2022-2023).The Baru Sahib PM disease nursery research farm is situated in a north Indian hilly state, Himachal Pradesh, and is known for climatic conditions favorable for the natural infection, growth, and the development of PM [88].During the PM disease development period, the temperature and relative humidity in Baru Sahib range between 15 and 22 • C (max. 30 • C) and 70 and 90%, respectively.The region is a natural hot spot for PM occurrence evidenced by data from the past 15 years (Dr.H. S. Dhaliwal, personal communication).
The experiments adhered to a randomized block design and were replicated twice.Random assignment of lines to each replication was carried out using the Fisher and Yates Random Table method [89].Fertilizers included 120 kg N, 60 kg P 2 O 5 , and 40 kg K 2 O. Except N, fertilizers were thoroughly applied at the time of sowing.The application of nitrogen was divided into three doses: half at sowing, one-fourth at the first irrigation (21 days after sowing), and the remaining one-fourth at the second irrigation (45 days after sowing).
A mixed population of the Bgt isolates was used in their natural state to create epiphytotic conditions in disease nurseries, which enhanced the severity of the disease.For this purpose, two rows of WL711, a highly PM-susceptible cultivar, were planted around the field, with an additional row seeded between the plots one month before sowing the WAMI panel in November each year.The WL711 was also sown alongside the experimental materials during main season sowing.Planting distances were maintained at 20 cm between rows and 5 cm between plants.
Disease severity was scored using a rating scale of 1 to 9, following the method outlined by Bennett and Westcott [90], with some modifications on randomly chosen five plants from each plot.At this scale, level 1 denotes a severity of less than 1%, while level 2 corresponds to a range of 1-5%.The severity increases gradually through the scale, with level 9 indicating a severity of over 90%.Additionally, there's a distinct note that severity scores falling within the range of 1 to 4 are denoted as resistance (Figure S1).The entire panel was rated once the powdery mildew reached its maximum expression, corresponding to a score of 9, observed on the susceptible check cultivar WL711.This approach enabled the tracking of disease progression over time and the identification of an appropriate growth stage exhibiting varying degrees of resistance or susceptibility to PM.

DNA Extraction and Genotyping
DNA extraction, genotyping of samples and data processing were performed as described previously [91].Briefly, DNA was extracted from fresh leaves of each line using the CTAB procedure, as outlined by Saghai-maroof et al. [92].Subsequently, genotyping was conducted at the USDA-ARS Small Grain Genotyping Center, Fargo, utilizing the Illumina 90K Infinium iSelect assay (Illumina Inc., San Diego, CA, USA) [43].The SNP calling process utilized the default clustering algorithm integrated into Genome Studio v2011.1 (Illumina Inc., San Diego, CA, USA), resulting in the identification of a total of 26,814 bi-allelic SNPs [91,93].To uphold data quality standards, SNPs characterized by a minor allele frequency (MAF) lower than 0.05 were omitted from the analysis, alongside monomorphic and low-quality SNPs.This meticulous filtration procedure led to the retention of approximately 21,132 polymorphic SNPs [91,94], which were subsequently harnessed for the purpose of conducting the GWAS in this study.

Descriptive Statistics Analyses
Pearson's correlation coefficients analysis was performed using the Agricolae (version 1.2-4) package of the R (version 4.0.3)software [95].Variance components were estimated using restricted maximum likelihood (REML) method implemented in META-R software v6.0.4 [96].Covariance parameter estimates were used as unbiased genetic variance component estimates for calculating broad-sense heritability (H 2 ) following Alemu et al. [46].

Population Structure, Kinship Matrix and Principal Components Analyses
Population structure (population structure matrix or Q matrix) were modeled using principal components analysis (PCA) utilizing the genotypic data belonging to a total of 21,132 high-quality SNPs as described earlier [98].The relatedness matrix (Kinship or K matrix) was calculated with R (version 4.0.3)software [95] employing the parameters suggested by VanRaden [99] and Yin et al. [100].The optimal numbers of PCAs were determined using the Bayesian information criterion (BIC) [101].The first two principal components were used to create a scatter plot that visualized the distribution of genotypes.

Genome-Wide Association Analyses
Genome-wide association analyses were conducted using a dataset of 21,132 highquality SNPs available from the CIMMYT, Mexico website (https://data.cimmyt.org/dataset.xhtml?persistentId=hdl:11529/10714; accessed on 24 July 2023) [91].Pairwise squared allele-frequency correlations (r 2 ) between SNP markers were calculated using the TASSEL (Trait Analysis by Association, Evolution, and Linkage) software version 5.2.91 with a sliding window size of 100.To assess LD between loci, r 2 values were plotted against genetic distance in cM.The LD decay curve was fitted using a smoothing spline regression line at the genome level, following the method outlined by Hill and Weir [45] and implemented in the R environment with a script previously employed by Marroni et al. [102].
Single-and multilocus models including the general linear model (GLM) [98], mixed linear model (MLM) [103], and the fixed and random model circulating probability unification (FarmCPU) [104] method were fitted to identify marker-trait associations (MTAs).Among these, FarmCPU was selected based on the best fit and used for the final GWAS.In FarmCPU, both fixed and random effects are incorporated to improve the accuracy of association mapping.It is a two-step procedure where the fixed effect model is initially fitted to control for population structure, and then the random effect model is used to further account for relatedness among individuals.The rMVP (R based Memory-efficient, Visualization-enhanced, and Parallel-accelerated Tool For Genome-Wide Association Study) software version 1.0.0 [100] was employed for the analysis.

Identification of Candidate Genes
To identify potential candidate genes, we focused on the most significant MTAs.The marker sequences associated with these MTAs were aligned with the wheat genome assembly IWGSC v.1 downloaded from the Ensembl database.Specifically, we analyzed a 200 kb window around each MTA marker (100 kb on each side) to extract highly significant and annotated candidate genes following Singh et al. [106] and Feng et al. [107].A 200 kb window is a reasonable size to capture regulatory hubs that may be involved in the expression of a candidate genes.It is also a large enough window to capture multiple genes, which can be helpful for identifying co-expressed genes that may be involved in the same biological pathway.For the gene ontology (GO) annotation of these candidate genes, we utilized IWGSC (http://www.wheatgenome.org,accessed on 5 August 2023), which offered the requisite information for annotating the GO terms associated with the identified candidate genes.

Conclusions
Incorporating the insights gleaned from the three-year investigation, a focused set of significant MTAs emerged, each holding the potential to drive advancements in resistance to PM disease in wheat.Notably, recurrent novel MTAs, consistently observed on chromosome 6B across 2020-2021 and 2021-2022 (Table 3), underscore the enduring influence of key genomic regions on PM resistance.By harnessing the knowledge embedded within these MTAs, breeders can refine their selection processes, develop tailored markers, and ultimately expedite the development of superior PM-resistant cultivars through precision-guided approaches.This study not only sheds light on the genetic architecture underpinning PM resistance but also paves the way for innovative and accelerated wheat breeding programs that align with the demands of sustainable and resilient agricultural practices.Moreover, we identified potential candidate genes involved in disease resistance mechanisms in wheat by analyzing significant SNPs and their associated gene sequences.

Supplementary Materials:
The following supporting information can be downloaded at: https://www.mdpi.com/article/10.3390/plants12223864/s1,Supplementary File 1.xlsx.Table S1: Details of candidate genes identified in genomic intervals of marker-trait associations detected during experimental years 2020-2021, 2021-2022 and 2022-2023.Table S2: Details of defense and host-pathogen-associated proteins encoded by candidate genes identified in genomic intervals of marker-trait associations detected during this study.Supplementary File 2.docx.Figure S1.Disease severity assessment scale.

Figure 1 .
Figure 1.Correlation coefficient, population structure and linkage disequilibrium (LD).(A).A graphical representation of correlation coefficient among three environmental years for powdery mildew disease severity.(B).Scatterplots showing the results of principal component analysis (PCA) conducted on genotypic data obtained from the WAMI panel containing 286 wheat accessions.The figure highlights the population structure of the WAMI panel, as revealed by the first two principal components (PC1 and PC2).(C).A scatter plot for r 2 values of pairwise markers showing genome-wide linkage disequilibrium decay in 286 genotypes of the WAMI panel assessed using the Hill and Weir formula[45].The red curve line is the smoothing spline regression model fitted to LD decay.The horizontal yellow line is the standard critical r 2 value of the genome (r 2 = 0.3), and the vertical blue line is the genetic distance (1.5 centimorgan, cM) at the intersect between the standard critical and the LD decay curve.The vertical green line is the genetic distance (3.0 cM) at which the LD half-decay (r 2 = 0.23, the horizontal purple line) intersects with the LD decay curve.

Figure 1 .
Figure 1.Correlation coefficient, population structure and linkage disequilibrium (LD).(A).A graphical representation of correlation coefficient among three environmental years for powdery mildew disease severity.(B).Scatterplots showing the results of principal component analysis (PCA) conducted on genotypic data obtained from the WAMI panel containing 286 wheat accessions.The figure highlights the population structure of the WAMI panel, as revealed by the first two principal components (PC1 and PC2).(C).A scatter plot for r 2 values of pairwise markers showing genomewide linkage disequilibrium decay in 286 genotypes of the WAMI panel assessed using the Hill andWeir formula[45].The red curve line is the smoothing spline regression model fitted to LD decay.The horizontal yellow line is the standard critical r 2 value of the genome (r 2 = 0.3), and the vertical blue line is the genetic distance (1.5 centimorgan, cM) at the intersect between the standard critical and the LD decay curve.The vertical green line is the genetic distance (3.0 cM) at which the LD half-decay (r 2 = 0.23, the horizontal purple line) intersects with the LD decay curve.

Figure 2 .
Figure 2. PM disease severity distribution, marker-trait associations and quantile-quantile (Q-Q) plot distributions.Histograms on the left (A,D,G) display the frequency distribution of BLUP estimated for PM disease severity for three consecutive years, 2020-2021, 2021-2022, and 2022-2023.Manhattan plots in the middle (B,E,H) illustrate the results of genome-wide association scans for PM resistance using BLUP of the 286 lines of the WAMI panel across 2020-2021, 2021-2022, and 2022-2023, respectively.The horizontal orange line represents FDR-adjusted p < 0.001.Associated SNP markers (above threshold) that are common among years are highlighted by star signs on the Manhattan plots.Colored bars under the Manhattan plot represent the density of markers.The Q-Q plots of p-values on the right (C,F,I) compare the observed −lg (P) for PM resistance to Bgt isolates in FarmCPU analyses to the expected values of −lg (P) under a uniform distribution.

Figure 2 .
Figure 2. PM disease severity distribution, marker-trait associations and quantile-quantile (Q-Q) plot distributions.Histograms on the left (A,D,G) display the frequency distribution of BLUP estimated for PM disease severity for three consecutive years, 2020-2021, 2021-2022, and 2022-2023.Manhattan plots in the middle (B,E,H) illustrate the results of genome-wide association scans for PM resistance using BLUP of the 286 lines of the WAMI panel across 2020-2021, 2021-2022, and 2022-2023, respectively.The horizontal orange line represents FDR-adjusted p < 0.001.Associated SNP markers (above threshold) that are common among years are highlighted by star signs on the Manhattan plots.Colored bars under the Manhattan plot represent the density of markers.The Q-Q plots of p-values on the right (C,F,I) compare the observed −lg (P) for PM resistance to Bgt isolates in FarmCPU analyses to the expected values of −lg (P) under a uniform distribution.

Table 1 .
The descriptive statistics and heritability for the PM disease severity in the wheat association mapping initiative (WAMI) panel containing 286 wheat accessions screened at Baru Sahib PM disease nurseries in Himachal Pradesh, India during three field seasons, spanning from 2020 to 2023.

Table 1 .
The descriptive statistics and heritability for the PM disease severity in the wheat association mapping initiative (WAMI) panel containing 286 wheat accessions screened at Baru Sahib PM disease nurseries in Himachal Pradesh, India during three field seasons, spanning from 2020 to 2023.

Table 2 .
Combined analysis of variance (ANOVA) table for PM disease severity in the wheat association mapping initiative (WAMI) panel containing 286 wheat accessions.