Genome-Wide Association Study Reveals Novel Genetic Loci for Quantitative Resistance to Septoria Tritici Blotch in Wheat (Triticum aestivum L.)

Septoria tritici blotch, caused by the fungus Zymoseptoria titici, poses serious and persistent challenges to wheat cultivation in Ethiopia and worldwide. Deploying resistant cultivars is a major component of controlling septoria tritici blotch (STB). Thus, the objective of this study was to elucidate the genomic architecture of STB resistance in an association panel of 178 bread wheat genotypes. The association panel was phenotyped for STB resistance, phenology, yield, and yield-related traits in three locations for 2 years. The panel was also genotyped for single nucleotide polymorphism (SNP) markers using the genotyping-by-sequencing (GBS) method, and a total of 7,776 polymorphic SNPs were used in the subsequent analyses. Marker-trait associations were also computed using a genome association and prediction integrated tool (GAPIT). The study then found that the broad-sense heritability for STB resistance ranged from 0.58 to 0.97 and 0.72 to 0.81 at the individual and across-environment levels, respectively, indicating the presence of STB resistance alleles in the association panel. Population structure and principal component analyses detected two sub-groups with greater degrees of admixture. A linkage disequilibrium (LD) analysis in 338,125 marker pairs also detected the existence of significant (p ≤ 0.01) linkage in 27.6% of the marker pairs. Specifically, in all chromosomes, the LD between SNPs declined within 2.26–105.62 Mbp, with an overall mean of 31.44 Mbp. Furthermore, the association analysis identified 53 loci that were significantly (false discovery rate, FDR, <0.05) associated with STB resistance, further pointing to 33 putative quantitative trait loci (QTLs). Most of these shared similar chromosomes with already published Septoria resistance genes, which were distributed across chromosomes 1B, 1D, 2A, 2B, 2D, 3A,3 B, 3D, 4A, 5A, 5B, 6A, 7A, 7B, and 7D. However, five of the putative QTLs identified on chromosomes 1A, 5D, and 6B appeared to be novel. Dissecting the detected loci on IWGSC RefSeq Annotation v2.1 revealed the existence of disease resistance-associated genes in the identified QTL regions that are involved in plant defense responses. These putative QTLs explained 2.7–13.2% of the total phenotypic variation. Seven of the QTLs (R2 = 2.7–10.8%) for STB resistance also co-localized with marker-trait associations (MTAs) for agronomic traits. Overall, this analysis reported on putative QTLs for adult plant resistance to STB and some important agronomic traits. The reported and novel QTLs have been identified previously, indicating the potential to improve STB resistance by pyramiding QTLs by marker-assisted selection.

Septoria tritici blotch, caused by the fungus Zymoseptoria titici, poses serious and persistent challenges to wheat cultivation in Ethiopia and worldwide. Deploying resistant cultivars is a major component of controlling septoria tritici blotch (STB). Thus, the objective of this study was to elucidate the genomic architecture of STB resistance in an association panel of 178 bread wheat genotypes. The association panel was phenotyped for STB resistance, phenology, yield, and yield-related traits in three locations for 2 years. The panel was also genotyped for single nucleotide polymorphism (SNP) markers using the genotyping-by-sequencing (GBS) method, and a total of 7,776 polymorphic SNPs were used in the subsequent analyses. Marker-trait associations were also computed using a genome association and prediction integrated tool (GAPIT). The study then found that the broad-sense heritability for STB resistance ranged from 0.58 to 0.97 and 0.72 to 0.81 at the individual and across-environment levels, respectively, indicating the presence of STB resistance alleles in the association panel. Population structure and principal component analyses detected two sub-groups with greater degrees of admixture. A linkage disequilibrium (LD) analysis in 338,125 marker pairs also detected the existence of significant (p ≤ 0.01) linkage in 27.6% of the marker pairs. Specifically, in all chromosomes, the LD between SNPs declined within 2.26-105.62 Mbp, with an overall mean of 31.44 Mbp. Furthermore, the association analysis identified 53 loci that were significantly (false discovery rate, FDR, <0.05) associated with STB resistance, further pointing to 33 putative quantitative trait loci (QTLs). Most of these shared similar chromosomes with already published Septoria resistance genes, which were distributed across chromosomes 1B, 1D, 2A, 2B, 2D, 3A,3 B, 3D, 4A, 5A, 5B, 6A, 7A, 7B, and 7D. However, five of the putative QTLs identified on chromosomes 1A, 5D, and 6B appeared to be novel. Dissecting the detected loci on IWGSC RefSeq Annotation v2.1 revealed the existence of disease resistance-associated genes in the identified QTL regions that are involved in plant defense responses. These putative QTLs explained

INTRODUCTION
Common wheat (Triticum aestivum L.) is the most widely cultivated and the major staple food crop in the world consumed by human, providing almost 20% of the total calories and 21% of protein demand globally (Arzani and Ashraf, 2017; International Wheat Genome Sequencing Consortium (IWGSC), 2018; Ramadas et al., 2019). By 2050, the world's human population is projected to reach nine billion, and we will need to increase wheat production by 70% to feed this projected growth (FAO, 2009;Ray et al., 2013;Marcussen et al., 2014). Hence, boosting the wheat harvest is very pertinent to achieve zero-hunger by 2050.
The deployment of genetic resistance is the most durable, economical, and environmentally friendly method to manage crop diseases like STB (Ghaneie et al., 2011;Mekonnen et al., 2019;Odilbekov et al., 2019). In particular, qualitative and quantitative types of resistance to STB have been reported in wheat (Arraiano and Brown, 2006;Arraiano et al., 2009). The former refers to a condition where one or few major Stb genes provide resistance to specific Z. tritici isolates (Brown et al., 2015). Quantitative resistance, on the other hand, results from the expression of many genes with minor effects and is generally not specific to isolates. As such, quantitative resistance is the most effective, durable, and preferred method to manage rapidly evolving wheat pathogens such as Z. tritici .
The resistance-breeding method used in Ethiopia is mainly conventional, making the crop-improvement program very slow and inefficient. Nowadays, the advent and application of modern genomic tools have revolutionized crop breeding by facilitating the precise identification, mapping, and introgression of genomic regions controlling useful agronomic traits, such as resistance, into new cultivars. To account for this, a genome-wide association study (GWAS) is a powerful approach to elucidating the genomic architecture of many traits . The development of high-throughput sequencing and bioinformatics technologies (Huang et al., 2017) has also enabled GWAS to scan single nucleotide polymorphisms (SNPs) associated with desirable traits at the whole-genome scale (Rafalski, 2010).
Genome-wide association studies have been successfully applied to many crop species (Xiao et al., 2017) such as maize (Rashid et al., 2018), rice (Huang et al., 2017), wheat (Kidane et al., 2017;Long et al., 2019;Cheng et al., 2020), and sorghum (Girma et al., 2019). In particular, this study design has been used in wheat to analyze several traits such as resistance to stripe rust Yao et al., 2019;Cheng et al., 2020), stem rust (Edae et al., 2015;Kankwatsa et al., 2017), Septoria tritici blotch (Kidane et al., 2017;Odilbekov et al., 2019), drought tolerance (Mathew et al., 2019), and other phenological characteristics, plus yield and yield-related traits (Jamil et al., 2019;Wang et al., 2019;Ward et al., 2019). While Ethiopia is the largest producer of wheat in sub-Saharan Africa, little is known about the resistance Ethiopian wheat cultivars have to STB, even though it is the most important disease economically. Thus, the objectives of this study were: (1) to determine the population structure, family relatedness, and level of linkage disequilibrium of the tested bread wheat association panel; (2) to elucidate the genomic architecture of adult plant resistance to STB; (3) to identify the SNP loci underlying yield, yield-related, and phonological traits in Ethiopian cultivars that could be useful in wheat breeding programs.

Association Mapping Panel
This study used an association panel of 180 bread wheat (Triticum aestivum L.) genotypes (Supplementary Table 1), of which 167 were obtained from the International Maize and Wheat Improvement Center (CIMMYT-Mexico) and 13 were commercial cultivars grown in Ethiopia. The 167 CIMMYT genotypes included 49 from the International Bread Wheat Screening Nursery (IBWSN), 56 from the International Septoria Observation Nursery (ISEPON), 14 from the High Rain Wheat Yield Trial (HRWYT), 34 from the High Rainfall Wheat Screening Nursery (HRWSN), 5 from an adaptation trial, 6 from the National Variety Trials (NVT), and the remaining 3 genotypes were from a preliminary variety trial (PVT).

Multi-Environment Trials
Field evaluations were carried out under natural STB infestation during the 2015 and 2016 main cropping seasons across three STB hotspots: the Holetta Agricultural Research Center (HARC) (9 • 3'N/38 • 30'E), Bekoji Agricultural Research Subcenter (7 • 32'N/39 • 15'E), and Kulumsa Agricultural Research Center (KARC) (8 • 02'N/ 39 • 15'E). The experimental design was an alpha lattice with two replications, six incomplete blocks, and 30 entries per sub-block per replication. The trial was sown by hand, with each entry planted in four rows of a 2.5-m length, 20-cm spacing between rows, and 40 cm between entries. The susceptible cv. "Lakech" was planted as a spreader row along the length of the blocks to create adequate disease pressure. The spaces between the blocks and replications were 1.5 m long. A seeding rate of 150 kg ha −1 and fertilizer rates of 100 and 75 kg ha −1 of N and P 2 O 5 , respectively, were used in all the experiments. Weeding was performed by hand three times each season.

STB Evaluation
Septoria tritici blotch disease severity (SDS) was estimated visually plot-wise by considering the percentage of necrotic leaf area (NLA) on the four uppermost infected leaves of 10-20 plants (Eyal and Levy, 1987) at three growth stages, namely, heading (SDH), medium milk (SDMM), and at maturity (SDM), using a double-digit 00-99 scoring scale (Eyal and Levy, 1987). The first digit (0-9) represented blotch development in terms of plant height (for instance, 5 if the disease reached the middle (50%) of the plant height, 8 for reaching the flag leaf, and 9 for reaching the spike), while the second digit stood for the disease severity as a percentage but in terms of 0-9 (1 = 10%, 2 = 20%, and 9 = 90 %). For each stage, Septoria disease severity percentage (SDS%) was computed from the 00-99 score using the following formula as described by Sharma and Duveiller (2007): where D1 and D2 are the first and the second digits of the doubledigit scores, respectively. The SDS values range from 0 to 100, where 0 indicates complete resistance and 100 indicates complete susceptibility (Kidane et al., 2017).
In addition, the Septoria progress coefficient (SPC) developed by Eyal and Ziv (1974) was computed to indicate the position of pycnidia relative to plant height according to the following equation: where SDH (Septoria disease height) is the maximum height (cm) above ground at which the pycnidia of the pathogen could be found on the plant at the maturity stage and PH is the average height of the genotype from the ground to the tip of its awn. The SPC coefficient indicates the position of pycnidia relative to plant height, regardless of pycnidial coverage, and allows for the comparison of the infection placements on cultivars with different plant statures. Furthermore, SPC values ranged from 0 to 1, where SPC = 0 means that there was no disease, while SPC =1 means that pycnidia were produced at the top of the plant (Eyal and Levy, 1987).

Other Agronomic Data Scoring
The phenotypic data that were recorded were heading date (HD, days to 50% heading), flowering date (FD = days to 50% flowering), grain-filling duration (GFD), maturity date (MD = days to 90% maturity), grain yield, hectoliter weight (HLW = kilograms per 100 liters of wheat), thousand-kernel weight (TKW = weight of 1,000 kernels, in grams), plant height (PH), number of spikelets per spike (SPS), number of kernels per spikelet (NKPS), and number of kernels per spike (NKS). These yield data were taken from the four rows of each plot and converted to kilograms per hectare (kg ha −1 ) at 12.5% moisture content using plot size as a factor. Plant height measurement was also carried out at physiological maturity from five randomly selected and tagged plants from the middle rows of each entry.

DNA Extraction and Genotyping by Sequencing
The wheat plants of the association panel were grown at the National Agricultural Biotechnology Research Center, Holetta under greenhouse conditions. The 2-week-old leaf samples were then collected into 96 deep-well sample collection plates, ovendried overnight at 50 • C, and sent to Integrated Genotyping Service and Support (IGSS) located at the Biosciences Eastern and Central Africa (BecA-ILRI) Hub in Nairobi, Kenya for highdensity genotyping by Diversity Arrays Technology sequencing (DArTseq TM technology). Furthermore, DNA extraction was carried out using the Nucleomag Plant Genomic DNA extraction kit (Macherey-Nagel GmbH & Co. KG, Duren, Germany). Afterward, extracted DNA quality and quantity were checked on a Thermo Scientific TM NanoDrop TM 2000 Spectrophotometers (Thermo Scientific TM , USA) and on 0.8% agarose gels. As a result, the extracted genomic DNA concentration ranged from 50 to 100 ng/µl. Whole-genome profiling was also carried out using the genotyping-by-sequencing (GBS) platform as described by Elshire et al. (2011). This method involved library construction following the DArTSeq complexity reduction method via the digestion of genomic DNA using ApeKI [a type II restriction endonuclease that recognizes a degenerate 5-bp sequence (GCWGC, where W is A or T)] and the ligation of barcoded adapters, which was also followed by the PCR amplification of adapter-ligated fragments. The libraries were then sequenced using single-read sequencing runs for 77 bases. The nextgeneration sequencing of the GBS library was also carried out using an Illumina HiSeq 2500 lane (Illumina, San Diego, CA, United States) following the protocol of the manufacturer.

Quality Control and SNP Calling
The technical quality of the sequencing was checked using a Sequencing Analysis Viewer. DArTSeq markers were scored using the DArTsoft14 software implemented in the KDCompute plug-in system developed by Diversity Arrays Technology (2017) (http://www.kddart.org/kdcompute.html) based on their alignment with the reference genome of the Chinese Spring Wheat RefSeq v1.0 [International Wheat Genome Sequencing Consortium (IWGSC), 2018], which was obtained from the International Wheat Genome Sequencing Consortium database (https://urgi.versailles.inra.fr/download/iwgsc/) at a minimum base identity of 90% and e-value of 5e-10. Two types of markers were scored, namely, SilicoDArT markers and SNP markers, which were both scored in a binary fashion (1/0), indicating the presence or absence of a marker in the genomic representation of each sample as described by Akbari et al. (2006). Marker quality was also maintained by removing monomorphic markers and those with lower call rates (>30% missing) and MAFs (minor allele frequencies) <5% using the ArTSoft14 software.

Phenotypic Data Analysis
We conducted an ANOVA for each location in each year using the SAS software version 9.2 (SAS Institute Inc., 2008) by considering genotype and the block as fixed and random factors, respectively. In an individual environment, the observed phenotypic response of the ith genotype in the jth replication and lth sub-block was computed using the following model: where y ijl = the observed phenotype, µ = the grand mean, g i = fixed effect of the i th genotype, γ j = effect of the j th replication, bl (j) = random effect of the l th block nested within the j th replication, and ε ijl = random error term. The ANOVA of all seasons and locations was executed by considering genotype as a fixed effect and the block, location, and year as random effects according to the following model: where Y ijklm = observed response of genotype m, replication k of block l of location j and year i; µ = grand mean; g m = fixed effect of genotype m; r ijk = effect of replication k in location j and year i; y ij = random effect of year i at location j that is ∼ normally and independently distributed (NID) (0, δ 2 y ); e j = random effect of location j that is ∼ NID (0, δ 2 e ); b ijkl = random effect of block l nested with replication k in location j and year i that is ∼ NID (0, δ 2 b ); (gy) im = random effect of the interaction between genotype m and year i that is ∼ NID (0, δ 2 gy ); (ge) jm = random effect of the interaction between genotype m and location j that is ∼ NID (0, δ 2 ge ); (ye) ij = random effect of the interaction between year i and location j that is ∼ NID (0, δ 2 ye ); (yeg) ijm = random effect of the three-way interaction of genotype m in location j and year i that is ∼ NID (0, δ 2 gey ); ε ijklm = random residual effect that is ∼ NID (0, δ 2 ε ).
The variance components were also computed. The broadsense heritability (H 2 ) within an environment was estimated for the traits from an ANOVA using the following formula: The broad-sense heritabilities across the environments were also estimated by the formula: H 2 = (δ 2 g)/(δ 2 g+δ 2 gy/y+δ 2 ge/l+δ 2 gye/yl+δ 2 ǫ/ylr) where δ 2 g is the genotypic variance, σ 2 gy is the genotype-by-year interaction variance, σ 2 ge is the genotype-by-location interaction variance, σ 2 gye is the genotype-by-year-by-location interaction variance, δ 2 e is the location variance, and l, r, and y represent the numbers of locations, replicates, and years, respectively. The percentage of heritability was categorized as low (<30%), moderate (30-60%), or high (≥60%) as described by Robinson et al. (1949). The relationship between agronomic traits was also determined by Pearson's correlation using the SAS software version 9.2 (SAS Institute Inc., 2008).

Population Structure Analysis
A population stratification of the association panel was visualized by principal component analysis (PCA) using the KDCompute plug in system version 1.0.1 (https://kdcompute.seqart.net/ kdcompute/plugins). Population admixture patterns were also determined using a Bayesian model-based clustering algorithm implemented in the STRUCTURE software v.2.3.4 (Pritchard et al., 2000). The STRUCTURE program was run with the admixture model, correlated allele frequencies, a burn-in period of 10,000, and 50,000 Markov Chain Monte Carlo (MCMC) replications after a burn in for hypothetical subpopulations K from 1 to 10 with 10 iterations. The optimum K value was predicted based on a study by Evanno et al. (2005) using STRUCTURE HARVESTER ver. 0.6.92 (Earl and von Holdt, 2012). A bar plot for the optimum K was determined using Clumpak beta version (Kopelman et al., 2015).

Genome-Wide Association Study
The association mapping of phenotypic traits with genome-wide scanned SNPs was conducted using the Genome Association and Prediction Integrated Tools (GAPIT) package (Lipka et al., 2012) in the R software (R Core Team, 2013). This GWAS was carried out for four Septoria disease traits, namely, SDSH, SDSMM, SDSM, and SPC, and some important agronomic traits such as the days to 50% heading (DH), days to 50% flowering (DF), grain filling duration (GFD), days to 90% maturity (DM), grain yield, thousand-kernel weight (TKW), and plant height (PH) in each individual environment; the study design also used the means across all environments [the best linear unbiased estimate (BLUE) values]. The analysis involved a total of 7,776 robust SNPs with a call rate of >70% and MAF of >5%. Missing SNP data were imputed using optimal impute ver. 1.0.0 in the KDcompute_plugin system based on the KNN imputation method. The marker distribution on each chromosome was determined using LD measure in R 2 ver.0.2.2 of the KDcompute_plugin. Pairwise LD measures (r 2 and P-value) between markers on each chromosome were also computed using TASSEL Ver. 5 (Bradbury et al., 2007). A genome-wide LD decay scatter plot was then produced by plotting the r 2 values against physical distance (bp) using the GAPIT software. Finally, r 2 = 0.2 was considered as a cutoff point for no LD between pairs of markers.
The GWAS was conducted using the fixed and random model circulating probability unification (FarmCPU) algorithm  implemented in the GAPIT R package (2.0) (Tang et al., 2016). The algorithm uses both fixedeffect and random-effect models iteratively to control spurious marker-trait associations due to population structure and family relatedness (Lipka et al., 2012). Furthermore, a kinship (K) matrix was computed using the method of VanRaden (2008). Principal components describing the population stratification were computed using R/GAPIT and iteratively added to the fixed part of the model. Quantile-quantile (Q-Q) plots generated from -log10 p-values were assessed visually to determine how well the model accounted for population structure and family relatedness among the study samples. Statistically significant marker-trait associations were declared using a false discovery rate (FDR)adjusted p < 0.05 as implemented in GAPIT. Furthermore, the Bonferroni correction rate at a significance threshold of p < 0.15 or -log10 (p-values) = 4.71 was also included in the analysis for comparison. Both the Q-Q and Manhattan plots were visualized using the R package qqman (Turner, 2014). The high-confidence candidate genes within the identified resistanceassociated regions were also extracted from the recently released IWGSC RefSeq Annotation v2.1 available on the URGI Seq repository (https://wheat-urgi.versailles.inrae.fr/Seq-Repository/ Annotations).

Adult Plant Responses to STB and Broad-Sense Heritability
The genotype effect was significant (p < 0.0001) for STB resistance at all the growth stages in all the test environments. Genotypic variance (σ 2 g) was the major contributor to STB resistance variability among the tested wheat genotypes. The Septoria disease severity traits also showed pseudo-normal distributions (Figure 1), indicating the quantitative nature of STB resistance in the tested wheat genotypes (Kidane et al., 2017). The analysis revealed that the STB infestation showed seasonal fluctuations, but that was still higher during the 2015 growing season (Table 1). Moreover, the disease severity showed an increasing trend from heading to the maturity stage. In each environment, mean SDS values at the heading and mid-maturity stages ranged from 18.2 to 31.2% and 21.7 to 37.6%, respectively, while the highest severity values were registered at Holetta in 2015. The mean disease severity at maturity and its vertical progress varied from 30 to 50.8% and 0.41 to 0.69, respectively, while they were the highest at Bekoji in the 2015 growing season. The lowest Septoria severity was recorded at Kulumsa in 2016. The broad-sense heritability for Septoria resistance in each environment ranged from moderate (H 2 = 0.58) to high (H 2 = 0.99) ( Table 2).
The combined ANOVA revealed that the effect of genotype, year, location, and their two-(genotype × year, genotype × location, and year × location) and three-way interactions (genotype × year × location) were significant for SDS traits except for the Septoria progress coefficient, where the effect of year was not significant ( Table 3). The analysis of pooled data revealed that genotypic variance (σ 2 g) was the highest for all the SDS parameters except SDSMM, where environmental effect was the highest ( Table 4). The broad-sense heritabilities of Septoria resistance traits showed that they were highly heritable (H 2 = 72-81%) ( Table 4) (Robinson et al., 1949). The phenotypic and genotypic coefficients of variation for SDS traits ranged from 32.4 (SPC) to 68.3% (SDSH) and 22.5 (SPC) to 41.4% (SDSM), respectively. At 5% selection intensity, the genetic advance for SDS traits ranged from 0.31 (SPC) to 35.23 (SDSM), while the magnitude of the expected genetic gains as a percent of the mean varied from 53.69% (SPC) to 102.38% for SDSH (Table 4).
Over all the environments, the average SDS of the individual wheat genotypes ranged from 5.3 to 39.8% at heading, 8.2 to 48.5% at mid-maturity, and 10.6 to 65.3% at maturity (Supplementary Table 2). The average SPC of the individual environments ranged from 0.37 to 0.79. The most resistant genotype at all the growth stages was G174, while G104 (39.8%), G76 (48.5%), and G127 (65.3%) were the most susceptible genotypes at the heading, mid-maturity, and maturity stages, respectively (Supplementary Table 2). A comparative severity analysis with the standard checkKing-bird (G40) and the mean performance of the released varieties also confirmed the presence of superior STB-resistant genotypes among the tested materials. Of the 180 tested genotypes, 56 (31%) at heading, 75 (42%) at mid-maturity, and 105 (59%) at maturity had numerically superior STB resistance compared with King-bird (Supplementary Table 2).
The top 5% best genotypes at maturity had 47.6-71% greater resistances than King-bird and 11.9-74.4 % greater resistances compared with the mean performances of the released varieties ( Table 5).
Pearson's correlation analysis of the means over all environments revealed that STB resistance traits were significantly negatively associated with important agronomic traits. Except for SPC, all the SDS traits showed non-significant and negligible negative correlations (r < −0.3) with plant height. Disease traits also showed little to negative associations with HD, FD, GFD, NKPS, and NKS. However, a significant weak negative association (−0.25 to −0.48) was observed between SDS traits and MD, grain yield, HLW, and TKW ( Table 6).

SNP Statistics
The Illumina HiSeq 2500 (Illumina, San Diego, CA, United States) sequencing failed to generate SNP data for two genotypes (9 and 95); hence, a total of 178 bread wheat genotypes were successfully DArTSeq genotyped. Initially, a total of 35, 672 SNPs were discovered, of which 31,052 (87%) were mapped to known chromosomal positions on the reference used and 828 (2%) of the sequences were mapped to unknown FIGURE 1 | Frequency distribution of some SDS traits. The combined data were from three locations and over 2 years field evaluation of 180 wheat genotypes. The right and left ends of the graphs indicate the highest and lowest affection classes, respectively. The combined Septoria disease severity at heading, mid-maturity, and maturity stages followed a normal distribution. The severity values increased from heading to mid-maturity, and then to maturity stages. Combined Septoria progress coefficient showed pseudo-normal distribution, confirming the quantitative nature of STB resistance in the tested wheat material. The x-axis represents the BLUE value of the study genotype. chromosomes in the reference. In contrast, approximately 11% (3,792) of the SNPs did not align to any of the chromosomes of the wheat reference genome. Furthermore, the discovered DArTSeq SNPs were not evenly distributed among the subgenomes, with the A, B, and D genomes accounting for 10,317, 10,979, and 9,756 SNPs, respectively (Supplementary Figure 1). Among the 21 wheat chromosomes, the highest (2,065) and the lowest (833) numbers of SNPs were mapped to chromosomes 7D and 4D, respectively (Supplementary Figure 1), and on average, each chromosome harbored about 1,479 SNPs. Maintaining SNPs with higher call rate (>70%) and MAF >0.05 resulted in 7,776 SNP markers, among which 87.3% had a known chromosome position in the wheat reference genome. Among the filtered SNPs, 2,410 were distributed on the A genome, 2,872 were distributed on the B genome, and 1,506 were distributed on the D genome. The remaining 988 SNPs were assigned to a   ); ***, very highly significant (p ≤ 0.0001); "sn" = non-significant at the α, 5% significance level. hypothetical chromosome "0" for the sake of analysis. Hence, 7,776 SNPs were used in downstream analyses, which included principal component analysis (PCA), population clustering, population structure, LD, and GWAS.

Population Structure Analysis
The STRUCTURE analysis indicated two sub-populations in the association panel (Figure 2A), where ∼43% (76) of the genotypes were assigned to cluster one and 57% (101) were assigned to cluster two. Additionally, the Clumpak result detected a greater degree of genetic admixture between the two sub-populations ( Figure 2B), where all the individual genotypes shared alleles inherited from both subgroups (Figure 2C), thus confirming the presence of close relationships among the study materials. Furthermore, the PCA results also suggested the presence of two sub-populations (Figure 3). The first two principal components  explained 65% (PC1 = 50% and PC2 = 15%) of the total variance contained in the data (Figure 3). With this, a scree plot, which was used to display the proportion of variation captured by each of the 10 principal components, also showed that the first two principal components (PC1 and PC2) explained the highest proportion of the total variation in the panel ( Figure 4A). Figure 4B represents the 3D plots of the first three principal components to depict the samples' relationship in space, the analysis also confirmed the presence of kinship in the association panel (Figure 4C), suggesting the importance  of using a powerful statistical GWAS model that accounts for the population structure and familial relatedness in the association study.

Linkage Disequilibrium (LD) Analysis
The linkage disequilibrium of alleles at different loci varied considerably across each chromosome and among chromosomes and sub-genomes (  0.0251) and 2D (r 2 = 0.211), respectively ( Table 7). The physical distance (bp) at which the LD decayed to the critical r 2 (0.2) value was used to determine the confidence interval for declaring the distinct QTL for each chromosome. Significant SNP markers from the same chromosome were also assigned to the same QTL if the distance between the significant markers was less than the critical physical distance.

STB Resistance
The GWAS identified 53 SNPs that were significantly (FDR < 0.05) associated with STB resistance at any growth stage. The report, however, only included the marker-trait associations (MTAs) that surpassed a Bonferroni-correction significance threshold of 0.15. Supplementary Table 3 presents the complete output of the GWAS results for STB resistance at the heading, mid-maturity, and maturity stages and for the Septoria progress coefficient. This table also reports allele identity, marker position, MAF, p-values, FDR-corrected q values, and the additive effects of the identified MTAs. The Q-Q plots demonstrating how well the used GWAS model accounted for population structure and kinship for STB resistance analysis are also presented in Supplementary Figure 2.
Among the 53 identified MTAs, 3 did not have chromosomal positions on the bread-wheat physical map  Table 4). Ten (18.9%) of the MTAs conferred STB resistance at heading, 4 (7.6%) were effective at midmaturity, 4 (7.6%) were effective at the maturity stage, and 35 (66.9%) of the MTAs were associated with a resistance to disease development as plant height increased (Supplementary Table 4). The percentage of phenotypic variance explained by the markers varied considerably, from 2.7% for the SDS measured at maturity at Bekojiin in 2016 to 13.2% for the severity data measured at the mid-maturity stage at the same location in 2015 (Supplementary Table 4). The proportion of phenotypic variation (R 2 ) explained by SDS MTAs at heading ranged from 2.9% for the allele 1195254|F|0-31:C>T-31:C>T on chromosome 3A to 11.1% for 1087857|F|0-41:T>C-41:T>C on chromosome 7D (Supplementary Table 4). Likewise, the R 2 for MTAs for SDS at mid-maturity, maturity, and SPC ranged from 8.6 to 13, 2.7 to 2.7, and 6 to 10.8%, respectively (Supplementary Table 4).

(Supplementary
The combined measure of SDS at the heading, mid-maturity, and maturity stages did not provide any significant associations at the used threshold. However, the combined measure of SPC identified eight MTAs at the stringent Bonferroni significance threshold on chromosomes 1B, 2D, 3A, 3B, 3D, 6B, 7B, and 7D, with one MTA that was unmapped on the bread wheat physical map (Supplementary Table 4,  Supplementary Figure 3).
The GWA scan for SDS at the individual-environment level identified considerable (45) MTAs conferring resistance to STB at different growth stages (Supplementary Table 4, Supplementary Figures 4, 5). The analysis for disease data measured in 2015 at Holetta identified six MTAs for STB resistance at heading on chromosomes 1D, 2A, 3A, 3D, 5A, and 7D, four MTAs effective for STB resistance at the mid-maturity stage on chromosomes 1B, 3D, and 7B, with one MTA with an unknown position on the bread wheat physical map at Holetta, and nine MTAs for SPC (six at Holetta and three at Kulumsa) on chromosomes 1A, 1B, 1D, 2B, 3A, 3D, and 7B (Supplementary Table 4, Supplementary Figure 4). However, no MTA was observed for SDS data measured at the maturity stage in the same year. Likewise, the association analysis for SDS data measured in 2016 identified 26 MTAs: 4 for STB resistance at heading at Bekoji on chromosomes 3A, 3D, and 7A, with 1 MTA with an unknown position on the bread wheat physical map, 4 for maturity stage resistance on chromosomes 1D, 4A, 6A, and 7D at the same location, 18 MTAs for SPC, which were mapped to chromosomes 1A, 1B, 2B, 2D, 3B, 5B, 5D, 6A, 6B, 7A, and 7D plus 1 MTA with an unknown position on the bread wheat physical map (Supplementary Table 4  Putative QTL for STB resistance were identified by combining the MTAs based on their genomic positions using a window of physical distance (in Mbp) determined through a pair-wise LD analysis of the genome-wide scanned SNPs. Supplementary Figure 6 presents a scatter plot of the genome-wise pairwise LD r 2 values between the SNPs on each chromosome against inter-marker physical distance. The MTAs falling on the same linkage group within the physical distance for LD decay specific for that chromosome were assigned to the same putative QTL. Hence, based on the LD criteria, the 53 markers were assigned to 33 putative QTLs (Table 8, Figure 5).
The functional association of the identified QTLs for STB resistance was further investigated by annotating genes found in the QTL regions on the recently released IWGSC RefSeq Annotation v2. The analysis resulted in several disease resistance-associated genes involved in plant defense systems (Supplementary Table 7). For instance, the highconfidence candidate genes TraesCS1A02G279300 on 1A, TraesCS1B02G332400 on 1B, TraesCS1D02G001700 on 1D, and TraesCS2A02G297500 on 2A are highly involved in systemic acquired resistance (SAR), which refers to the long-lasting, broad-spectrum resistance of plants to pathogen infections. Specifically, the high-confidence gene detected near qSTB.08 on chromosome 2A (TraesCS2A02G297500) controls mitogenactivated protein kinase (MAPK) cascades, which are involved in signaling multiple defense responses of plants against pathogen attacks (Meng and Zhang, 2013).

MTAs for Agronomic Traits
The combined measures of the agronomic traits, such as days to heading, days to flowering, days to maturity, grainfilling duration, grain yield, and 1,000-kernel weight, did not provide any MTAs at the stringent Bonferroni significance threshold used except for plant height, which resulted in one MTA on chromosome 7A at the 514.43 Mbp position (Supplementary Table 6). However, the dissections of these traits in the individual environments in separate years identified considerable MTAs at the significance threshold utilized (Supplementary Figure 7). In particular, a GWA scan for days to heading in 2015 at Holetta provided six MTAs on chromosomes 1A, 5A, 5B, 6A, 6B, and 7B (Supplementary Table 6). The same trait measured in 2016 at Kulumsa identified three significant (FDR < 0.05) SNPs on chromosomes 2B, 5D, and 6A (Supplementary Table 6). Additionally, the total phenotypic variance explained by the SNPs for this trait ranged from  Table 6). The GWA scan for days to flowering in the individual environments identified six MTAs for the data measured in 2015 at Holetta on chromosomes 1A, 5A, 5B, 6A, 6B, and 7B (Supplementary Table 6, Supplementary Figure 8). The detected markers could explain 22.8-25% of the total phenotypic variation for days to flowering, while the SNP markers 1000134|F|0-15:T>C-15:T>C on chromosome 6B and 1204551|F|0-57:C>T-57:C>T on chromosome 1A accounted for the lowest and highest phenotypic variations in days to flowering in the association panel, respectively (Supplementary Table 6). Likewise, the association analysis for the days to maturity data measured in 2016 at Bekoji identified 15 MTAs pointing to nine putative QTLs on 1A, 2A, 3A, 3B, 5B, 6D, and 7A (Supplementary Table 6, Supplementary Figure 9). Three of these significant SNPs, however, were not mapped on the bread wheat physical map. The detected markers explained 11.6-15.6% of the total phenotypic variation for days to maturity, while the lowest and highest values were reported for the SNPs 2278215|F|0-18:A>G-18:A>G on 2A and for 7332831|F|0-9:T>C-9:T>C on 1A, respectively (Supplementary Table 6).
Moreover, a GWA scan on grain-filling duration data measured in 2016 revealed 15 MTAs, among which 11 were identified from the data collected at Holetta on chromosomes 1B, 2B, 3A, 3B, 6D, and 7D, plus 1 MTA with an unknown position on the bread wheat physical map (Supplementary Figure 9). Four of the significant associations for this trait were obtained from data measured at Kulumsa on chromosomes 2B, 3B, 5B, and one MTA with an unknown chromosomal position on the bread wheat genome. The total phenotypic variation for GFD explained by the SNPs ranged from.5% for the allele on chromosome 3B (5325155|F|0-26:G>T-26:G>T) to 7.9% for the SNP marker on 3A at 250.82 Mbp (5325155|F|0-26:G>T-26:G>T) (Supplementary Table 6).
The association analysis for pooled plant height data identified 1 MTA on chromosome 7A (Supplementary Figure 9) and 24 for the data measured in the individual environments (Supplementary Table 6). The plant height data measured in 2015 at Bekoji resulted in 4 MTAs on chromosomes 1A, 5A,7A, and 7B and 10 MTAs for Kulumsa on chromosomes 1B, 2A, 5B, 5D, 6B, and 7D (Supplementary Figures 7, 8).
The study revealed that grain yield and yield-related attributes measured over all the environments did not provide any significant SNPs. However, their association analysis based on mean values in the individual environments identified numerous MTAs. The grain yield measured in 2015 at Kulumsa identified 10 significant MTAs on chromosomes 1A, 3A, 3D, 4A, 5A, 5D, and 7D, with one SNP that was not mapped in the bread wheat genome (Supplementary Table 6,  Supplementary Figure 8). Similarly, a GWA scan for TKW measured in 2016 at Kulumsa provided one MTA pointing to a putative QTL on 7D at 79.52 Mbp (Supplementary Table 6,  Supplementary Figure 9). Furthermore, the total phenotypic variations explained by the significant SNPs for yield and yieldrelated traits varied considerably based on traits and markers, with the lowest being 0.5% for grain yield on 3D (1045011|F|0-60:T>A-60:T>A) and the highest being 15.7% for TKW on 7D (4262368|F|0-28:A>G-28:A>G) (Supplementary Table 6).
The analysis revealed that some putative QTLs identified by SDS data overlapped with those determined for agronomic data. For instance, the putative QTL determined for plant height measured in 2016 at Kulumsa on chromosome 6B was comapped with qSTB.28, which was identified for combined SPC. Similarly, the putative QTL identified on 7B for plant height measured at Bekoji in 2015 was co-mapped with qSTB.31, which was also identified for combined SPC. In addition, the putative QTLs mapped on chromosome 7D for plant height measured in 2015 and 2016 at Kulumsa and for GFD and TKW data measured in 2016 at Holetta and Kulumsa, respectively, were co-mapped with qSTB.32, which was identified for the SDS data measured at the maturity stage at Bekoji and for pooled SPC data.

DISCUSSION
The phenotypic evaluation revealed significant genetic variability for STB resistance among the tested wheat genotypes, thus confirming the availability of relevant alleles for future breeding and improvement. The observed high broad-sense heritabilities within the individual environments (H 2 = 0.58-0.99) and across the environments (H 2 = 0.72-81) indicated a strong genetic signal in the data, which can be used for improving STB resistance through selection. Similarly, a high broad-sense heritability (H 2 = 0.78) for STB resistance was reported for European bread wheat varieties in Germany (Muqaddasi et al., 2019) and Tunisia (H 2 = 0.55) (Berraies et al., 2014). The present field evaluation results confirmed that STB infestation is significantly influenced by year, location, and all interaction effects, thus confirming the importance of multiple locations and years for germplasm evaluations at disease hotspots to identify durable and stable STB-resistant genotypes.
The correlation analysis revealed that Septoria disease ratings had negligible negative correlations with plant height, indicating that tallness only had a weak contribution for reducing STB infections (Muqaddasi et al., 2019). The lack of or slight negative correlations of SDS traits with the agronomic traits HD, FD, GFD, NKPS, and NKS and the moderately negative correlation with MD indicate that genotypes with late phenology could escape STB with slightly reduced infection. Moreover, the significant negative correlations of STB infection with yield and yieldrelated traits such as HLW, TKW, and KN could most likely be due to the fact that the infection of the flag and second leaves at the grain-filing stage could significantly influence the photosynthesis process, and, thus, result in reduced grain yield. Negative associations of SDS with days to flowering, days to maturity, number of seeds per spike, and thousand-grain weights were also reported for Ethiopian durum wheat (Kidane et al., 2017).
The STRUCTURE and principal component analyses revealed population stratifications and admixtures, suggesting the need to use a powerful statistical model in the association analysis that controls for spurious marker-trait associations. The analyses suggested two sub-groups in the population. The very powerful statistical model used in the analysis, FarmCPU, sufficiently accounted for population stratification, familial relatedness, and marker effects, which consequently reduced the confounding effects that could result in false-positive MTAs. This was confirmed by visualizing the Q-Q and Manhattan plots. Similar indistinct population stratifications, higher admixtures, and weak population sub-structuring were reported among 371 European wheat genotypes based on 35k and 90k SNP marker arrays (Muqaddasi et al., 2019).
Like previous reports, this study confirmed the unequal distribution of the SNP markers among wheat genomes, where most were harbored by the A (10,317) and B genomes (10,979), while fewer SNPs (9,756) were harbored by the D genome (Berkman et al., 2013;Edae et al., 2015;Rahimi et al., 2019). This variation most likely resulted from the evolutionary and domestication history of the crop (Dvorak et al., 2006;Jordan et al., 2015), where the D genome had less time to accumulate mutations.
These analyses revealed that the LD between the markers and genes contributing to STB resistance declined to r 2 < 0 within a physical distance of 1. 26-105.61 Mbp in all the chromosomes, with an overall mean of 31.44 Mbp. This is much lower than the average physical distance (69.1 Mbp) for LD decay in Ethiopian durum wheat at the critical threshold of r 2 = 0.2 (Alemu et al., 2021). The marker distances at which the LD decayed across the older sub-genomes (A and B) were relatively lower than those for the D sub-genome, most likely because of the long evolutionary history of the A and B genomes as compared with the D genome. Furthermore, the LD between alleles can decay because of a number of factors such as selection, recombination, the mating system, genetic drift, mutation, and/or population relatedness (Stich and Melchinger, 2010). Hence, it is likely that the shorter selection history of the D sub-genome did not allow linkage breakdown due to the recombination that occurs between SNPs located at longer physical distances.
The GWAS analysis identified 53 MTAs pointing to 33 QTL for STB resistance and 82 MTAs for agronomic traits where markers had a FDR p ≤ 0.05 and a Bonferroni correction significance threshold of 0.15. The number of SDS MTAs identified in this study was significantly lower than that in the findings of Rahimi et al. (2019), who reported 313-394 MTAs for an Iranian bread wheat association panel. However, this number was still substantially higher than that in the report of Kidane et al. (2017), who identified 35 significant associations for an Ethiopian durum wheat panel. Only seven QTLs for SPC were identified in an analysis of the mean over environments, while none was detected for SDS. Kidane et al. (2017) also reported no QTLs for SDS in an analysis of means.
More QTL were noted in the analysis of data from individual environments, although none was detected in more than four of the six environments and 24 of the 33 were detected in just one environment. The failure to detect QTL effects over environments could be due to the seasonal specificity of QTL effects, disease pressure, or the use of a very stringent FDR level that controls type I errors but leads to increased type II errors, e.g., declaring a QTL not significant when it actually is. The 2015 growing season at Holetta was characterized by extended and heavy rainfall that resulted in the highest STB natural infestations across all the growth stages. Additionally, in this growing season, 12 individual QTLs were detected (five for SPC and seven for SDS), 3 of which affected both SPC and SDS. In contrast, the heavy rainfall and prolonged moisture experienced at Bekoji in 2016 produced the highest SDS, while 12 QTLs were detected, of which 6 were for SPC, 4 for SDS, and 2 that affected both traits, using the data that were obtained. Therefore, climatic conditions, such as persistent crop moisture and prolonged heavy rain, favor the successful infection and spread of the disease throughout the crop canopy (Fones and Gurr, 2015). Furthermore, no QTLs for SDS were detected for Kulumsa in either year. Although a total of 11 QTLs for SPC were detected in the same environment, none was repeated over the years. The different climatic conditions may have caused the later onset of the disease in both growing seasons. However, one QTL, qSTB.21, had the most repeatable effect and was significant for SPC overall for both SPC in two environments SDS at the heading and mid-maturity stages in 2015 at Holetta.
The identification of defense-related candidate genes, such as TraesCS1A02G279300, TraesCS1B02G332400, TraesCS1D02G278400, TraesCS2B02G233600, TraesCS2A02 G297500, TraesCS2D02G497400, TraesCS2D02G506300, and TraesCS4A02G341300, in the vicinity of the significant markers indicates the possible functional association of the detected QTL regions in plant defense systems against pathogen infections. For instance, the translations of the genes found in qSTB.02 (TraesCS1A02G279300) and qSTB.04 (TraesCS1B02G332400) on chromosomes 1A and 1B, respectively, are involved in the jasmonic-acid and ethylene-dependent systemic-acquired resistance of plants to pathogen infections. This systemic acquired resistance (SAR) is a broad-spectrum, long-lasting resistance acquired after the initial localized infection of plants by pathogens (Lawton et al., 1995). Furthermore, the gene TraesCS2A02G297500 found in qSTB.08 on chromosome 2A controls MAPK cascades, which are involved in signaling multiple defense responses, including the biosynthesis/signaling of plant stress/defense hormones, reactive oxygen species (ROS) generation, stomatal closure, defense gene activation, phytoalexin biosynthesis, cell wall strengthening, and hypersensitive response (HR) cell death. Moreover, most of the genes in proximity to the detected significant markers are inferred to be involved in salicylic acid (SA) biosynthesis, with SA being an important plant hormone that is best known for mediating host responses upon pathogen infection (Lefevere et al., 2020).
In this study, we discovered some STB resistance QTL that appear to be novel. These include the putative QTLs on chromosomes 1A (qSTB.01-3), 5D (qSTB.26), and 6B (qSTB.28) that explained >5% of the genetic variations, suggesting their relevance for wheat resistance breeding against STB. To the best knowledge of the authors, none of the known major STB resistance genes published in existing literature have been mapped to these regions of the wheat chromosomes; therefore, these QTLs could be considered novel.
Moreover, the study revealed that some of the putative STB resistance QTLs were co-located with QTL for agronomic traits. For instance, the putative QTLs derived from plant height measured in 2016 at Kulumsa on chromosome 6B (R 2 = 11.36) and in 2015 at Bekoji on chromosome 7B (R 2 = 8.88) were co-mapped with qSTB.28 and qSTB.31, which were identified for combined SPC. Likewise, the putative QTLs mapped on chromosome 7D for grain-filling duration (R 2 = 4.69) and 1,000kernel weight measured in 2016 (R 2 = 0.54) at Holetta and Kulumsa, respectively, were co-mapped with qSTB.32, which was identified for the SDS data measured at the maturity stage at Bekoji and for the pooled SPC data. Furthermore, STB is the most destructive foliar disease in Ethiopia. Hence, the infection of the flag and second leaves, which contributes most to photosynthetic assimilates at the grain-filling stage (King et al., 1983;Muqaddasi et al., 2019), can result in the substantial loss of grain weight and yield. This is consistent with the findings of Kidane et al. (2017), who reported the co-mapping of putative QTLs for 1,000kernel weight with SDS data. Moreover, the putative QTL on chromosome 1A identified for grain yield measured at Kulumsa in 2015 (R 2 = 0.79) was co-mapped with qSTB.03, which was identified for the SPC measured in the same environment. It is, therefore, expected that the vertical progression rate of the disease could affect grain yield by influencing grain filling and the number of seeds produced per spike.
In this study, most of the QTLs identified for agronomic and phenological traits did not overlap with those detected for SDS traits, likely because of the lack of common genetic effects for STB resistance and these traits. Many of the correlations of STB traits with agronomic traits were non-significant and had negligible to weak negative coefficients, indicating that the traits were independent. However, some level of co-localization was observed for the putative QTL for days to heading and days to flowering measured at Holetta in 2015 on chromosome 6B (R 2 = 22.82), with the putative QTL qSTB.28 on 6B being identified for the SPC measured at Bekoji in 2016 and for the pooled data.

CONCLUSIONS
In this study, the genetic architecture of adult-plant resistance to STB was explored in bread wheat using high-density, genome-wide SNP markers and multi environment-derived phenotype data. The analysis revealed that the association panel possessed considerable STB resistance alleles that could be deployed to improve wheat resistance to the prevailing Z. tritici populations in Ethiopia. Several genotypes with better resistance than the moderately resistant check King-bird were identified. Furthermore, the GWAS identified 33 putative QTLs, which were associated with 53 SNPs. Most (24) of the QTLs were detected in just one environment, suggesting the presence of resistance gene/genes effective against location-specific Z. tritici races. The detected QTLs also explained 2.7-13.2% of the total phenotypic variance for STB resistance. Several disease resistance-associated gene/s were also identified in the proximity of the detected SNPs, which can be targeted in efforts to understand the actual causative genes at the associated loci. Additionally, most of the detected putative QTLs shared similar chromosomal positions with previously reported genes and QTLs. Among these detected alleles, five were potentially novel, accounting for >5% of STB resistance. However, the effects of these QTLs need to be validated before being deployed in MAS. Finally, we conclude that the identified stably resistant wheat genotypes and the identified QTLs can be deployed in wheat breeding programs for the development of durable and broad-spectrum-resistant varieties against Z. tritici.

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are publicly available. This data can be found here: Dryad data repository doi: 10.5061/dryad.7sqv9s4s4.

AUTHOR'S NOTE
We confirm that the manuscript has been read and approved by all the named authors and that there are no other persons who satisfied the criteria for authorship but are not listed. We further confirm that the order of authors listed in the manuscript has been approved by all of us.

AUTHOR CONTRIBUTIONS
TM was involved in the conceptualization, methodology, statistical analysis, and writing of the original draft. CS was involved in genotyping and editing. TH was involved in data curation and visualization. BA contributed material, performed the field experiments, and edited. CZ was involved in investigation and validation. DL was involved in reviewing. SG was involved in editing. KT was involved in methodology, funding acquisition, project administration, and supervision. All the authors approved the final version of the manuscript.

FUNDING
This study was supported financially by the Ministry of Innovation and Technology (formerly Ministry of Science and Technology) of the Federal Democratic Republic of Ethiopia.