Genome-wide association mapping in wheat reveals novel QTLs and potential candidate genes involved in resistance to septoria tritici blotch


 Septoria tritici blotch (STB) caused by Zymoseptoria tritici is one of the most important foliar diseases of wheat that can lead to significant yield losses worldwide. In this study, we analyzed STB resistance pattern of 185 diverse bread wheat genotypes at the seedling stage in order to confirm known resistance genes or identifying novel putative QTLs conferring resistance to STB. The phenotyping data using 10 Z. tritici isolates showed that large genetic variance exists for STB infection in wheat germplasm. Genome-wide association mapping (GWAS) using CMLM algorithm identified 37 quantitative-trait loci (QTL) for resistance to STB on 17 chromosomes. Most of QTLs overlapped with known STB resistance genes, while GWAS scan revealed several newly identified QTLs on chromosomes 4A, 5D and 6D for resistance to STB. Genome annotation of marker-trait association (MTA) against wheat reference genome revealed 29 MTAs corresponded to the putative genes, which their functional descriptions identifiers included diseases resistance proteins in plants. These MTAs and novel QTLs associated with STB resistance found here, can be used in future wheat breeding programs to recombine different loci for durable resistance to STB.


Introduction
Common wheat (Triticum aestivum L.) is the most important cereal crop in the world and plays important role in the diets of humans and livestock. Average global wheat production is about 766.4 million tons in 2019 (http://www.fao.org/worldfoodsituation/csdb/en). This makes wheat the third important crop in terms of production after maize and rice. Global wheat production can be negatively in uenced by abiotic and biotic stresses. Septoria tritici blotch (STB) caused by Zymoseptoria tritici is one of the most destructive fungal diseases of wheat worldwide (Kema et al. 1996;Hardwick et al. 2001). The fungus causes expanded necrotic lesions from early-emerging leaves to ag leaves, and the overall damage caused by STB can result in grain yield losses up to 50% under optimal environmental conditions (Mehrabi et al. 2006;Goodwin 2007;Kema and van Silfhout, 1997;Suffert et al. 2011). Under epidemic conditions, fungicide application is required to control STB, but this strategy is not adequately effective due to quick adaptation of the pathogen to fungicides by mutation (Torriani et al. 2009;Mohammadi et al. 2017;Kema et al. 2018). In addition, the application of fungicides has severe threats to human health and the environment. Therefore, the characterization of new resistance sources and the development of resistant wheat cultivars is the most economical and environmentally-friendly approach and fundamental strategy in breeding programs for sustainable agriculture and food security (Talebi et al. 2010;Mehrabi et al. 2015). Up to now, 22 major genes and 89 QTLs for resistance to Z. tritici have been reported and mapped on the wheat genome along with their closely linked markers (reviewed by Brown et al. 2015;Yang et al. 2018).
Both qualitative and quantitative resistance has been reported for resistance to STB in wheat. Qualitative resistance is usually controlled by major genes that confer complete resistance and follows by the gene-for-gene model (Brown et al. 2015). This model has been reported for the rst time on Z.tritici isolate IPO323 and cv. Flame (Kema et al. 2000;Brading et al. 2002). Quantitative resistance is controlled by many minor genes and widely reported in wheat cultivars at both seedling and adult growth stages (Arraiano and Brown, 2006;Chartrain et al. 2004;Goodwin, 2007).
Although the qualitative resistance is effective and important for controlling STB, the rapid evolution and adaptation of Z.tritici populations result in overcoming resistance in most improved cultivars (Cowger et al. 2000;Stukenbruke et al. 2007;Makhdoomi et al. 2015;Muqaddasi et al. 2019).
Therefore, in order to increase the durability of wheat resistance to STB, combining qualitative and quantitative resistance genes is required to preserve resistance effectiveness over time (Brown et al. 2015;Vagndorf et al. 2017). Availability of high-de nition genotyping using genotype-bysequence (GBS) technologies enables breeders to identify resistance genes in diverse germplasm that can be employed for the development of new resistant cultivars using the gene pyramiding approach (Kidane et al. 2017;Vagndorf et al. 2017;Arraiano and Brown, 2017;Muqaddasi et al. 2019).
The genetic architecture of resistance to STB has been mostly evaluated in different bi-parental populations for detection of large-effect genes or quantitative trait loci (QTLs) (Chartrain et a. 2004(Chartrain et a. , 2005aSimon et al. 2005;Tabib Ghaffary et al. 2011. Genome-wide association mapping (GWAS) using a large number of markers with high genome coverage is a powerful tool for detecting the resistance loci associated with diverse germplasm possessing natural variation of resistance genes (Bartoli and Roux, 2017). Various GWAS studies have been performed on diverse wheat germplasm (Mirdita et al. 2015;Vagndorf et al. 2017;Kidane et al. 2017;Gerard et al. 2017;Muqaddasi et al. 2019;Arraiano and Brown, 2017). Most previous GWAS studies for revealing QTLs against STB have been done for detecting resistance loci against a mixture of isolates under natural eld conditions. Therefore, GWAS analysis for speci c-isolate resistance using globally diverse pathotypes of Z. tritici is of great interest for breeders to identify new genes/loci that can be used in wheat breeding programs.
The present study relies on genotyping and phenotyping of a globally diverse panel of 185 wheat genotypes that were evaluated at the seedling stage against 10 Z. tritici isolates (collected from Iran, Algeria, Turkey, France and The Netherlands). DNA ngerpring of the wheat genotypes was performed using high-throughput DArTseq technology (SilicoDArT and SNP markers) aiming to identify and localize possible novel QTLs associated with isolate-speci c resistance against Z. tritici in wheat.

Plant material and evaluation of STB infection
The 185 wheat genotypes comprised 64 Iranian improved cultivars, 103 accessions from globally diverse origins, and 18 wheat differential lines with known STB genes (Supplementary Table S1). Phenotypic evaluation of wheat genotypes against 10 Z. tritici isolates (Table 1) was described in detail by Mahbubi et al. (2020). Phenotyping of wheat genotypes was done under greenhouse conditions at the seedling stage as described previously Makhdoomi et al. 2015). Brie y, STB isolates were pre-cultured in yeast-extract glucose (YG) liquid medium, and then these pre-cultures were used to inoculate 100 ml YG media per isolate. YG media were placed in an orbital shaker (set at 125 rpm) and incubated at 18 °C for 5-6 days. The propagated yeast-like spores were collected and their concentrations were adjusted to 10 7 spores/ml. Five seeds of each genotype were sown in plastic pods and 10-day old seedlings were spray-inoculated using hand sprayers. Inoculated plants were kept in dark plastic bags for 48 h at 20-25 ºC and then transferred to a greenhouse with the environmental condition of 18-22 ºC and >85% humidity (Kema et al. 1996). After 3-weeks, the plants were scored visually based on the percentage of necrotic leaf area covered by pycnidia as described before (Kema et al. 1996;Ghaneie et al. 2012). Experiments were set up with a randomized complete block design with three replications.
Phenotypic data analysis Data were normalized using arcsin square root-transformation (Sokal and Rolhf, 1995). Normalized data were analyzed using a linear mixed model (LMM), in which the isolates and wheat genotypes were the main effects and their two-way interaction as xed effects. The effects of replication (block) and experiment were considered as random effects. Broad sense heritability estimates were calculated for best linear unbiased estimates (BLUE) values from all STB isolates. BLUE values were then used for cluster analysis, correlation analysis and also to perform GWAS.
Cluster analysis of wheat genotypes according to pycnidia coverage data by each isolate was performed using the un-weighted pair-group method (UPGMA) and the dissimilarity matrix was measured using the Ward's method implemented in PAST software (Hammer et al. 2001). Pearson correlation coe cients among 10 STB isolates were calculated for the percentage of leaf covered by pycnidia in 185 wheat genotypes.
Genotypic data DNA was extracted from seedling plants of each genotype using the CTAB protocol (Lasner et al. 1989). The quality and quantity of DNA were checked by spectrophotometer, and DNA concentration was adjusted to 100 ng/µl. DNA samples were plated in 96-plex and shipped to DArT Pty Ltd, Canberra, Australia, and genotyped using the DArTseq technology (Sansaloni et al. 2011;Alam et al. 2018). The detailed methodology for the generation of DArTseq markers (SNP and SilicoDArT) was described before (Egea et al. 2017). We received a total of 94535 (54309 SilicoDArT and 40225 SNP) markers, which were polymorphic across 185 wheat genotypes. Marker loci with unknown chromosome positions (based on genome assembly) were removed from the analysis, and the remaining markers were ltered using a minimum minor allele frequency (MAF) of 0.01 in TASSEL v.5.2.37 software (Bradbury et al., 2007). Wheat is a self-pollinated crop and we assumed that all genotypes are homozygous. Therefore markers showing heterozygous were indicated as missing and markers with >10% missing were excluded. In total, 21773 (15856 SilicoDArT and 5917 SNP) distributed across the 21 chromosomes were maintained for analysis (Supplementary Table S2).
Linkage disequilibrium and population structure analysis Linkage disequilibrium (LD) for DArTseq markers was implemented in TASSEL v.5.2.37. The critical r 2 -value was determined by root transforming the unlinked r 2 -values and taking the 95th percentile of the distribution as the threshold beyond which LD is likely caused by genetic linkage (Nielsen et al. 2014;Monostori et al. 2017). The graphical LD decay was imputed by the GAPIT R package (VanRaden, 2008;Lipka et al., 2012).
Population structure was performed in STRUCTURE 2.1 based on an admixture model (Evanno et al., 2005). Model run as the burn of 10000 iterations followed by 10000 Markov Chain Monte Carlo iterations was set for accurate parameter estimates. The optimal value of K ranging from 1 to 10, with three independent runs. Principal coordinate analysis (PCoA) and cluster analysis among the wheat genotypes estimated in DARwin ver. 5.0 software using the Unweighted Neighbor-Joining (UNJ) algorithm.

Genome-wide association analysis
Genome-wide association mapping (GWAS) analysis for resistance to 10 Z.tritici isolates in a panel of 185 wheat genotypes was conducted in R package Genome Association and Prediction Integrated Tool (GAPIT) (Lipka et al. 2012) using all 21773 DArTseq markers. Association analyses were done using different models like GLM, MLM, CMLM and Super-MLM models. It has been reported that the MLM model is better than GLM by controlling false positive and also MLM able to control false positive disappear in GWAS analysis for complex traits in a panel of genotypes that show extensive genetic divergence (Yu et al. 2006;Kaler et al. 2020), but in some cases, MLM model may lead to false negatives and consequently weakness of marker-trait associations (MTAs) (Zhang et al. 2010). Compressed mixed linear model (CMLM) improves the statistical power of analysis. Finally, all models were compared for false positives and false negatives based on the Q-Q plots. CMLM showed the best results and nal GWAS analysis was performed in this model.
In order to detect the marker-trait association (MTA) in GWAS, the signi cance threshold for MTA was set to P = 0.05 after applying the false discovery rate (FDR) correction. Phenotypic data for all isolates gave models that did not cross the false discovery rate threshold, and therefore, the arbitrary threshold of signi cance of -log10(P) ≥ 3.0 corresponding to a P-value of 0.001 was accepted as signi cant MTA (Vagndorf et al. 2017;Muqaddasi et al. 2019). Map positions of signi cant MTAs were determined according to their genetic positions in a high-resolution DArT-seq consensus map (version 4.0), provided by Dr. Andrzej Kilian (Diversity Arrays Technology Pty Ltd, Canberra, Australia) and bread wheat IWGSC RefSeq v1.04 with BLAST+ v2.7.1 (Camacho et al. 2009). Overlapping markers on the same chromosome for resistance to different STB isolates were considered to tag a single QTL if their positions were closer than 10 cM and showed LD r 2 >0.3. Then for comparison of the QTLs identi ed in this study, the position of each QTL was projected onto the wheat SSR consensus map (Somer et al. 2004;Maccaferri et al. 2015) for crossreference with previous SSR maps. Each QTL considered new if its position was ≥ 10 cM from previously reported STB resistance genes or QTLs.
Candidate genes for individual MTAs were identi ed by aligning the physical position of markers to the sequence of the wheat genome assembly IWGSC v.1.0 (https://plants.ensembl.org/Triticum_aestivum/Info/Index). High-con dence annotated genes were retrieved from a 10 kb window of left and right of each identi ed MTA.

Results
Phenotypic data analysis Analysis of disease severity data showed signi cant differences (P 0.001) among wheat genotypes for resistance to Z. tritici (Table 2). Genotype × isolate interaction was highly signi cant (P 0.001) and indicated the differences in wheat genotypes for their responses to Z. tritici isolates. Heritability values for disease severity were high for all isolates. The high heritability values indicated that there was a limited replication variation relative to genotypic variation for all isolates. This is supported by high Pearson correlation coe cients for disease severity between all isolates (Table 3). The Pearson correlation coe cient between STB isolates ranged from 0.26 (IPO02166 and IPO92034) to 0.90 (IPO323 and IPO86013) with an average value of 0.61 (Table 3). In total, 239 isolate-speci c resistances were found among all interactions (n=1850), of which 183 showed disease severity P≤5% (highly resistance), 56 were disease severity 5<P≤10% (resistance) (Supplementary Table S1).
Cluster analysis and principal component analysis using omitted data from the percentage of leaf area covered by pycnidia, grouped wheat genotypes in three distinct clusters (Fig 1). The rst cluster comprised 63 wheat genotypes, of which 19 genotypes were Iranian cultivar and breeding lines, 7 genotypes with previously known Stb genes (Bulgaria, Israel 493, Cs-Synthetic, Sha r, Es-Federal, M6 and Balance) and the remaining genotypes originated from diverse origins. Resistance spectra of genotypes for this cluster ranged from 30.98% (IPO323) to 56.61% (IPO02166) ( Table 4). Cluster-II contained 85 genotypes, with 32 Iranian cultivars, one differential cultivar (Curtot) and 52 genotypes from diverse origins. All genotypes in this cluster showed a highly susceptible pattern to most of the isolates and the mean disease severity ranged from 45.14% (IPO323) to 66.57% (IPO02166) ( Table 4). Cluster III comprised 37 genotypes, which most of the differentials with known Stb genes (Veranopolis, Tadinia, Kavkaz-K4500, TE9111, Salamoni, Arina, Riband and M3) grouped in this cluster. In general, genotypes belonging to this cluster showed low mean disease severity ranged from 7.76% (IPO323) to 43.59% (IPO02166) ( Table 4). Cluster-III comprised of 37 genotypes, including 13 Iranian cultivar and breeding lines, eight differentials (Veranopolis, Tadinia, Kavkaz-K4500, TE9111, Salamoni, Arina, Riband and M3) and 16 landraces from different origins. In general, genotypes grouped in this cluster had a low mean disease severity ranging from 7.76% (IPO323) to 43.59% (IPO02166) ( Table 4).

Population structure and linkage disequilibrium
Extensive genotyping on 185 wheat genotypes resulted in 21773 (15856 SilicoDArT and 5917 SNP) markers. The Unweighted Neighbor-Joining cluster analysis (Fig 2a) and Bayesian model-based structure analysis (Fig 2b) grouped wheat genotypes into four distinct subpopulations. Subpopulation 1 (50 genotypes) consisted of 17 Iranian landrace genotypes and other genotypes that were mainly from Western Asia and Eastern Europe (Turkey, Tajikistan, Romania, Hungary). Sub-population 2 (29 genotypes) consisted mostly of genotypes from diverse origins and two Iranian cultivars. Sub-population 3 (86 genotypes) comprised most of the Iranian improved genotypes (45 genotypes) and eight wheat differentials for resistance to STB. Sub-population 4 (20 genotypes) consisted of landraces accessions from different origins. In general, there was no strong signi cant relationship between cluster grouping and origin of wheat genotypes, most probably due to the ancient and recent international exchange of germplasm, while the Iranian improved cultivars mainly grouped in the same cluster. Several of the genotypes used in this study have been utilized as parental lines or have the same pedigree background, so a mixture of origins was observed in all clusters. Nevertheless, a clear distinction on the above-mentioned four main subpopulations was clearly observed based on molecular data analysis. Principal component analysis (PCA) was used to con rm the results of population structure and also showed a distinct pattern of subpopulations. The rst two PCs represented 66.72% of the total variation (Fig 2c). A comparable result similar to cluster analysis and PCA was also observed by the heatmap plot of kinship matrix where four distinct clusters were identi ed (Fig 3). Different subpopulations showed different resistance levels for most of the isolates. Subpopulation 1 (mostly of Iranian landraces and landraces from West Asia) has the highest susceptibility (mean DS=48.2), followed by the subpopulation 2 (mean DS=43.4) and subpopulation 3 and 4 ( mean DS are 43.3 and 40.4, respectively) (Fig 4). These associations between population structure and STB resistance indicated that corrections for population structure were required for association mapping analysis.
In LD analysis the square of the correlation coe cient of alleles between loci (r 2 ) was not signi cant for most of the pairwise comparisons, whereas out of 1048575 marker pairs only 297546 (28%) intra-chromosomal pairs showed a signi cant level (P<0.001) of LD. Marker pairs on the genome-B showed a higher number of signi cant pairs in comparison to the genome-A and genome-D. Mean and critical r 2 values were 0.09 and 0.16, respectively. LD declined with increasing physical genetic distances between markers and r 2 value falling below the critical value over distances of 1.6 kb (Fig 5).

Genome-Wide Association Analysis
Among different GWAS models, CMLM model shows reliable results and presented low spurious associations. For the analysis, phenotypic data (percentage of necrotic leaf area covered by pycnidia) and 21773 mapped DArTseq markers on 185 wheat genotypes were used. The highest number of markers mapped on A genome (8031) and B genome (9537) compared to those mapped on D genome (4205). Association analysis was performed separately for each isolate. A total of 67 signi cant MTAs (P<0.001) were identi ed for STB resistance on 17 chromosomes (Table   5). Manhattan plots for the association between markers and STB responses in different isolates were displayed in Supplementary Figure 1. The QTLs identi ed by different isolates but located at overlapping genomic region on a chromosome were considered a single QTL and assigned the same name using the nomenclature Qstb-iaufollowed by the number of QTL in chromosome order (Table 5; Supplementary Figure S2), and nally 37 QTLs were detected on 17 chromosomes (Table 5). Overall, most QTLs showed small effects and the proportion of phenotypic variation explained by each MTA for disease severity (% of leaf area covered by pycnidia) ranged from 7.28 to 13.74% (Table 5).
The position of these QTLs was compared to the position of mapped QTLs and known genes reported in previous studies (Summarized in Supplementary Figure 2). Results showed that 28 QTLs identi ed in overlap regions or very close to genomic regions in previously reported QTLs/genes. Out of 37 QTLs, ten QTLs were identi ed as novel QTLs for resistance to STB. Two QTLs on chromosomes 1A and 1B (Qstb-iau-2 and Qstb-iau-3, respectively) were mapped far from any reported known STB genes or QTLs on these chromosomes. Qstb-iau-12 on chromosome 2D and another three QTLs (Qstb-iau-17, Qstb-iau-18 and Qstb-iau-19) on chromosome 4A are localized far from previously mapped QTLs and genes. Four other QTLs (Qstb-iau-22, Qstb-iau-23, Qstb-iau-26, Qstb-iau-32) also were mapped on chromosomes 5B, 5D and 6D. These QTLs mapped in novel genomic regions as compared to previously known genes/QTLs for resistance to STB.

Putative candidate gene identi cation
All the MTAs associated with resistance to STB isolates identi ed by CMLM method were mapped to the wheat physical genome (described in the materials & methods section). The physical reference genome of Chinese Spring cv. was used to survey the genes in the anking regions on each MTA (IWGSC RefSeq v1.04 with BLAST+ v2.7.1). We restricted the gene survey to a 10 Kb region, covering 5 Kb regions toward the left and right side of the MTA markers to survey the genes more precisely. Out of 68 MTAs, 39 MTAs showed the sequences on wheat reference genome without characterized molecular function. However, 29 MTAs corresponded to the putative genes, which their functional descriptions identi ers included the most famous disease resistance proteins in plants. Some examples are UDP-glucosyltransferase on chromosome 2D, P-loop-NTPase on chromosome 2D and 6D, leucine-rich repeat on chromosome 2B and 6B and othe well-known genes like as Cytochrome P450, Alkalinephosphatase, Oxidoreductase, Patatin-like phospholipase and Haem peroxidase (Table 5).

Discussion
Novel sources for resistance to Z. tritici Zymoseptoria tritici is one of the most important foliar diseases in many wheat-growing areas, including Europe, Northern America and Asia (Hardwick et al. 2001;Mehrabi et al. 2015). The use of genetic resistance is the most appropriate strategy to control the disease. However, the rapid adaptation of Z. tritici populations leads to a quick breakdown of resistance and. Thus, continuous characterization and utilization of new sources of resistance in breeding programs are prerequisites (Abrinbana et al. 2010;Ghaneie et al. 2012). Iran is one of the primary centers of origin of wheat, and it is proposed that the co-evolution of wheat and Z. tritici occurred in this region. Therefore, the characterization of resistant wheat genotypes using isolates from this region is likely required to increase the durability of resistance (Ghaneie et al. 2012;Makhdomi et al. 2015;Aghamiri et al. 2015). We have recently studied the interactions of 185 wheat genotypes against 10 Z. tritici isolates from different sources (Mahboubi et al. 2020). Most wheat genotypes showed a susceptible pattern to all isolates. High broad-sense heritability suggested that the resistance variation is heritable, which are in agreement with previous studies on septoria resistance in different wheat germplasm in both seedling Among the isolates, IPO323 showed the highest number of incompatible interactions (n=35), while IPO02166 (originated from Iran) showed a high level of aggressiveness on wheat genotypes. Most of the wheat differentials possessing known Stb genes were susceptible to most of Z. tritici isolates, which are in agreement with previous reports for the ineffectiveness of known Stb genes against Z. tritici populations Makhdoomi et al. 2015;Mahboubi et al. 2020). Interestingly, two genotypes (ER-M-92-20 and IPK40740) and four differentials (Kavkaz-K4500, Arina, Riband and M3) were resistant to all isolates. These genotypes could be of interest as resistance sources that contain resistance genes or a combination of diverse yet-unknown Stb genes. Previous reports indicated that Stb6 was ineffective against studied isolates (Mahboubi et al. 2020). Therefore, it can be concluded that the high resistance pattern of these genotypes should be due to the presence of Stb15 or new unknown resistance genes. Besides these highly resistant genotypes, ve genotypes (IPK45227, IPK26116, IPK41079, IPK16452 and IPK40793) showed resistance to nine isolates and can be used as valuable resistance sources in wheat breeding programs. Among all genotypes, M3 showed highly resistance (immune) responses to all isolates. This genotype contains Stb16 and Stb17 (Tabib Ghaffary et al. 2012;Mahboubi et al. 2020). Therefore, it can be concluded that this gene still is effective against STB, which is consistent with previous (Hosseinnezhad et al. 2014;Makhdoomi et al. 2015). The resistant genotypes identi ed in this study are likely novel sources of resistance, which can be used in breeding programs for the development of modern wheat cultivars.

QTL validation and alignment to previously reported STB genes and QTLs
In line with our previous study, we used GWAS analysis to identify novel QTLs against Z. tritici isolates. This approach enables breeders to enhance crop genetic improvement by incorporating suitable QTLs into wheat breeding programs (Ibrahim et al. 2020). To this aim, DArTseq markers were successfully used to genotype a globally diverse wheat germplasm. The use of high-density markers with broad genome coverage in GWAS improved the accuracy of identi ed QTLs for resistance to STB, which is a highly quantitative disease trait with a minor contribution of each QTL (Mirdita et al. 2015;Muqaddasi et al. 2019). Overall, we found 37 QTLs for resistance to STB that were located on 17 chromosomes ( Table 5). The phenotypic variation explained by each QTL mapped on different chromosomes was relatively low (R 2 ≤0.14), suggesting that the resistance to this pathogen follows a highly quantitative nature, which is consistent with previous reports (Kidane et al. 2017;Arraiano and Brown, 2017;Yates et al. 2019).
The precise comparisons of these QTLs with known QTLs were di cult due to different populations, isolates and markers used elsewhere. However, using consensus wheat maps, it was possible to compare QTLs with the mapped chromosomal location of previously known genes/QTLs. Most of the QTLs identi ed in this study were localized in adjacent regions with known QTLs that have previously been identi ed (Goudemand et al. 2013;Brown et al. 2015). For example, the QTLs mapped on chromosomes 3B (Qstb.iau-14), 5A (Qstb.iau-21) and 5B (Qstb.iau-24) were mapped in close chromosomal position of Stb14 (Cowling, 2006), Stb17 (Tabib Ghaffary et al. 2012) and Stb1 (Adhikari et al. 2004b). In addition, QTLs on chromosomes 6A and 7B associated with resistance to IPO86013 may also represent Stb15 (Arraiano et al. 2007b) and Stb8/Stb13 (Adhikari et al. 2003;Cowling, 2006), respectively (Supplementary Figure 2). Other identi ed QTLs in this study on chromosomes 2B, 5A, 7A and 7D also may represent Stb9 (Chartrain et al. 2009), Stb17 (Tabib Ghaffary et al. 2012), TmStb1 (Jing et al. 2008) and Stb4/Stb5 (Adhikari et al. 2004a;Arraiano et al. 2001), respectively. Furthermore, we found several QTLs in close position of Meta-QTLs for resistance to Z. tritici isolates at both seedling and adult stages reported by Goudemand et al. (2013). Therefore, QTLs identi ed in this study can be used as valuable sources for introgression of these QTLs into advanced wheat lines (Odilbekov et al. 2019).
In addition to previously known QTLs, we found several novel QTLs. On chromosome 5D, a new QTL (Qstb.iau-26) was identi ed for resistance to IPO86013 isolate. This QTL is localized distantly from recently reported QTLs on this chromosome (Yates et al. 2019). This region on chromosome 5D is also known as the introgressed region from Aegilops umbellulata and Ae. tauschii into wheat and is an important chromosomal location having resistance genes against leaf rust, stripe rust, soil-borne mosaic virus and powdery mildew (Bansal et al. 2020;Mohler et al. 2020;Liu et al. 2020;Fu et al. 2014). Additionally, we found two novel QTLs on chromosome 4A (Qstb.iau-17 and Qstb.iau-18) and a QTL on chromosome 6D (Qstb.iau-32), which are distantly localized far from the known STB genes (Stb7/Stb12 and Stb18 mapped on chromosomes 4A and 6D, respectively). Therefore, we assume these QTLs could be novel QTLs associated with resistance to STB. In addition, we detected novel genomic regions on chromosomes 1A, 1B, 2D and 5B, in which their chromosomal positions were not overlapped with previously known QTLs/genes. Detection of these novel QTLs provided useful information that could be used to track favorable alleles for developing wheat cultivars resistant to STB. This knowledge can be used for generation of new allelic combination through cross between novel sources for resistance to STB (Riaz et al. 2020). Lastly, some QTLs identi ed in this study were co-located with known STB genes or previously mapped Meta-QTLs for resistance to STB at the seedling and adult stages in the eld (Goodwin et al. 2015;Goudemand et al. 2013). This suggests that these QTLs remain effective as durable sources for resistance to STB against isolates from diverse countries (Riaz et al. 2020).

Annotation of disease resistance genes
We characterized 29 MTAs corresponded to the putative genes functionally related to disease resistance in plants (Table 5). Intervals on chromosome 2A contained genes that were enriched for cytochrome P450, papain-like cysteine protease, zinc nger protein and chloramphenicol acetyltransferase-like protein. Cytochrome P450s (CYPs) belong to the oxidoreductases class of enzymes and their potential roles in xenobiotic metabolism are well characterized in plant species for resistance to various biotic and abiotic stresses (Yan et al. 2016;Pandian et al. 2020). It has been shown that a pathogen-induced CYP82C2 gene and other possible CYPs play important roles against the bacterial canker pathogen (Pseudomonas syringae) (Rajniak et al. 2015). Papain-like cysteine proteases (PLCPs) control key processes and are a central hub in plant immunity (Richau et al. 2012;Misas-Villamil et al. 2016). F-Box Protein CPR1 is well-known for its effect on plastid function in the absence of pathogens (Hedtmann et al. 2017). Secondary metabolism plant glycosyltransferases (UGTs) are a large family of plant genes known for their putative function during plant-pathogen interactions (Langlois-Meurinne et al. 2005). He et al. (2018) carried out a genome-wide analysis of family-1 UDP glycosyltransferases in wheat and their results showed that DP-glycosyltransferase enzymes (UGTs) were involved in detoxi cation as well as FHB resistance by glycosylating DON into DON-3-glucoside (D3G). Other intervals corresponded to some putative genes enriched for proteins containing leucine-rich repeat protein, anti-microbial extrusion protein, protein kinase and Haem peroxidase (Table 5), which are well-known genes for their signi cant roles in plants for resistance against biotic stresses (Kushalappa et al. 2016;Andersen et al. 2018;Han, 2019).

Conclusion
Resistance to STB is genetically complex and thus identi cation and comparison of QTLs for multiple isolates can be informative for designing appropriate breeding programs (Goodwin, 2007;Eriksen et al. 2003;Mirdita et al. 2015;Poppe et al. 2015;Mirzadi Gohari et al. 2015;Haueisen et al. 2019). Association mapping is a powerful genetic tool to identify chromosomal locations and QTLs related to disease resistance. This strategy has been used to detect QTLs against Z. tritici. However, the majority of these studies were conducted based on eld trials under the epidemic condition (Vagndorf et al. 2017;Arraiano and Brown 2017;Moqaddasi et al. 2019 ) or plants inoculated by limited isolates (Odilbekov et al. 2019). In this study, we used a large number of Z. tritici isolates (n=10) originated from different countries against a panel of worldwide diverse wheat genotypes to perform a deep genome association study. In addition, we used known sequence-based markers (SilicoDArT and SNP) annotated to the wheat reference genome to nd the possible function of identi ed QTLs/genes. Although, most of the QTLs identi ed in this study co-localized with previously known STB QTLs/genes, several novel genomic regions associated with resistance to multiple Z. tritici isolates were identi ed. The QTL localized on chromosome 5D confers resistance to IPO86013 isolate and QTLs on chromosomes 4A and 6D confer resistance to multiple isolates. Interestingly, based on physical mapping of QTLs against available reference genome sequence, we characterized several candidate genes involved in plant defense mechanisms against pathogens. These genes are of interest and their exact roles in STB resistance remain to be functionally analyzed in the future.
Identi cation of wheat genotypes possessing new QTLs/alleles is an important step in breeding programs towards introgression of these QTLs into new modern cultivars. Molecular and functional characterization of these QTLs/genes eventually will enhance our understanding of how resistance is achieved and sheds light on biochemical mechanisms underlying resistance against STB.

Declarations
Acknowledgements: The authors acknowledge the Iranian Seed and Plant Improvement Institute (SPII) for their kind support in germplasm preparation and diseases phenotyping experiments.
Author contributions: MM contributed in disease phenotyping and wrote the manuscript. RT and RM planned and designed the research, performed in analysis and wrote the manuscript. AMN performed the analysis. MM and GK read, edited and approved the manuscript.
Funding: Partial nancial support was received from Islamic Azad University, Sanandaj Branch, Iran.  Table 5. Summary of the septoria tritici blotch resistance quantitative trait loci identi ed against 10 Zymoseptoria tritici isolates in the panel of 185 wheat genotypes.

Figure 2
The neighbor-joining cluster analysis (a), determination of the optimal value of K and population structure analysis (b) and principal coordinate analysis (c) of 185 wheat genotypes based on DArTseq markers Figure 3 Heatmap plot of kinship matrix displaying relationships of 185 wheat genotypes based on DArTseq markers.

Figure 4
Different sub-population of 185 wheat genotypes panel that showed a different level of septoria tritici blotch disease severity (measured using the percentage of leaf area covered by pycnidia)