Genome-wide association mapping reveals a rich genetic architecture of stripe rust resistance loci in emmer wheat (Triticum turgidum ssp. dicoccum)

Key message SNP-based genome scanning in worldwide domesticated emmer germplasm showed high genetic diversity, rapid linkage disequilibrium decay and 51 loci for stripe rust resistance, a large proportion of which were novel. Abstract Cultivated emmer wheat (Triticum turgidum ssp. dicoccum), one of the oldest domesticated crops in the world, is a potentially rich reservoir of variation for improvement of resistance/tolerance to biotic and abiotic stresses in wheat. Resistance to stripe rust (Puccinia striiformis f. sp. tritici) in emmer wheat has been under-investigated. Here, we employed genome-wide association (GWAS) mapping with a mixed linear model to dissect effective stripe rust resistance loci in a worldwide collection of 176 cultivated emmer wheat accessions. Adult plants were tested in six environments and seedlings were evaluated with five races from the United States and one from Italy under greenhouse conditions. Five accessions were resistant across all experiments. The panel was genotyped with the wheat 90,000 Illumina iSelect single nucleotide polymorphism (SNP) array and 5106 polymorphic SNP markers with mapped positions were obtained. A high level of genetic diversity and fast linkage disequilibrium decay were observed. In total, we identified 14 loci associated with field resistance in multiple environments. Thirty-seven loci were significantly associated with all-stage (seedling) resistance and six of them were effective against multiple races. Of the 51 total loci, 29 were mapped distantly from previously reported stripe rust resistance genes or quantitative trait loci and represent newly discovered resistance loci. Our results suggest that GWAS is an effective method for characterizing genes in cultivated emmer wheat and confirm that emmer wheat is a rich source of stripe rust resistance loci that can be used for wheat improvement. Electronic supplementary material The online version of this article (doi:10.1007/s00122-017-2957-6) contains supplementary material, which is available to authorized users.


Introduction
Stripe rust (yellow rust) of wheat, caused by Puccinia striiformis f. sp. tritici (Pst), is a threat for wheat production in temperate areas throughout the world. Since the year 2000, stripe rust epidemics have been more common in the wheat growing areas previously believed to be too warm for widespread epidemic development as a consequence of the migration of temperature-adapted populations (Chen et al. 2010;Hovmøller et al. 2016). Populations of Pst with increased virulence and aggressiveness continue to spread and have defeated many major resistance genes and even combinations of these genes present in wheat varieties (Milus et al. 2009;Wan and Chen 2014). Although more than 70 officially designated yellow rust (Yr) genes and 100 temporarily named Yr genes or quantitative trait loci (QTL) (Rosewarne et al. 2013) have been identified in cultivated wheat and its wild relatives, susceptibility of widely grown cultivars is a global concern and nearly US$1 billion are lost to Pst annually (Beddow et al. 2015). Effective Yr resistance deployed in wheat varieties across most wheat growing regions of the world is needed to address the growing threat of stripe rust and minimize ongoing losses.
Cultivated emmer wheat (Triticum turgidum ssp. dicoccum), a self-pollinated allotetraploid species (2n = 4x = 28, AABB) with hulled grain, was one of the first crops domesticated in the world (Zohary and Hopf 2000). It was domesticated from wild emmer wheat (Triticum turgidum ssp. dicoccoides), which is also the progenitor of the two widely cultivated wheat species: tetraploid durum wheat (Triticum turgidum ssp. durum, 2n = 4x = 28, AABB) and hexaploid bread wheat (Triticum aestivum ssp. aestivum, 2n = 6x = 42, AABBDD) (Zeist et al. 1991;Luo et al. 2007;Faris et al. 2014). Previous studies have revealed that wild emmer wheat exhibits promising genetic diversity for resistance against Pst (Zaharieva et al. 2010) and direct hybridization between emmer wheat and cultivated wheat is highly fertile (Grama et al. 1984;Peng et al. 1999;Uauy et al. 2005), making it a promising genetic resource for enhancement of Pst resistance in cultivated wheat species. Several Yr resistance genes have been transferred from emmer wheat to cultivated wheat including Yr15 (Gerechter-Amitai et al. 1989), Yr30/Sr2 (McFadden 1930), YrH52 (Peng et al. 1999), Yr35/Lr52 (Dadkhodaie et al. 2011) and Yr36 (Uauy et al. 2005;Fu et al. 2009). The emmer wheat-derived pleiotropic adult plant resistance (APR) gene Sr2/Yr30 is arguably one of the most important resistance genes used in wheat breeding (McIntosh et al. 1995), especially given its importance in breeding for durable resistance to highly virulent 'Ug99' races of the stem rust pathogen (Singh et al. 2005). Despite the potential of emmer wheat to contribute novel alleles to wheat improvement, large collections of emmer wheat have been scarcely characterized using contemporary genomics resources and statistical genetics methods, which has hindered the exploitation of genetic variation in emmer wheat for cultivar improvement.
Most Yr genes and QTL were mapped by using bi-parental mapping populations that detect genetic variants between two, or at most several, parental accessions. Recently, linkage disequilibrium (LD)-based genome-wide association (GWAS) has attracted much interest, mostly in hexaploid bread wheat, as a method complementary to traditional bi-parental mapping to identify novel loci responsible for rust resistance in global germplasm collections. In wheat, GWAS investigation methods have been used for the first time to identify rust-resistance genome signals by Crossa et al. (2007) in hexaploid wheat and by Maccaferri et al. (2010) in tetraploid wheat. Since then, the method has been widely used to identify leaf rust (Kertho et al. 2015), stem rust (Yu et al. 2012;Letta et al. 2013;Laidò et al. 2015) and stripe rust resistance, to verify the effect of previously discovered Yr genes and QTL, and to assess the distribution of resistance loci in global germplasm (Tadesse et al. 2013;Zegeye et al. 2014;Maccaferri et al. 2015b;Naruoka et al. 2015;Bulli et al. 2016;Pasam et al. 2017).
A major advantage of GWAS approaches over bi-parental linkage mapping is that phenotype-genotype relationships in existing germplasm collections or natural populations may be explored for hundreds to thousands of accessions simultaneously, hence avoiding the resource intensive development of bi-parental genetic mapping populations. Association mapping enables researchers to use more genetic diversity and historical recombination events, which may result in higher mapping resolution than linkage mapping with a few hundred progeny (Zhu et al. 2008;Brachi et al. 2011) when high-density sequence-based markers are available (Wang et al. 2014b).
Association mapping has limitations, particularly as detection power depends heavily on the allele frequency distributions in a population. Rare alleles, with a frequency of ~5% or less, are not efficiently discovered by association mapping, unless their effects are large and can be validated through other information (Rafalski 2010). In addition, if samples in the population are not genetically homogeneous, population structure may cause spurious correlations between markers and traits resulting in elevated false-positive rates (Zhao et al. 2007). Recently, several statistical models have been developed to effectively control population structure as well as genetic relatedness among individuals in GWAS studies, including the mixed linear models (MLM) (Yu et al. 2005), compressed MLM (cMLM) (Zhang et al. 2010b), efficient mixed-model association (EMMA) (Kang et al. 2008), Settlement of MLM Under Progressively Exclusive Relationship (SUPER) model (Wang et al. 2014a) and fixed and random model circulating probability unification (FarmCPU) (Liu et al. 2016). LD decay rate, which relates to the mapping resolution of association analysis, is another important feature that determines the effectiveness of association mapping (Flint-Garcia et al. 2003). A faster LD decay rate in a population often implies a higher mapping resolution as long as a large number of genetic markers are available (Zhu et al. 2008). Compared with cultivated durum and bread wheat, emmer wheat populations may harbor more genetic variants and show a higher mapping resolution because populations have a longer history of recombination and have not been subject to the same extent of domestication bottlenecks (Salamini et al. 2002;Zaharieva et al. 2010). Sela et al. (2014) observed a rapid LD decay in a wild emmer wheat population of 128 accessions, which decayed to its background level within 1 cM. But until 1 3 now, no study has targeted the whole-genome LD pattern in cultivated emmer wheat.
In this paper, we employed a panel of worldwide cultivated spring emmer wheat accessions, representing a large portion of the genetic diversity in the cultivated spring emmer wheat gene pool. Our objectives were to assess the genetic diversity, population structure and LD pattern of cultivated emmer wheat using high-density genome-wide SNP markers, identify stripe rust resistant/susceptible accessions, while exploring the genetic architecture of stripe rust resistance loci by GWAS experiments.

Association mapping panel
A collection of 196 cultivated spring emmer wheat accessions, representing a large portion of the worldwide genetic diversity in the cultivated spring emmer wheat gene pool, was assembled from germplasm bank accessions originally collected and stored at the Institute of Biosciences and Bioresources, National Research Council of Italy (Volpe et al. 2005), Bari, Italy. The accessions were initially grown in observation rows at the University of Bologna experimental farm in Cadriano (44°33′N, 11°21′E, 32 m asl), Bologna, Italy, and then subjected to a 3-year single seed descent (SSD) purification step in 2009-2011. At the end of the third year of head-row advancement, the resulting selections (one to three for each original accession, depending on the phenotypic variability observed in the original accessions) where thus considered as genetically homogeneous and therefore amenable to subsequent genetic and phenotypic characterization. The selections were subjected to Illumina Infinium iSelect 90k wheat array genotyping and phenotyping for responses to Pst infection. After filtering for genotypic data quality and redundancy of genotypes, 176 accessions were further retained for population structure and kinship analyses (Supplemental File 1). Briefly, the genotypes can be grouped based on their geographic origins: group I: 22 accessions selected in Western European countries (Austria, Germany and United Kingdom); group II: 36 accessions derived from Southern Europe (Italy, Romania and Spain); group III: 28 accessions collected from Eastern European countries (Bulgaria, Georgia, Hungary, former Yugoslavia, Russia, Ukraine); group IV: 42 accessions cultivated in Eastern Africa (Ethiopia and Kenya); group V: one accession derived from North Africa (Morocco); group VI: five accessions selected in Southern Asia (India); group VII: two accessions bred in Eastern Asia (China); group VIII: 33 accessions collected from the Southern Levant and Western Asia (Afghanistan, Armenia, Iran, Syria and Yemen); and group IX: three accessions cultivated in Northern America (Canada and the USA) and group X: four accessions with unknown origins.  Liu et al. (2017a, b): PSTv-37 predominated in all the environments apart from CLF15, and PSTv-52 was another important race detected in all environments except from MTV15.

Stripe rust evaluation at the adult plant stages in the field
At each location, approximately 150 seeds of each accession were planted in a 0.5-m row, spaced 0.3 m apart without replicates. The position of the 196 accessions was random in different field trails. To insure the uniformity of Pst inoculum across the trial, spring wheat cultivar 'Avocet Susceptible' (AvS) was planted every 20 rows as spreader rows, and another susceptible spring wheat variety, 'Lemhi', was planted around each plot as borders. Before planting, standard practices for fertilization and weed control common to the region were applied for field management. Two disease phenotypes were evaluated: infection type (IT) based on a 0-9 numerical scale and disease severity (SEV) rated as a percentage of disease area covering flag leaves (0-100%) (Line and Qayoum 1992). The stripe rust response in six environments was evaluated three times after the susceptible check 'AvS' had IT from 7 to 9 and SEV from 70 to 100%. The highest IT and SEV were used as the final phenotypes for each accession. Heading date (HD) and plant height (PH) traits for each accession were also recorded. Since the evaluation of Pst resistance may be influenced by flowering time, the accessions with very late HD were excluded when analyzing field nursery response to Pst races.

Seedling resistance evaluation
Five Pst races  that are most common in US Pst populations, and one Pst race (PSTv-125) collected in Italy (Liu et al. 2017a), were used for seedling resistance screening of the emmer wheat panel under controlled conditions (Supplemental Table 1). PSTv-37 is widely distributed across the US 1 3 (Wan and Chen 2014), PSTv-14 and PSTv-40 are common in the western US and PSTv-51 had the broadest spectrum of virulence. PSTv-18 is the least virulent race identified in the US so far, which is avirulent on all 18 Yr single-gene differentials. PSTv-125 was the most common race collected from both hexaploid and tetraploid wheat in Italy in 2014.
Seedling resistance to Pst was tested in controlled greenhouse conditions. Three seeds of each accession and susceptible check 'AvS' were planted in 96-well trays. Six Pst races (Supplemental Table 1), representing various virulence patterns, were chosen to test seedlings of this panel for allstage resistance. Seedlings were grown under a diurnal temperature cycle of 20-24 °C and 15-18 °C at night with 16 h photoperiod for 10-14 days. At the two-leaf stage, plants were inoculated in a 4 °C dark dew chamber with 100% relative humidity for 24 h and then placed in a growth chamber. A core set of Pst differentials containing 18 Yr genes (Wan and Chen 2014) was included in each test to confirm the race identity. Infection type was recorded 16-20 days after inoculation. For each accession, three plants were recorded as one IT score if their reactions were the same. If reactions were different, we re-tested that accession. In addition, we re-tested the accessions with IT scores of 0-2 in order to validate their resistance reactions to the corresponding Pst races.

Molecular profiling
Genomic DNA was extracted from a bulk of 25 one-weekold seedlings per accession using the DNeasy 96 Plant Kit (Qiagen GmbH, Hilden, Germany). The 196 accessions were genotyped using the Illumina ® iSelect 90K wheat SNP assay (Wang et al. 2014b) at the USDA-ARS Biosciences Research Laboratory, Fargo, ND. SNP marker genotype calling was performed using the GenomeStudio v2011.1 software package (Illumina, San Diego, CA, USA). The dataset filtering was carried out based on the following criteria: (1) markers showing residual heterozygosity were entered as missing values; (2) markers with less than 10% missing data and accessions with less than 20% missing data were retained; (3) markers with minor allele frequency (MAF) greater than 10% were retained; (4) one representative accession was retained for groups of accessions with genetic similarity equal to 1. A total of 176 lines and 5106 polymorphic SNPs with known map positions based on a high-density consensus map of tetraploid wheat (Maccaferri et al. 2015a) were used for population structure and kinship analyses. Molecular markers barc8 (Yaniv et al. 2015) and wMAS000005 (Distelfeld et al. 2006) that are widely used as markers for Yr15 and Yr30/Sr2, and KASP_IWA6121 and KASP_IWA4096 (Naruoka et al. 2016) that flank Yr5, were profiled on this emmer wheat panel (Supplemental Table 2).

Statistical analysis of phenotypic data
The broad sense heritability (H 2 ) of IT and SEV were calculated by a random effect model using PROC MIXED COV-TEST in SAS v.9.3 (SAS Institute Inc., Cary, NC, USA) treating genotype, environment and genotype × environment interaction as random factors. The equation of H 2 is as follows (Piepho and Möhring 2007): where σ G 2 is the variance component of genotypes; σ 2 E × G is the variance of the interaction of environment and genotype; σ e 2 is the variance of residual; n is the number of environments and r is the number of replicates per environment.
Using the PROC ANOVA in SAS v.9.3, a one-way ANOVA and Tukey's mean comparison tests were performed for IT and SEV among STRUCTURE-based subpopulations at a significance level of P < 0.05, where subpopulation was treated as fixed effect.
PROC CORR in SAS v.9.3 was used to compute the Pearson's correlation coefficients. Correlations of IT and SEV among environments were computed to investigate the consistency of Pst responses across different environments.
Best linear unbiased predictors (BLUP) for each accession across environments was computed for traits including field-based Pst responses (IT and SEV), HD and PH, using PROC MIXED procedure in SAS v.9.3 (Wolfinger et al. 1997). In the PROC MIXED procedure, phenotypes of check 'AvS' were considered to be a fixed effect while environment and the interaction of environment and genotype were random effects.

Genetic diversity and LD analyses
POWERMARKER v3.25 (Liu and Muse 2005) was used to calculate two genetic diversity parameters: genetic diversity and polymorphism information content (PIC). Genetic diversity is the probability of two randomly chosen alleles from the population being different (Weir and Cockerham 1996). PIC estimates the detection power and informativeness of the SNP markers (Botstein et al. 1980). The formula is where p i is the frequency of the ith allele of the jth marker, n is the number of allele at the jth marker.
LD was measured as the squared correlation coefficient (r 2 ) between SNP pairs in the same chromosome using JMP ® Genomics 6.0 (SAS Institute Inc., Cary, NC, 2012). All 5106 mapped SNP markers were used in this LD analysis. To visualize the local LD patterns and overall LD decay rate, the intra-chromosomal pairwise r 2 estimates were plotted against inter-marker genetic distance (cM) (Fang et al. 2013;Rexroad and Vallejo 2009) and a non-parametric regression curve was fitted with smoothing spline, lambda = 10,000 1 3 function in JMP ® Genomics 6.0. The confidence interval of QTL was defined based on the intersection of the smoothing spline curve with the pairwise LD baseline (r 2 = 0.3) (Maccaferri et al. 2015b;Bulli et al. 2016).

Association analyses
To minimize false-positive marker-trait associations (MTAs), population structure and genetic relatedness were factored into the GWAS analyses. Only SNP markers with known map position were loaded into HAPLOVIEW v4.2 (Barrett et al. 2005) to choose representative SNPs (haplotype tag-SNP). A total of 3361 tag-SNPs were chosen for calculating the marker-based kinship matrix (K matrix) based on the tagger function set at r 2 = 1.0, and 1422 tag-SNPs were further selected for calculating population structure (Q matrix) analysis using a setting of r 2 = 0.5. The population structure for this panel was determined by the Bayesian model-based clustering method in the software STRUCTURE v.2.3.4 (Pritchard et al. 2000;Falush et al. 2003). The parameter settings of STRUCTURE were: admixture model of population structure with correlated allele frequencies, 50,000 length of burn-in period, 100,000 number of Markov Chain Monte Carlo (MCMC) replications after burn-in, hypothetical subpopulations k setting from 1 to 10 with five iterations. STRUCTURE HARVESTER (Earl 2012) collated the output that was generated by STRUC-TURE. An ad hoc statistic ∆k, on the basis of the rate of change in the logarithm of likelihood [ln P(D)] between successive k values, was used to detect the number of groups that best represented this population (Evanno et al. 2005).
Kinship coefficients of all pairs of accessions were calculated based on the identity-by-state (IBS) genetic similarity (Kang et al. 2008). Association mapping was conducted in the R package Genome Association and Prediction Integrated Tool (GAPIT) (Lipka et al. 2011) using all 5106 mapped polymorphic SNPs suitable for GWAS (MAF >0.10). Different association models were tested separately for the entire panel for seedling resistance and the subset panel for field resistance: (1) the general linear model (GLM) corrected for population structure using STRUCTURE membership (thereafter, Q GLM); (2) MLM (Yu et al. 2005) incorporating kinship as the random component (thereafter, K MLM); (3) MLM with Q matrix as the fixed factor and K matrix as the random factor (thereafter, Q+K MLM). Model selection was made based on the Bayesian information criterion values (BIC) computed by GAPIT (Supplemental Table 3). The K MLM model was chosen for the seedling response MTA analyses since it showed the best BIC values for the majority of Pst races (analysis based on all 176 accessions selected for GWAS). For the field-based adult plant responses, GWAS was conducted on a further subset of 109 accessions selected for narrow HD (14 day range). Thus Q and K matrices were re-calculated based on the genotypic data of 109 accessions. Model comparison was conducted and the K MLM model was chosen due to the best BIC values across the majority of environments (Supplemental Table 3). To control the impact of HD and PH on the level of stripe rust resistance in the field, HD and PH were included as covariates in the K MLM model for each of the six environments and the BLUP data across these six environments. Liu et al. (2017a) characterized 93 loci associated with seedling and APR to Pst races in a worldwide collection of elite durum wheat. Straightforward map-based comparison of significant SNPs detected in this study with SNPs detected in elite durum wheat panel was performed by using a highdensity SNP based on the tetraploid wheat consensus map (Maccaferri et al. 2015a). Maccaferri et al. (2015b) constructed an integrated map including 225 previously reported Yr genes and QTL, SNPs, simple sequence repeat (SSR), diversity arrays technology (DArT) and express sequence tag (EST) markers. The map distances in the integrated map are based on the hexaploid wheat consensus map reported by Wang et al. (2014b), which is different from the tetraploid wheat consensus map (Maccaferri et al. 2015a) used in this study. The map positions of identified QTL in the tetraploid wheat consensus map and integrated map were reported in Supplemental Table 5, and the map positions of identified QTL in the integrated map were utilized to compare with the positions of previously published Yr genes in the integrated map, based on common shared markers. The use of alignment results is approximate and caution should be exercised because of the inherent limitations in the consensus map and integrated map.

SNP marker quality, genomic distribution and genetic diversity
In total, 196 accessions of cultivated emmer wheat were genotyped using the 90K (81,587 actual markers) Illumina Infinium iSelect wheat SNP array. A total of 30,559 polymorphic SNPs with high-quality genotype calling were generated, of which 20,730 SNPs were mapped to the cultivated emmer wheat genome based on a highdensity consensus map of tetraploid wheat (Maccaferri et al. 2015a). After filtering the markers with the criteria of MAF > 0.10 and missing data frequency <0.10, 5106 polymorphic SNP markers with known map positions were retained. Removing accessions with more than 20% missing marker frequency and genetic similarity equal to one, 176 genotypes were used for the following analyses. Marker distribution, MAF, number of alleles of each SNP marker, genetic diversity and PIC of each chromosome were calculated and presented in Supplemental Table 4.
In this population of cultivated emmer wheat, the B genome (3009 SNPs) had more polymorphic markers than the A genome (2097 SNPs), but differences in genetic diversity (B genome = 0.2976; A genome = 0.2802) and PIC (B genome = 0.2374; A genome = 0.2251) values were not statistically significant based on paired t tests. Polymorphic marker distributions ranged from 243 to 564 markers per chromosome, with an average density of 365 markers per chromosome. Chromosome 1B and 2B had the largest number of markers, and higher genetic diversity and PIC values than the genome-wide averages, while chromosome 3A had the smallest number of markers, and lower genetic diversity and PIC values. The overall genetic diversity and PIC values for this cultivated emmer wheat panel were 0.2895 and 0.2321, respectively.
To understand genetic diversity based on area of origin, the proportion of polymorphic loci, genetic diversity and PIC were calculated for each of the 10 origin groups (Table 1)

Population structure, kinship and LD
Population structure was analyzed by Bayesian-based structure analysis using 1422 non-redundant markers (LD r 2 = 0.5) across the whole genome. A clear peak at k = 3 revealed that the optimal population structure for further analyses was to divide this panel into three main subgroups. Subpopulation one, two and three contained 45, 78 and 53 accessions, respectively. The majority of accessions (71%) collected from Africa (Ethiopia, Kenya and Morocco) were grouped into subpopulation one ("African-subpopulation"), while 92% of accessions from Southern Europe and 77% of accessions from Western Europe were grouped into subpopulation two ("European-subpopulation"). All accessions from Eastern Asia and 82% of accessions from Western Asia contributed to subpopulation three ("Asian-subpopulation").
Genetic relatedness was also calculated using the IBS coefficient among all pairs of individuals. The IBS allelesharing matrix ("kinship matrix") clearly mirrored subpopulation ancestry as distinguished by the population structure analysis (Fig. 1). Only two accessions in STRUCTUREbased subpopulation two were clustered into IBS-based subpopulation three. Accessions in each subpopulation were highly related, and the mean IBS coefficients within the three subpopulations were 0.86, 0.70 and 0.82, for subpopulations one, two and three, respectively.
The mean genome-wide LD pattern of the 176 cultivated emmer wheat genotypes is graphically displayed by scatter and box plots of pairwise LD r 2 over genetic distance in cM (Fig. 2). The box-plot of LD r 2 distribution (Fig. 2b) showed that the median LD r 2 of the completely linked marker pairs was 0.91 with interquartile ranging from 1.00 to 0.31. In the next marker group (tightly linked markers separated <1.0 cM), the median LD r 2 dramatically decreased to 0.20, indicating an average of 78.02% LD decay rate within a 1.0 cM genetic distance. For pairs of markers with genetic distances of 1.0-5.0 cM, the median LD r 2 reduced to 0.11. Linkage disequilibrium r 2 started to decrease slowly after pairwise marker genetic distances greater than 5.0 cM. The average LD r 2 dropped from 0.1 to 0.06 when genetic distances of pairwise markers increased from 5.0 to >50.0 cM. Similar LD decay trends were determined when fitting the smoothing spline curve into the scatter-plot of LD r 2 distribution (Fig. 2a). The baseline (r 2 = 0.3) intersection with the smoothing spline curve was at 2.1 cM, which was used to determine the QTL coverage regions with marker distance confidence intervals of ±2.1 cM from the peak of the significant associations.

Seedling response to six Pst races
Seedling responses showed highly significant differences (P < 0.0001) among accessions within each race. The IT distributions of all Pst races summarized in Fig. 3 were skewed toward susceptible scores (IT = 7-9), except for the most avirulent race PSTv-18. For races PSTv-14, PSTv-37, PSTv-40, PSTv-51 and PSTv-125, susceptible reactions (IT = 7-9) were observed in 51, 81, 55, 61 and 66% accessions of this emmer wheat panel, respectively. A higher proportion (60%) of accessions exhibited highly resistant reactions (IT = 0-2) to PSTv-18 than the other races, while intermediate (IT = 3-6) and susceptible reactions to PSTv-18 were 17 and 23%, respectively. Five accessions in subpopulation one were resistant to all six Pst races, and 14 accessions in subpopulation two were susceptible to all races (Supplemental File 1). The effect of population structure was significant for responses to Pst races. Accessions in subpopulation one (African germplasm) are more resistant than those in subpopulation two (European germplasm) and three (Asian germplasm) in responses to all Pst races tested in this study except for PSTv-18 (Table 2). The mean IT of accessions in subpopulation two and three are significantly (P < 0.05) higher than that of accessions in subpopulation one. The Red to white to blue lines indicate reactions to Pst changed from susceptibility to intermediate to resistance to stripe rust. Black lines indicate the accessions that displayed winter habit or had very late heading date that excluded from the GWAS analysis at the adult stages. The order of individuals in a-c were arranged according to their IBS-based genotypic distance accessions in subpopulation two were more susceptible to PSTv-18 than accessions in subpopulation one and three. The mean IT of subpopulation two was significantly (P < 0.05) different from subpopulation one and three, but the mean IT of subpopulations one and three did not show a significant difference.

Field response to Pst in six field environments
Accessions in this cultivated emmer wheat panel showed substantial diversity in morphological and physiological traits under field conditions, including substantial variations in HD (a 32-day range). Eleven accessions in subpopulation two displayed winter habit in our field screening nurseries, and 55 accessions (53 accessions in subpopulation two) had very late HD, associated to IT and SEV = 0. In order to reduce the influence of late HD on the analysis and interpretation of field stripe rust response evaluation, 109 accessions were selected from the original collection of 176 accessions based on HD consistently comprised within a 14-day range.
The IT and SEV of 109 accessions at the adult stage were presented in Supplemental File 1. The distributions of IT and SEV of the six environments were all skewed toward resistance, and the cultivated emmer wheat population exhibited a higher level of resistance in CLF15 and WHT14 compared to the other four environments (Fig. 4). Based The heritability values for IT and SEV, calculated across the six environments, were 0.921 and 0.919. The high heritability indicated limited environmental variation compared to genotypic variation across the six environments. Pearson's correlation coefficients of IT and SEV (Table 3) recorded in the six environments were all highly significant (P < 0.0001) with the range from 0.60 to 0.77 for IT and 0.57 to 0.82 for SEV. The correlation coefficients between IT and SEV within the same environment were moderate to high, ranging from 0.50 to 0.91. The five accessions that were resistant to all six Pst races at the seedling stages were also resistant based on their evaluation over six field environments at the adult stages.

Association analyses for race-specific resistance
A mixed linear model with K as a covariate was determined to be the optimal model based on the BIC values. The ITs of 176 accessions in response to races PSTv-14, PSTv-18, PSTv-37, PSTv-40, PSTv-51 and PSTv-125 were used to identify race-specific resistance MTAs. A total of 57 SNP markers within 37 distinct loci on all 14 chromosomes were identified to be significantly associated (P < 0.001) with seedling resistance. The phenotypic variance (R 2 ) explained by each of these MTAs ranged from 3.2 to 7.5% (Table 4).
The second group included resistance loci significantly associated with only one Pst race. YrTtd-6BL and YrTtd-7AL.2 showed significant association with US race  and YrTtd-7BL were significantly associated with resistance to US race PSTv-37. Ten loci were detected to be associated with resistance to US race . Three loci, YrTtd-3AL, YrTtd-3BL.2 and YrTtd-5AL.1, showed significant association with resistance to US race PSTv-51. PSTv-18, the least virulent Pst race, was used to identify the Yr resistance loci that were effective against relatively avirulent Pst races but have been defeated by more virulent Pst races. Two loci, YrTtd-4BS and YrTtd-5BL.2, were detected in this population conferring resistance to PSTv-18. Two loci, YrTtd-2BS and YrTtd-2BL.2, were exclusively associated with resistance to the virulent race of Italian origin, PSTv-125.

Association analyses for resistance in field environments
To identify MTAs across field environments effective against contemporary US Pst populations at the adult plant stage, association tests were carried out separately for each of six field environments and also using BLUPs generated from data combined across the six environments. Significant (P < 0.05) Pearson's correlations were observed between IT and HD (r = 0.24), IT and PH (r = 0.21), SEV and HD (r = 0.45), and SEV and PH (r = 0.41). Therefore, HD and PH were incorporated into MLM + K model as two covariates.
In total, 18 loci represented by 24 SNPs were significantly associated (P < 0.001) with IT and/or SEV in the field nurseries. The phenotypic variance explained by individual loci ranged from 7.2 to 12.5%. Without considering the interaction effects, the cumulative phenotypic variation explained by all 18 loci was 56.7% for BLUP_IT and 66.1% for BLUP_SEV. Four of these loci were also significantly associated with resistance in the seedling tests. These all-stage resistance loci were YrTtd-4BL.1, YrTtd-4BL.2, YrTtd-7AS and YrTtd-7AL.2 highlighted in Table 4. The cumulative phenotypic variance explained by these loci was 10.1% for BLUP_IT and 11.2% for BLUP_SEV. Among these seedling resistance loci effective under field conditions, YrTtd-4BL.1, YrTtd-4BL.2 and YrTtd-7AS were significantly associated with resistance in multiple field environments. The remaining 14 loci on chromosomes 1A, 1B, 2AS, 2BS, 3BS, 6BS, 7AL and 7BL were significantly associated with adult plant responses in the field, but not seedling responses, and may represent APR loci (Table 5). QYrTtd-7AL.1 tagged by IWB25121 and QYrTtd-7AL.2 tagged by IWA501 are mapped to chromosome 7AL with a genetic distance of 3.6 cM. Haplotype analysis showed these two QTL are not in LD (r 2 = 0.11), and QYrTtd-7AL.1 is associated with SEV while QYrTtd-7AL.2 is associated with IT, indicating that they are likely distinct QTL.

Distributions of Pst resistance alleles in subpopulations and their influences on Pst reaction
The seven resistance loci that are effective against multiple races and/or consistently effective for field resistance were selected for exploring the frequencies of resistance-associated alleles present in STRUCTURE-based subpopulations ( Table 6). The African subpopulation had higher frequencies of associated SNP alleles for YrTtd-2BL.1, YrTtd-4AL.1, YrTtd-4BL.1, YrTtd-6AS.1 and YrTtd-7AS. The favorable alleles for the majority of resistance loci detected in the seedling tests in subpopulation two (European materials) and three (Asian materials) had lower frequencies than that of subpopulation one. Exceptions are YrTtd-7AL.2 and YrTtd-4BL.2, whose frequencies of favorable SNP-associated alleles are much higher in subpopulation two and subpopulation three, respectively. The lower frequencies of resistance-associated alleles in accessions of subpopulation two and three might explain their lower levels of resistance to Pst races PSTv-14, PSTv-37, PSTv-40, PSTv-51 and PSTv-125 at the seedling tests (Table 2).

SNP genome properties and genetic diversity revealed by 90K SNP array
The rapid advancement in next-generation sequencing in the last few years has enabled the development of genome-wide SNP arrays for the polyploidy wheat genome. In this study, the currently available 90K SNP wheat array (Wang et al. 2014b) was used to genotype the cultivated emmer wheat accessions. In spite of ascertainment bias when using the 90K SNP array on cultivated emmer wheat and the presence of D-genome markers, we were still able to detect 7756 polymorphic SNPs after data quality filtering steps, of which 5106 were positioned on a tetraploid wheat consensus map. A total genome coverage of 2542 cM, with an average marker density of 2.0 SNP/cM and an average inter-marker distance of 0.5 cM/SNP were found in this population. The 4.2 cM confidence interval of loci contained over eight polymorphic SNPs suitable for GWAS, on average. Therefore, the wheat 90K SNP array provided sufficient polymorphic markers for this cultivated emmer wheat population to carry out genetic diversity, population structure, LD and association mapping analyses.
Although many populations have been assembled to study the genetic diversity of different species of wheat (Autrique et al. 1996;Barcaccia et al. 2002;Maccaferri et al. 2003;Zhang et al. 2010a), this is the first using the high-density 90K SNP array on cultivated emmer wheat. The accessions in this study were sampled from nearly all major emmer wheat growing regions, and comprehensively sampled gene pools of spring emmer wheat, as recently    Badaeva et al. (2015). Therefore, information on genetic diversity and LD of this population can offer relevant guidelines for assembling association mapping populations and designing new breeding strategies for exploiting and utilizing cultivated emmer wheat genetic resources to improve durum and bread wheat cultivars. The genome-wide 90K SNP markers revealed an average genetic diversity index of 0.2895, PIC of 0.2321 and MAF of Table 4 Putative resistance genes that were significantly (P < 0.001) associated with responses to races   0.2168 for this cultivated emmer wheat population. These values are slightly elevated compared to previously reported estimates of cultivated durum and bread wheat populations also using genome-wide SNP markers. Other researchers reported a genetic diversity value of 0.2280 and PIC of 0.1888 among a worldwide germplasm collection of 150 durum wheat genotypes (Ren et al. 2013b), and PIC of 0.20 among 478 US and CIMMYT spring and winter bread wheat cultivars (Chao et al. 2010). The higher genetic diversity of cultivated emmer wheat may be attributed to its longer cultivation in a large range of Table 5 Putative quantitative trait loci (QTL) significantly (P < 0.001) associated with response to stripe rust in field environments  a Putative Yr loci that have significant association with seedling response to multiple Pst races are given in bold. Yr loci that have significant association with seedling and field responses are marked by asterisks b Chromosome positions of putative Yr loci are based on the tetraploid wheat consensus map (Maccaferri et al. 2015a) c Confidence intervals are determined by confidence intervals of ± 2.1 cM from the peak of the significant associations d Tag-SNPs are the most significant SNP in the confidence interval of putative genes e Resistance-associated alleles of tag-SNPs are highlighted in bold f Significant SNPs associated to the same Pst races as the tag-SNP and falling into the confidence intervals of putative genes g Phenotypic variations explained by the tag-SNPs h Pst races that SNP markers are significantly (P < 0.001) associated with 1 3 eco-geographical conditions, a longer history of recombination and less domestication bottlenecks than the later-domesticated durum wheat and bread wheat (Zaharieva et al. 2010;Badaeva et al. 2015). Lower genetic diversity (0.1841) and PIC (0.1530) values were also reported in a collection of 200 wild emmer wheat lines from Israel and Turkey (Ren et al. 2013a), which may be explained by a more narrow range of geographical distribution of this wild emmer wheat panel than our cultivated emmer wheat panel sampled globally.

Linkage disequilibrium in cultivated emmer wheat
The effectiveness of GWAS experiments heavily depends on LD decay rates regulated by the genetic distance between loci, the number of generations since it arose (Mackay and Powell 2007) and, consequently, the density and genomewide coverage of the molecular polymorphisms used in the MTA tests. Based on 5106 SNP markers that yielded 1892,402 intra-chromosomal SNP pairs to calculate the LD r 2 values, the genome-wide r 2 values abruptly declined to 0.2 within 1 cM in this cultivated emmer wheat panel. Sela et al. (2014) observed a similar LD decay rate in a wild emmer wheat population. When compared with cultivated durum wheat and bread wheat, LD decayed more rapidly in this cultivated emmer wheat collection. Previous studies reported that the genome-wide LD r 2 decayed to 0.2 within a distance of 5-20 cM in the elite bread wheat populations (Chao et al. 2007;Naruoka et al. 2015;Zegeye et al. 2014;Zhang et al. 2010a). Long-range LD was found at locus pairs with an intra-chromosomal distance more than 50 cM using 70 SSR markers in durum wheat (Maccaferri et al. 2005). The faster LD decay rate observed in this panel indicates a higher mapping resolution is possible when conducting GWAS on cultivated emmer wheat than durum and bread wheat.

Map-based comparison of the significant all-stage resistance loci with previously published Yr genes
In this emmer wheat panel, a total of 37 loci were significantly associated with all-stage resistance to Pst races detected in the seedling tests. YrTtd-1AS,  and YrTtd-7BL co-locate with resistance loci characterized from a GWAS study of elite durum wheat (Liu et al. 2017a) based on the tetraploid consensus map (Supplemental Table 5).
Only YrTtd-4AL.1, YrTtd-5BL.1 and YrTtd-6BL were associated with the same Pst races as these resistance loci tagged by IWB31333, IWB48863 and IWB55752, respectively, which were detected in the elite durum wheat panel. Haplotype analyses indicate that YrTtd-6BL and IWB55752 are likely to be the same loci because they are in high level of LD (r 2 = 0.8), but YrTtd-4AL.1 and IWB31333 on chromosome 4AL, YrTtd-5BL.1 and IWB48863 on chromosome 5BL are most likely different genes since the LD r 2 values between these two pairs of markers are less than 0.3. Further map-based comparison of significant all-stage resistance loci with previously reported Yr genes were conducted by using the integrated genetic map (Maccaferri et al. 2015b, Supplemental Table 5). On chromosome arm 1BS, YrTtd-1BS mapped within the gene region of Yr64 and very close to Yr65. Yr64 and Yr65 were resistant to all the tested races including PSTv-40,  and   (Cheng et al. 2014) used in the present study. However, YrTtd-1BS was significantly associated with resistance to PSTv-40. This result illustrated that YrTtd-1BS is likely neither Yr64 nor Yr65, and indicated this region of chromosome 1B harbors several alleles or a cluster of Yr genes (Cheng et al. 2014).
On chromosome 4A, YrTtd-4AL.1 resided in the confidence interval of resistance gene Yr60 (Herrera-Foessel et al. 2015). YrTtd-4AL.1 was significantly associated with response to race PSTv-37, PSTv-40 and PSTv-51 in the current panel, whereas Yr60 conferred a moderate level of resistance to Pst isolate Mex96.11. It is possible that YrTtd-4AL.1 represents the same resistance gene as Yr60, but seedling evaluation for both donor plants with the same Pst races and/or allelism testing are needed to further investigate their relationships. Moreover, Maccaferri et al. (2015b) reported that IWA1034 was significantly associated with APR in multiple environments, but was not detected by GWAS analysis of resistance to races  in a worldwide collection of hexaploid spring bread wheat. Since two APR QTL Prins et al. 2011) have already been identified in this chromosome region, IWA1034 may be closely linked with both all-stage resistance genes and APR genes.
In the proximal region of chromosome arm 4BL, YrTtd-4BL.1 was significantly associated with responses to PSTv-37 and PSTv-125 in the seedling tests and field responses in MTV14 and SPM15 in this cultivated emmer panel, but not associated with responses to PSTv-40. YrTtd-4BL.1 overlapped with Yr50  and Yr62 ). However, YrTtd-4BL.1 is most likely not Yr62; Yr62 is an APR gene. Yr50 confers race-specific all-stage resistance, putatively derived from Thinopyrum intermedium in Chinese winter wheat cultivar 'CH223'. Lu et al. (2014) conducted the seedling test on 'CH223' and confirmed Yr50 was highly resistant to  and . Given the wild relative origin and the different reactions of YrTtd-4BL.1 and Yr50 to the same Pst races, it is unlikely that YrTtd-4BL.1 and Yr50 represent the same gene.
The other 22 significant resistance loci detected in the seedling tests in this GWAS study appeared to be novel since they are located on the chromosome regions where no race-specific all-stage resistance genes have been mapped before, apart from IWB63861 (YrTtd-6AS.1) because map position for this marker is unknown in the hexaploid integrated map.
YrTtd-4BS and YrTtd-5BL.2 were only significantly associated with the least virulent race PSTv-18. Therefore, YrTtd-4BS and YrTtd-5BL.2 should be considered relic loci with limited effectiveness against contemporary Pst populations. Four resistance loci  and YrTtd-7AS) on chromosomes 2BS, 2BL, 4BL and 7AS were significantly associated with resistance with Italian race PSTv-125. Among them, YrTtd-2BS and YrTtd-2BL.3 were exclusively associated with resistance to PSTv-125. On the other hand, YrTtd-4BL.2 and YrTtd-7AS were more broadly effective than other resistance loci detected in this study because they were detected with three or more Pst races in the seedling tests and in fields of multiple environments. Hence, YrTtd-4BL.2 and YrTtd-7AS are effective all-stage resistance genes to be introgressed into adapted durum and bread wheat accessions.
In addition, all-stage resistance genes Yr5 and Yr15 are known to be effective against all Pst races identified so far in the United States (Murphy et al. 2009;Wan and Chen 2014) and have been incorporated into numerous bread and durum wheat cultivars (Peng et al. 2000;Sun et al. 1997). In this cultivated emmer wheat population, no significant associations were detected close to the putative positions of Yr5 and Yr15. In addition, two flanking markers KASP_IWA6121 and KASP_IWA4096 (Naruoka et al. 2016) and a tightly linked SSR marker barc8 (Yaniv et al. 2015) were used to screen the panel for the presence/absence of Yr5 and Yr15, respectively; none of the accessions carried the expected marker alleles for Yr5 or Yr15.

Field resistance loci and their relations to previously reported stripe rust QTL
It is not surprising that four all-stage resistance loci (YrTtd-4BL.1, YrTtd-4BL.2, YrTtd-7AS and YrTtd-7AL.2) conferred resistance to stripe rust in the field because Pst race PSTv-37 that was used in the seedling test prevailed among five field screening nurseries except CLF15, and YrTtd-4BL.1, YrTtd-4BL.2 and YrTtd-7AS were detected to be associated with resistance to PSTv-37.
A total of 14 loci that were significantly associated with adult plant response in the field but not seedling response were characterized in this emmer wheat population, but only QYrTtd-3BS.1 showed consistent associations with field resistance in multiple screening nurseries at the significance level of P < 0.001. When applying a less stringent significance level of P < 0.005, QYrTtd-1BL, QYrTtd-2AL, QYrTtd-2BS.2, QYrTtd-3BS.2 and QYrTtd-6BS.1 were significantly associated with resistance in at least two environments. The other eight loci that showed significant associations within single environments should be deployed cautiously since they are unstable in different environments and susceptible to some Pst races.
Comparing the APR loci putatively identified in this study with those identified in an elite durum wheat GWAS panel (Liu et al. 2017a), QYrTtd-2AL (tagged by IWB4635) mapped within the confidence interval of QYrdurum-2AL (tagged by IWA1008). However, most of the SNP markers in the haplotype of QYrdurum-2AL were not polymorphic in the emmer wheat panel which made it difficult to establish the relationship between QYrTtd-2AL and QYrdurum-2AL.
Further map-based comparisons of significant APR loci in this study with known APR loci were performed on the basis of the consensus map generated by Maccaferri et al. (2015b). Seven additional resistance loci in this emmer wheat panel overlapped with known resistance QTL (Supplemental Table 5). QYrTtd-1A, spanning the centromere of chromosome 1A, maps to the confidence interval of APR QTL QYr.sun-1A. QYr.sun-1A was derived from the Australia cultivar 'Janz' (Bariana et al. 2010), and had racespecific resistance to Pst race 110 E143A+ but susceptibility to 134 E16A + . QYrTtd-1A is likely race-specific since it is identified at a significance level of P < 0.005 in only one environment (CLF15). Allelism testing would be needed to determine the relationship of QYrTtd-1A and QYr.sun-1A.
Previously reported APR QTL, QYr.cau-1BS, mapped to the confidence interval of QYrTtd-1BS on chromosome arm 1BS. QYr.cau-1BS was characterized from the wheat breeding line 'AQ24788-83' that is the progeny of a double cross consisting of four stripe rust-resistant Chinese winter wheat landraces (Quan et al. 2013). However, QYr.cau-1BS contributed to prolonged latency for the appearance of uredinia but not reduced IT and SEV, hence, QYr.cau-1BS and QYrTtd-1BS may not be the same QTL.
On the short arm of chromosome 2B, two loci detected in the field tests in this GWAS overlapped with known APR QTL. QYr.inra-2BS derived from 'Renan' (Dedryver et al. 2009) and QYrst.orr-2BS.1 derived from 'Stephens' (Vazquez et al. 2012) mapped in the same region where QYrTtd-2BS.1 was identified in the current study. 'Stephens' was a soft white winter cultivar grown in the Pacific Northwest, USA and has a combination of all-stage and APR resistance. Stephens conferred resistance in the PNW region for more than 30 years and in recent years this resistance has been mostly been explained by APR loci. Five APR loci have been detected from 'Stephens' by Vazquez et al. (2012). QYrst.orr-2BS.1, peaked at wPt-5738, was a minor QTL that was unstable across environments. Similarly, QYrTtd-2BS.1 in the current study was only detected in CLF15 but not the other environments. QYrTtd-2BS.2 falls into the confidence intervals of six previously reported APR regions including QYrlo.wpg-2BS_Louise (Carter et al. 2009), QYrid.ui-2B.2_IDO444 , QYr.caas-2BS_Pingyuan 50 (Lan et al. 2010), QYr-2B_Opata 85 (Boukhatem et al. 2002), QYr.tam-2BL_TAM111 (Basnet et al. 2014) and QYr. ucw-2B_UC1110 (Lowe et al. 2011). Additionally, three race-specific all-stage resistant genes YrKK , Yr41 (Luo et al. 2008) andYr27 (McDonald et al. 2004) were also clustered here. This region is very close to the centromere of chromosome 2B, where recombination frequencies are greatly reduced. Therefore, allelism tests for evaluation of the complex relationships among these Yr loci are difficult in this chromosome region.
The emmer wheat-derived pleiotropic APR locus Sr2/ Yr30 is located on the distal part of chromosome 3BS. No MTAs in this emmer wheat panel were identified on the genetic region of Sr2/Yr30. Profiling with the diagnostic marker wMAS000005 of Sr2/Yr30, only four accessions ("MG5287/1" from Iran, "MG5297/1" from Bulgaria, "MG5442" from the United States and "MG15518" from Syria) were detected to carry this particular allele in this panel.
A haplotype of five SNPs (IWB124, IWB5813, IWB25636, IWB32812 and IWB50708) was located near the centromeric region of chromosome arm 3BS and was designated as QYrTtd-3BS.2 in our study. Results of haplotype analysis showed a high level of LD among these SNPs with pairwise r 2 equal to 0.90. Based on the integrated map, the confidence interval of QYrTtd-3BS.2 overlaps with the region of QYrco.wpg-3BS.2 which represents a minor APR locus detected from the soft white winter wheat cultivar 'Brundage' (Case et al. 2014).
On the long arm of chromosome 7A, QYrTtd-7AL.2 resided within the QTL region of QYr.cim-7AL that was identified from a stripe rust susceptible cultivar 'Avocet' by Rosewarne et al. (2012). QYr.cim-7AL, flanked by marker wPt-2260 and wPt-2501, is a minor APR locus in cultivar 'Avocet', that showed high susceptibility to stripe rust in the field. But QYrTtd-7AL.2 identified in our study is a major effect QTL that explained 9.1% of phenotypic variation for IT in SPM15.
QYrTtd-7BL (tagged by IWB72249) detected on chromosome arm 7BL in this study overlapped with the marker interval between wgp5175 and barc32, designated as Yr59 by Zhou et al. (2014). Yr59 was characterized from the Iranian spring wheat landrace 'PI 178759' and confers HTAP resistance. Allelism testing would be needed to verify whether QYrTtd-7BL is Yr59 or not.

Sources of resistance and susceptibility to Pst races
Breeding for Pst resistance in cultivated durum and bread wheat continues to be impeded by the lack of resistance sources effective against contemporary Pst races in elite genetic backgrounds. The relatives of these wheat species are promising resistance sources to provide various resistance genes (Uauy et al. 2005;Yin et al. 2006;Marais et al. 2009;Liu et al. 2013). In the current study, we evaluated the reactions of 176 cultivated emmer wheat accessions to four contemporary Pst races of the United States  and one Italian race (PSTv-125) at the seedling stage. About 9, 10, 14, 4 and 8% of accessions were resistant to respectively. Among them, accessions "IDG392", "MG5346/1", "MG5362", "MG5530" from Ethiopia and "MG5312" from Afghanistan were resistant to all races in seedling tests and displayed moderate to high levels of resistance in adult plants across six environments. Accordingly, these accessions may be good donor parents for wheat breeding. Conversely, 14 accessions from Italy, the United Kingdom, Spain, Romania, Germany, Yugoslavia and Iran (Supplemental File 1), were identified as ideally susceptible parents for geneticists to construct bi-parental mapping populations. These accessions were highly susceptible to all tested races at the seedling stage. At the adult stages, 13 of them exhibited winter habit or very late heading dates. The remaining one accession ("MG29292" from former Yugoslavia) was susceptible in field nurseries and would be a good susceptible parent for future studies.
In this cultivated emmer population, four accessions ("MG5287/1" from Iran, "MG5297/1" from Bulgaria, "MG5442" from the United States, and "MG15518" from Syria) were detected to carry the resistance allele of Sr2/Yr30 as well as resistance-associated alleles for QYrTtd-1A, QYrTtd-2AL,  and QYrTtd-7AL.1. These accessions show moderate to high resistance at the adult stages but susceptibility to most tested races at the seedling stage. Hence, they may serve as donor parents of APR for tetraploid wheat breeding.

Conclusion
Cultivated emmer wheat, as an important genetic resource for durum and bread wheat improvement, has attracted the attention of wheat researchers in recent years. Our study reports the presence of valuable genetic variation for multiple Pst races in seedling tests and in field nurseries. The cultivated emmer wheat accessions that showed all-stage resistance or APR could be crossed with adapted durum and bread cultivars to enhance stripe rust resistance, particularly those accessions carrying improved resistance alleles at novel loci. High-density, whole-genome SNP markers revealed a high level of genetic diversity and relatively fast LD decay, indicating that cultivated emmer wheat is suitable for GWAS studies aimed at directly identifying MTAs closely linked to the causal loci. In the GWAS analyses, 37 loci significantly associated with all-stage resistance and 14 loci that were associated with field resistance were found. More than half of these were mapped to chromosome regions where no stripe rust resistance genes have been reported before, indicating that a large reservoir of useful alleles currently underexploited in the modern breeding germplasm is present in domesticated emmer. This result highlights the presence of novel Pst resistance loci in cultivated emmer wheat that could be pyramided into durum and bread wheat cultivars by recurrent selection or marker-assisted selection. Finally, the relatively large number and small phenotypic effects of loci detected in this cultivated emmer wheat population demonstrated the complex genetic basis of stripe rust resistance in wheat, suggesting that genomic selection might be a more comprehensive strategy towards accumulation of beneficial alleles in breeding and pre-breeding populations than conventional marker-assisted selection.
Author contribution statement WL managed and performed all the experiments, analyzed the data, interpreted the results, and wrote the manuscript. MP supervised the experiments, advised data analyses and revised the manuscript. MM advised data analyses and revised the manuscript. XC, RT, DP and GL provided materials and revised the manuscript.