Genome-Wide Association Study Reveals Novel Quantitative Trait Loci Associated with Resistance to Multiple Leaf Spot Diseases of Spring Wheat

Accelerated wheat development and deployment of high-yielding, climate resilient, and disease resistant cultivars can contribute to enhanced food security and sustainable intensification. To facilitate gene discovery, we assembled an association mapping panel of 528 spring wheat landraces of diverse geographic origin for a genome-wide association study (GWAS). All accessions were genotyped using an Illumina Infinium 9K wheat single nucleotide polymorphism (SNP) chip and 4781 polymorphic SNPs were used for analysis. To identify loci underlying resistance to the major leaf spot diseases and to better understand the genomic patterns, we quantified population structure, allelic diversity, and linkage disequilibrium. Our results showed 32 loci were significantly associated with resistance to the major leaf spot diseases. Further analysis identified QTL effective against major leaf spot diseases of wheat which appeared to be novel and others that were previously identified by association analysis using Diversity Arrays Technology (DArT) and bi-parental mapping. In addition, several identified SNPs co-localized with genes that have been implicated in plant disease resistance. Future work could aim to select the putative novel loci and pyramid them in locally adapted wheat cultivars to develop broad-spectrum resistance to multiple leaf spot diseases of wheat via marker-assisted selection (MAS).


Introduction
Wheat (Triticum aestivum L. subsp. aestivum) is the major staple for more than 35% of the world's population [1]. In particular, the demand for wheat has been always high in the developing world. To meet the projected global food demand by 2050 and alleviate poverty [2,3], the pace of wheat improvement must accelerate. However, wheat production faces numerous threats, especially climatic changes and onset of severe plant disease epidemics, which significantly reduce both yield and grain quality [2,3,4,5]. Among plant diseases, bacterial leaf streak (BLS) [caused by Xanthomonas translucens pv. Undulosa [6], Tan spot [Pyreno-  [7,8,9,10,11]. These diseases can cause up to 50% yield reduction under conditions conducive to disease development [12,13].
We previously identified a set of Diversity Arrays Technology (DArT) loci associated with resistance to BLS, PTR races 1 and 5, SB, and SNB [43,44,45]. A major limitation using DArT markers is that these markers are not uniformly distributed across wheat genome. In addition, only relatively few polymorphic markers were identified that could be used for analysis. More recently, genome-wide SNP markers have been used to uncover multiple targets for wheat improvement [46]. We hypothesize that these newly developed SNPs would be associated with loci conferring resistance to multiple leaf spot diseases of wheat and could be used to validate the QTL identified previously using DArT markers [43,44,45]. Using sequence data associated with the significant SNPs, it is now possible to postulate the potential biological function related to resistance. To discover new allelic diversity and loci underlying resistance to major leaf spot diseases, and accelerate MAS, we characterized 528 diverse spring wheat accessions from the USDA-ARS National Small Grains Collection (NSGC) using 9K SNP wheat chip and GWAS analysis.

Association mapping panel
The association mapping panel utilized for GWAS consisted of 528 hexaploid spring wheat accessions. These accessions were locally-grown landraces originated from 55 countries in six continents and held by the NSGC in Aberdeen, ID. Where necessary, accessions were advanced by single plant selection. To identify QTL associated with STB resistance, we performed two independent experiments in controlled growth chambers as described previously [15]. A highly aggressive test isolate of Zymoseptoria tritici (Ma04-9-4) from North Dakota [3,4,47] was selected and inoculum was prepared as described previously [47]. Seven-week-old plants were spray-inoculated (20 ml inoculum per pot) with a hand sprayer and immediately transferred into a mist chamber with 100% relative humidity at 24uC for 72 to 96 h. Test plants then were transferred to a growth chamber programmed for a 22/18uC diurnal temperature and a 16-h photoperiod. Flag leaves were assessed 21 to 28 days post inoculation based on percentage of leaf area of necrotic lesions containing pycnidia [4]. The disease severity rating scale was from 0 to 100%, where the accessions with scores #30% severity were classified as resistant and those with scores.30% were classified as susceptible. Analysis of variance (ANOVA) was conducted using Statistical Analysis System (SAS) software version 9.3 (SAS Institute, Cary, NC). To identify novel QTL, the STB phenotypic data from this study, and the BLS, PTR races 1 and 5, SB and SNB phenotypic data from the previous studies [43,44,45] were utilized for GWAS and performed individually.

SNP marker data
Briefly, the 528 spring wheat landraces were genotyped at the Regional Genotyping Laboratory, USDA-ARS, Fargo, North Dakota using the Illumina iSelect beadchip assay for wheat having 9,000 SNPs. To avoid monomorphic and low-quality SNPs, data was sorted using Genome Studio software [46]. Nearly 5,634 informative SNPs were selected and used for GWAS. Missing data were imputed using the FastPHASE [48] with default settings. Markers with minor allele frequencies (MAF) less than 0.05 were removed from the data set in subsequent analysis, since the power of association in these alleles were low [49].

Population structure and relatedness
Population structure was analyzed using two methods. The principal components (PC) were estimated in SAS 9.3 using the Princomp procedure. The principal components were further used for GWAS. The population structure was also estimated using STRUCTURE.2.3.4 [50]. The admixture model with a burn-in was 100,000 and 500,000 iterations was used for each run. The subpopulations tested range from 1 to 15 and five runs for each K value were performed. The optimum number of subpopulations was determined by the Wilcoxon two sample test as described by Rosenberg et al. [51] and Wang et al. [52]. The Delta K approach used structure harvester [53] and the Wilcoxon test compared the posterior probabilities of two successive sub-populations (k1 vs. k2, k2 vs. k3, k3 vs. k4, and so on) using the NPAR1WAY procedure in SAS. The smaller k value in a pairwise comparison for the first non-significant Wilcoxon test was chosen as the best number of subpopulations [54,55,56]. These results were further used to interpret the geographic distribution of the landraces.

Genome-Wide association analysis
We employed four regression models: Naïve, PC, Kinship, and PC+Kinship. Among these, the Naïve model that did not account for population structure and relatedness and the regression model with only PC were analyzed in SAS 9.3. Models with only kinship and a combination of both PC and kinship were analyzed in Gemma 0.92 [33]. The number of PCs explaining 50% of the cumulative variation were used in the regression model to control for population structure and the kinship matrix estimated as a center matrix using Gemma 0.92 [57] was used to control for population relatedness. The underlying regression equation for the association mapping analysis is y = Xa+Pb+In+e where, y is a vector of phenotypic values, a is the fixed effect for the candidate marker, b is a vector of fixed effects regarding population structure, X is the vector of genotypes at the candidate marker. P is a matrix of the principal components, n is a vector of the random effects pertaining to co-ancestry; I, is an identity matrix, and e is a vector of residuals. The variances of the random effects are assumed to be Var(n) = 2KVg and Var(e) = IV R , where K is the kinship matrix that defines the degree of genetic covariance between a pair of individuals, Vg is the genetic variance and V R the residual variance [58]. Among the four models for each trait, a best model was selected based on the smallest Mean Square Difference (MSD) between the observed and expected p-values [54], since the random marker p-values follow a uniform distribution [59].

Marker-Trait associations
Association between SNPs and disease resistance traits were considered significant if the p-value was #0.001 [60,61]. To detect significant markers for each trait, the phenotypic variation (R 2 ) was calculated using a simple regression equation implemented in General Linear Model (GLM) procedure in SAS 9.3. The least squares means of the alleles of significant markers were estimated using the GLM procedures with six principal components as covariates in the model. In addition, stepwise regression implemented in the SAS REG procedure was used to estimate the combined variation explained by the markers. A significant pvalue of 0.05 was necessary for both marker and model for stepwise inclusion of the marker in REG procedure in SAS 9.3. This approach identified the major markers within a QTL excluding markers in LD of the QTL. In addition, this approach includes only the markers from major QTL masking the effects of minor QTL and has the advantage of considering the correlations and interactions between the QTL [62]. This subset explains the most phenotypic variation similar to variation explained by all markers together. A subset of markers is also more easily used for MAS compared to the entire set of markers. The adjusted consensus map developed using seven parental crosses that has 7,497 markers mapped on 21 chromosomes (represented as 25 LG) [46], were used to position the QTL on the wheat genome.

SNP marker annotations
The sequences of the significant markers available for 9K SNP wheat chip [46] were blasted against gene models of Brachypodium distachyon [63], Oryza sativa [64], and Sorghum bicolor [65] that are available at phytozome.net. The search was limited to the top hit with an E-value cut off of at least 1E-10. Further, we determined if the significant marker was in the coding or noncoding region. If the marker was in coding region, the substitution was designated as synonymous (no change in amino acid) or nonsynonymous substitution (change in amino acid).

Allelic combinations
Allelic combination refers to the combination of the marker alleles that effect the changes in phenotype. To discover allelic combinations, we employed SNPs in stepwise regression and calculated the mean and standard deviation of the phenotype. Based on the cut-off scale for resistance and susceptibility, allelic combinations were further used to identify resistant sources.

Genome-wide linkage disequilibrium
Linkage disequilibrium (LD) is the square of the correlation coefficient (r 2 ) between markers. To investigate the extent of LD across the wheat genome and the markers that have a position on the consensus map [46], r 2 between intra chromosomal SNP markers was estimated using SAS 9.3. We plotted the intrachromosomal r 2 values against the genetic distance, using a nonlinear regression in SAS [66]. The distance at which the LD decays to 0.7 was considered as the critical distance up to which a QTL region extends.

Phenotypic diversity
Association mapping panel exhibited substantial phenotypic diversity for all leaf spot diseases investigated ( Table 1). The distribution of phenotypes ranged from susceptible to resistant across the 528 accessions (Table S1). ANOVA revealed that the interactions between the two STB experiments were not significant (p#0.05), suggesting that the results of both experiments were independent. The Bartlett's chi-square (x 2 ) value was 5.1 and the associated p value with 1 degree of freedom was 0.04. Therefore, data from homogenous experiments were pooled and used for GWAS. Similar procedures were used for BLS, PTR race 1, PTR race 5, SB, and SNB [25,26,27]. Pair-wise comparison of the Pearson correlation coefficient values across diseases were significant (p,0.05) except for PTR race 1 and SNB, SNB and BLS, and SNB and SB (Table 2).

SNP marker data
Of the 5,634 SNPs obtained from the 9K SNP wheat chip, 80% had known chromosomal locations [46]. Heterozygotes accounted for 0.1% of the SNPs and were converted to missing data to estimate for sequencing errors. Approximately, 1% of the missing data was imputed using a likelihood approach. In total, 4,781 SNPs with minor allele frequencies greater than 5% were used for subsequent analysis. Of the total 4781 polymorphic markers, 94.37% (4512 SNPs) had known chromosome positions. The mapped markers were not evenly distributed across wheat genome. Our analysis revealed 45.47%, 43.79% and 5.10% of the SNPs were distributed on wheat genome A, B, and D, respectively ( Figure 1). Nearly 430 (8.99%) polymorphic SNPs were on chromosome 2B. In contrast, only 11 (0.23% of the total) polymorphic SNPs were detected on chromosome 4D.

Population structure
The Bayesian based clustering approach implemented in STRUCTURE revealed the presence of six subpopulations evaluated using the Wilcoxon test ( Figure 2). A majority of the individuals have a membership coefficient (q i ).0.7 to be assigned to a subpopulation revealing a strong population structure among individuals with little admixture [67,68]. However, these populations could not be assigned based on the geographic regions, because each of these populations identified from STRUCTURE analysis had wheat accessions from Africa, Asia, Europe and the Americas.
Based on the MSD for the four regression models tested, the regression model that has only Kinship was considered best for PTR race 1, PTR race 5 and SNB (Table 3; Figure S1). Similarly, mixed model containing PC and Kinship was considered best for BLS, SB and STB ( Figure S1; Table 3).

Marker -Trait associations and annotations
Eight SNP markers were significantly (p,0.001) associated with resistance to BLS and detected on chromosomes 1A, 2B, 3A, 5A, 5D and 6B (Table 4, Figure 3). The phenotypic variation ranged from 1.9 to 7.6% (Table 4). Of the eight significant SNPs identified, five were associated with a gene model (Table 5). Among five changes, three were non-synonymous on chromosome 5A, 5D and 6B and two changes were synonymous on chromosome 1A (Table 5). Of the eight significant SNPs identified, four fit into a stepwise regression and explained 14.3% of the phenotypic variation (Table 6). These SNPs belong to four QTL regions on chromosomes 1A, 5A, 5D, and 6B.
Eight SNPs were significantly (p,0.001) associated with resistance to PTR race 1 and located on chromosomes 2B, 4B, and 7A (Table 4; Figure 3). The phenotypic variation explained by these SNPs ranged from 1.0 to 6.5% (Table 4). The phenotypic mean difference between the alleles for the significant SNPs ranged from 0.36 to 0.67. Of the eight significant SNPs identified, three were associated with a gene model (Table 5) with two being non-synonymous (chromosome 4B and 7A) and one being a synonymous change (chromosome 2B) ( Table 5). Among the eight significant SNPs, five fit into a stepwise regression and explained Table 1. Statistical properties of major leaf spot diseases analyzed in this study. 13.8% of the phenotypic variation (Table 6). These SNPs belong to five QTL regions and were present on chromosomes 2B, 4B and 7A. For PTR race 5, six SNPs were significantly (p,0.001) associated with resistance and located on chromosomes 2A, 3A, 3B, and 6A. The phenotypic variation ranged from 2.5 to 5% (Table 4; Figure 3) while the phenotypic mean difference between the alleles for the significant SNPs ranged from 0.22 to 0.4. Among the six significant SNPs detected, four were associated with a gene model (Table 5). Of these four changes, three were nonsynonymous (chromosome 2A, 3A and 3B) and one was synonymous (chromosome 6A). For the six significant SNPs, four fit into a stepwise regression and explained 13.2% of the phenotypic variation (Table 6). These markers belong to four QTL regions and were present on chromosomes 2A, 3A, 3B, and 6A. Eleven SNPs were associated with SB resistance and detected on chromosomes 1B, 5A, 5B, 6B and 7B (Table 4, Figure 3). The phenotypic variation explained by the eleven SNPs ranged from 0.1 to 5.8% (Table 4). The location of the QTL identified by SNP markers wsnp_JD_c12281_12555386 and wsnp_Ku_c44362_51657973 (R 2 ranged from 0.2 to 0.4) could not be assigned to a chromosome because the map location are unknown. The phenotypic mean difference between the alleles for the significant markers is ranged from 0.01 to 0.4. Of the eleven significant markers identified, six were associated with a gene model (Table 5). Of these six changes, five were non-synonymous on chromosome 1B, 5A, and 5B and one was synonymous located on chromosome 5B (Table 5). Among the 11 significant markers, four fit into a stepwise regression and explained 8.1% of the phenotypic variation (Table 6). These markers were located on four QTL regions on chromosomes 1B, 5B, and 6B.
In total, eight SNPs were significantly (p,0.001) associated with SNB resistance and detected on chromosomes 2D, 3A and 5B (Table 4; Figure 3). The phenotypic variation explained ranged from 0.01 to 14.5% (Table 4). The phenotypic mean difference between the alleles for the significant markers ranged from 0.38 to 0.90. Of the eight significant SNPs detected, five were associated with a gene model (Table 5) and were non-synonymous. Of the eight significant SNPs, four fit into a stepwise regression (Table 4) and explained 28.3% of the phenotypic variation (Table 6). These SNPs belong to three QTL regions.
Seven SNPs were significantly (p,0.001) associated with resistance to STB and detected on chromosomes 3B, 6B and 7B (Table 4, Figure 3). These belong to four QTL regions based on   Table 4). The phenotypic mean difference between the alleles for the significant SNPs ranged from 14.98 to 20.95. Of the seven significant SNPs detected, six were associated with a gene model (Table 5). Of these six changes, five were non-synonymous on chromosome 3B, 6B and 7B and two changes were synonymous located on chromosome 3B (Table 5). Among seven significant SNPs, four fit into a stepwise regression and explained 19.5% of the phenotypic variation (Table 6). SNPs from four QTL regions on chromosomes 3B, 6B, and 7B were included in the stepwise regression model.

Linkage disequilibrium and allelic combinations
Based on a LD cutoff of 0.7 (correlation is 60.83), the critical LD was defined at a distance of ,4 cM ( Figure 4). Significant markers within a 4 cM interval can be defined as a single QTL. The numbers of allelic combinations for each disease ranged from 9 to 21 based on markers included in the stepwise regression (Table 7, Figure 5). Except for PTR race 5, all had at least one resistant allelic combination. Of the 528 spring wheat landraces analyzed, nearly 40% were susceptible to all leaf spot diseases and 45.4% were resistant to at least one disease (Table S3).
Importantly, four wheat landraces were resistant to at least three leaf spot diseases (STB, SNB and BLS) and 65 landraces were resistant to at least two of the diseases.

Discussion
In spring wheat, few sources of broad-spectrum resistance to major leaf spot diseases are available. Due to this limitation, tremendous efforts have been made in the past decades to identify and introduce new sources of resistance from wild tetraploid wheat, such as emmer (T. diccoccum), Persian (T. cathalicum) and Polish (T. polanicum), and other wheat-alien species derivatives [69,70,71]. Recently, it has become feasible to rapidly test for thousands of SNP markers [39,46]. In the present study, we analyzed association between disease resistance and SNPs from an association mapping panel of 528 spring wheat landraces. Our data indicate that spring wheat landraces exhibit considerable phenotypic and molecular variation, possibly due to the diverse genetic background of the accessions. We identified 32 SNPs significantly associated with loci conferring resistance to major leaf spot diseases. To validate broader applicability SNPs and GWAS, we also sought to verify resistances that were previously detected using DArT markers [43,44,45]. The higher marker density utilized in the present study enabled validation of previous findings [25,26,27]. Indeed, most of the loci detected previously by DArT markers also were identified by SNPs. Many of the SNPs significantly associated with QTL were found to be co-localized with candidate genes for plant defense and host plant resistance to several important diseases. This study provides a first step towards pyramiding resistance loci from these donors via MAS, which will enhance the genetic diversity for resistance in modern wheat germplasm and facilitate accelerated breeding to develop broadspectrum resistance to manage leaf spot diseases of spring wheat.
The high-throughput SNP genotyping array and a high-density map developed previously have enabled GWAS to identify putative QTL associated with disease resistance [46]. In contrast to bi-parental mapping, association mapping provides us an excellent opportunity to analyze a larger pool of wheat accessions to uncover QTL. Furthermore, an association mapping panel with high MAFs, low LD, and limited population structure are ideal to perform association mapping analysis [49]. To investigate genetic structure, we performed a comprehensive GWAS to discover and localize QTL in spring wheat landraces of diverse geographic origin. Our data analysis supports hypothesis that spring wheat mapping panel did not have a strongly defined population structure and in addition the panel had a low LD, thus making  Table 3. Mean square difference for four models used to identify the best regression model to discover single nucleotide polymorphisms (SNPs) and leaf spot disease resistance trait -marker associations.  it an excellent source for AM of multiple leaf spot diseases of wheat. Population structure using structure analysis suggested two to six sub-populations. Previously, population structure was assessed for the same population using DArT markers [45], which were developed from coding regions. However, the SNPs tested here are from across wheat genome with various levels of selection pressure, presumably leading to a better picture of the population structure. When population structure is present, using the structure to control for spurious associations is imperative in association mapping analysis [72]. However, the relationship between individuals varied greatly, which might be due to unique selection pressures in each of the diverse environment where the accessions originated. These factors should be adjusted appropriately in an association analysis to control false positives [49].
One of the prerequisites for GWAS is LD, which is the nonrandom association of alleles at separate loci located on the same chromosome. The marker density needed for achieving a reasonable mapping resolution is highly related to the distance at which LD declines with genetic or physical distance. The amount of LD decay also varies with different crops. For example, LD decayed within 20 to 30 cM in rice [73,74] and 10 to 40 cM in wheat [75,76] based on different samples and marker systems used. If the LD distance is too large the QTL extends to 10 cM, making it difficult to identify the significant genes within the QTL region. Such large LD distance is possible if the population has narrow genetic diversity. In the present study, the extent of LD decay was about 4 cM. In general, if the LD is low, more markers are necessary. Although we found several SNPs associated with resistance to wheat leaf spot diseases, marker coverage was not distributed evenly across the genome, suggesting that the present study may have been unable to detect QTL in the genomic regions with low marker density.  We deployed 4,781 SNPs to perform GWAS and identified 48 significant SNPs associated with resistance to major leaf spot diseases of wheat. These constituted 32 QTL and explained the phenotypic variation (R 2 ) ranging from 0.6 to 14.5%. As expected, the phenotypic variation effect was low compared to previously reported QTL detected from bi-parental mapping [9,10,14,25,77]. The large differences in the explained phenotypic variations of QTL reported in linkage versus association mapping could be due to the number of recombinant events under study. A stepwise regression was used to find the subset of markers and QTL that can have a masking effect on minor effect markers and QTL [62,78]. This procedure will limit the markers for MAS. With this approach, we were able to limit the number of loci to 25 markers and 24 QTL regions. Some of the QTL detected in this study may have already been previously identified in other studies. However, the positions cannot be related precisely due to the use of different linkage maps and markers for each of these studies.
A high level of diversity of wheat accessions can help us to better understand the genetics of resistance and to identify novel genomic regions linked with resistance genes. Due to low marker coverage of the D genome, we were unable to identify genetic regions within this genome associated with BLS, PTR 1, PTR 5, SB and STB disease resistance. In addition, none of the QTL identified were common among the different diseases.
Four QTL were identified for BLS resistance on chromosomes 1A, 5A, 5D, and 6B. Of these four QTL, those on 1A and 6B were also detected in association analysis using DArT markers [44]. The remaining QTL were mapped in novel genomic regions where no QTL were previously reported. The QTL responsible for PTR race 1 resistance were mapped to chromosomes 2B, 4B, and 7A. Some major and minor QTL regions on 2B and 7A were previously identified [14,56]. The QTL on 2B also was detected in AM analysis using DArT markers [45]. Ptr ToxB toxin sensitivity gene Tsc2 is located on chromosome 2B [14]. Similarly, another major QTL, QTs-ksu-2B, has also been mapped to chromosome 2B [79]. In addition, 2B has a QTL responsible for resistance to a novel PTR isolate [56]. Chromosome 7A, where the present study found the QTL responsible for resistance to PTR race 1, also has QTL responsible for resistance to novel PTR isolates [56]. QTL associated with resistance to PTR race 5 were identified on chromosomes 2A, 3A, 3B, and 6A. The genomic region 6A also was detected in a GWAS analysis using DArT markers [45]. Similarly, QTL identified in the genomic region 2A coincide with the host selective toxin insensitivity QTL QTs.fcu-2A [15], which conferred resistance to all known races of PTR tested. None of the other three QTL identified were mapped previously, and thus were considered novel. The present study further confirmed that the PTR-wheat pathosystem is complex and that targeting toxin insensitivity gene alone will not inevitably lead to PTR resistance [15,17].
Of the five genomic regions (chromosomes 1B, 5A, 5B, 6B, and 7B) identified for SB resistance, 1B, 5B and 7B were mapped previously [18,19,79]. The genomic region 7B also was detected in a previous study [44]. Likewise, one minor QTL with phenotypic variation effect (R 2 ) of 15.1% was detected on chromosome 1B [80] and one major QTL, Qsb.bhu-5B, was mapped on chromosome 5B [18,19]. No previous evidence was observed for resistance to SB on chromosome 5A and thus this region appears to be a novel.
Three genomic regions (2D, 3A and 5B) had major and minor QTL responsible for SNB resistance [20,22,43]. The genomic regions 2D and 5B were detected previously [43]. The QTL on 2D and 5B were previously identified with the major genes Snn2, responsible for sensitivity to SnTox2, and tsn1, responsible for Table 5. Cont. sensitivity to SnToxA [81,82]. Three QTL responsible for resistance to STB were detected on 3B, 6B, and 7B, where no QTL have been detected previously. GWAS also was able to detect SNPs associated with QTL that were identified previously with DArT markers. For example, genomic regions 1A and 6B for BLS, 2B for PTR race 1, 6A for PTR race 5, 7B for SB, and 2D and 5B were detected using both SNP and DArT markers. However, some of the QTL detected in the previous studies were not found in this study. One possible explanation was that the mapping populations used to develop the consensus DArT and SNP maps were different, thus making it difficult to compare linkage maps. However, SNP markers were able to detect additional novel QTL which were not identified by DArT markers. This result could be expected since the SNP markers were more dispersed across the wheat genome compared to DArT markers, which tend to cluster and show low density in the centromeric regions and D genome of wheat [83]. SNP markers for GWAS may be more robust and cost-effective for QTL discovery than are DArT markers. Several major genes or QTL responsible for resistance to PTR race 1, SB, and SNB were detected, which also were previously identified from conventional bi-parental mapping. GWAS can dissect the putative genes responsible for controlling phenotype [54,84] and ann in silico approach was used to probe for such genes. The SNP marker sequences were blasted against database with coding regions of Oryza sativa, Sorghum bicolor, and Brachypodium distachyon and genes associated with plant disease resistance were identified. In addition, some of the gene models identified either have no known function or may not be involved in plant defense to pathogens. The sequences of SNPs associated with QTL for resistance to BLS show similarity to sequences coding for chaperone DnaJ-domain superfamily protein, ATP-citrate lyase A-3, and MAK10 homologue. The gene that encodes the chaperone DnaJ-domain super-family protein might play a critical role in biotic and abiotic stress response [85]. This gene was over-expressed in soybean revealing its vital role in cell death and disease resistance [85]. Another gene encodes ATPcitrate lyase and it may be involved in phytoalexin formation and was up-regulated in hot pepper leaves when challenged by a pathogen [86,87].    Some QTL identified for resistance to PTR race 1 and PTR race 5 were found in the same genomic regions where known functional genes that are up-regulated in response to biotic or abiotic stress have been reported. For example, the SNP sequences that are linked to QTL responsible for PTR race 1 disease resistance are related to sequences coding for multidrug resistanceassociated protein 5 and protein kinase superfamily protein with a octicosapeptide/Phox/Bem1p domain and are non-synonymous with these genes. Similarly, sequences for SNP markers that are non-synonymous and associated with resistance to PTR race 5 in the present study are related to sequences coding for glycosyl hydrolase family 10 protein/carbohydrate-binding domain-containing protein, a heat repeat-containing protein, and oxidative stress 3. For example, haloacid dehalogenase phosphatases were found in glycation repair by direct dephosphorylation of phosphoglycated proteins or DNA or by averting the intracellular concentrations of the phosphorylated aldoses from reaching toxic levels [88]. Other important genes that encodes Ca 2+ -ATPase are directly or indirectly involved in several functions including processing of proteins in the secretary pathway, transport of Mn 2+ , and adaptation to salt stress [89]. The other gene that encodes a multidrug resistance-associated protein assisted in transporting the oxidized form of glutathione, a function essential in redox signaling activated by reactive oxygen species (ROS) in plant reactions to pathogen attack [90]. Yet another gene that encodes a protein kinase super-family protein that was associated with biotic and abiotic stress in plants [91,92,93]. A calcium dependent protein kinase has been reported to be the major component of innate immunity signaling pathways and some of the receptor-like protein kinases have been associated with plant defense responses.
The sequences from SNPs that are non-synonymous and are linked to QTL responsible for resistance to SB are related to cysteine protease 1 precursor, a NagB/RpiA/CoA transferase-like superfamily protein, a Calcium-dependent lipid-binding (CaLB domain) family protein, and an O-fucosyltransferase family protein. The roles of cysteine proteases and protease inhibitor genes in the regulation of programmed cell death in plants have been well-documented [94]. The calcium-dependent lipid-binding (CalB domain) family protein gene is concerned with transducing various stress signals to alter stress-regulated gene expression [95,96].
The QTL responsible for resistance to SNB are possibly related to alpha-N-acetylglucosaminidase family/NAGLU family, thymidylate synthase 1, DHHC-type zinc finger family protein, MYB family transcription factor and all of these are non-synonymous changes. Of the several encoded genes for SNB, the zinc finger family protein plays a major role in plant disease resistance and has been shown to be highly unregulated and responsible for early defense responses against E. amylovora infection in apple [97]. The sequences of SNP markers linked to QTL conferring resistance to STB were found to be related to sequences for different disease resistance gene groups such as a nucleotidediphospho-sugar transferases superfamily protein, a cofactor assembly of complex C cofactor assembly of complex C, and transmembrane proteins 14C. Of the genes encoding for STB, the transmembrane protein plays an important role in disease resistance [98]. In particular, the Arabidopsis NDR7 gene (contains two putative transmembrane domains) was essential for the expression of resistance to bacterial and fungal pathogens mediated by several R gene products [98]. Although we discovered several SNPs associated with novel QTL, functional analysis of the selected genes involved in host plant resistance needs further investigation.

Implications for wheat disease resistance breeding
The genome-wide analysis of SNP markers in spring wheat landraces provided a basis for comprehensive analysis of QTL resistance to the major leaf spot diseases. We discovered potentially novel QTL and further confirmed a number of major and minor QTL detected in previous association analyses using DArT markers and bi-parental mapping approaches. Resistance to each leaf spot disease of wheat appears to be controlled by relatively high numbers of QTL. Pyramiding putative resistant alleles for resistance to several diseases had been successfully utilized in various crops via MAS [99,100,101]. The spring wheat landraces used in the present study harbor multiple putative resistant alleles, which can be useful for MAS breeding. We identified 32 QTL associated with resistance to the major foliar diseases of wheat and markers identified using stepwise regression and the allelic combinations would be good candidates for further marker validation work. For example, SNP markers wsnp_Ex_c10596_17293363, wsnp_CAP11_rep_c4157_1965583, wsnp_Ex_c5998_10513766 and wsnp_Ex_rep_c67164_65655648 can be used for MAS while developing wheat cultivars resistance to BLS. Similarly, wsnp_BF473744B_Ta_2_2, wsnp_Ex_ c19772_28771627, wsnp_Ex_rep_c67159_65649966, wsnp_Ex_ c9971_16412345, and wsnp_Ex_c9971_16412270 can be used for PTR 1, while wsnp_Ex_c2887_5330426, wsnp_Ra_c44141_ 50623811, wsnp_Ex_c2920_5385184 and wsnp_Ex_rep_ c67468_66068960 can be used for PTR 5. Likewise, wsnp_Ex_c24700_33953160, wsnp_Ex_rep_c70120_69069789, wsnp_Ku_c50354_55979952, wsnp_Ku_c20701_30355248, and wsnp_Ex_c15785_24157360 can be used for SB, wsnp_ Ex_c23239_32477458, wsnp_Ku_c9269_15583444, wsnp_ CAP11_c318_261649, and wsnp_Ku_c40334_48581010 can be used for SNB, and markers wsnp_Ex_c12220_19528388, wsnp_RFL_Contig4792_5787180, wsnp_Ex_rep_c106072_ 90285324, and wsnp_JD_c646_966400 can be used for MAS while pyramiding QTL in wheat cultivars resistance to STB. One approach for validation would be to develop near-isogenic lines (NILs) in different genetic backgrounds via MAS backcrossing and evaluating them in multi-location field trials to confirm the efficacy of these QTL. Further, the broader effects of these QTL can be determined by testing NILs against multiple leaf spot pathogens. MAS breeding can be performed at an allelic level by combining several putative resistance QTL in a cultivar. In the present study, at least 15 spring wheat landraces had QTL for resistance to five of the six diseases tested. Based on the breeding target, wheat landraces from the present study could be selected as parents. For example, spring wheat landraces PI624606 and PI422235 were the most resistant accessions for majority of the leaf spot diseases tested in this study. These resistance sources could be crossed with commercial cultivars that are susceptible to various diseases. Progeny could be selected with both superior commercial traits and the markers for various disease resistant QTL. Finally, progeny with the highest number of putative resistance QTL could be further selected for testing disease resistance in multi-environments. This strategy may enable the development of cultivars with stable resistance to multiple leaf spot diseases of spring wheat. Figure S1 Comparison of QQ plots for different association models for major wheat leaf spot diseases. Observed vs. expected P values are shown for (A) Bacterial leaf streak (BLS), (B) Pyrenophora tritici-repentis race 1 (PTR race 1) (C) Pyrenophora tritici-repentis race 5 (PTR race 5), (D) Spot blotch (SB), (E) Stagonospora nodorum blotch (SNB), (F) Septoria tritici blotch (STB) using four different models with different corrections of co-founding factors (see Materials and Methods). Based on MSD for the four regression models tested, a regression model that has only Kinship was considered best for resistance to PTR race 1, PTR race 5, and SNB and mixed model containing PC and Kinship were considered best for resistance to BLS, SB and STB. (PDF) Table S1 Lists of wheat accessions along with their origin and disease reactions to multiple leaf spot diseases. R and S determine if the genotype is resistant or susceptible based on the raw score of the wheat accessions. (XLSX)