Plant Material

Nine Kenyan spring wheat cultivars (‘Kenya Fahari’, ‘Gem’, ‘Kenya Kudu’, ‘Kulungu’, ‘Kenya Ngiri’, ‘Kenya Paka’, ‘Pasa’, ‘Kenya Popo’, and ‘Romany’) that were released in the period 1964–1989 were selected for crossing based on the high level of APR exhibited by these lines during screening in the Njoro stem rust nursery in Kenya (Table 1). ‘Ada’, a recent hard red spring wheat variety released by the University of Minnesota [32], also exhibits moderate resistance to the Ug99 race group and a high level of resistance against N. American Pgt races (Table 1). These ten lines were crossed to the stem rust susceptible line ‘LMPG-6’ to develop RIL populations (F₆ or more inbred) via the single seed descent method at the University of Minnesota (S1 Fig). The number of inbred lines in the ten populations ranged from 55 to 110 (Table 1). Hereafter, the lines Kenya Fahari, Kenya Kudu, Kenya Ngiri, Kenya Paka, and Kenya Popo are referred to as Fahari, Kudu, Ngiri, Paka, and Popo, respectively.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Table 1. Origin, pedigree, and stem rust reaction of parent lines used to develop the NAM population.

https://doi.org/10.1371/journal.pone.0155760.t001

Stem Rust Evaluation

The ten RIL populations, comprised of 852 lines, were evaluated for APR to stem rust in four environments: St. Paul, MN, USA during May–August 2012 & 2013 (referred as StP12 and StP13 in the text); South Africa during October 2012–January 2013 (referred as SA12); and Njoro, Kenya during June–October 2013 (referred as Ken13). The stem rust data collected on all populations in all four environments are available in S1 File.

In the Njoro nursery, lines were planted in an augmented design with the susceptible check ‘Red Bobs’. Each line was sown in double rows 70 cm long and 20 cm apart. On each side and in the middle of the plots, a twin-row of susceptible spreader wheat cultivar ‘Cacuke’ was sown. The field was also surrounded by a border of several spreader rows comprised of susceptible wheat varieties that were artificially inoculated using a bulk inoculum of Pgt urediniospores collected at the Njoro field site. The wheat stem rust differential lines with known stem rust resistance genes indicated that the predominant, if not only, race present in the nursery was TTKST (avirulence/virulence formula on the wheat stem rust differential panel: Sr36, SrTmp/Sr5, Sr6, Sr7b, Sr8a, Sr9a, Sr9b, Sr9d, Sr9e, Sr9g, Sr10, Sr11, Sr17, Sr21, Sr24, Sr30, Sr31, Sr38, SrMcN) [33].

In the St. Paul 2012 environment, lines were planted in 2 m long single rows with 20 cm between the rows. The populations were planted in an augmented design with 4 check varieties ‘Oklee’ [34], ‘Thatcher’ [35], ‘Tom’ [36], and ‘Verde’ [37] planted after every 30 entries. In the St. Paul 2013 environment, lines were planted in hill-plots with 20 cm distance between the hills. The same checks were planted in the same manner as in the 2012 season. In both St. Paul environments, the lines were surrounded by a mixture of susceptible lines ‘Morocco’, ‘Max’, and ‘Little Club’ planted perpendicular to the lines on all sides. To initiate disease, spreader rows were syringe-injected at the jointing stage with a mixture of N. American stem rust races MCCFC (isolate 59KS19), QFCSC (isolate 03ND76C), QTHJC (isolate 75ND717C), RCRSC (isolate 77ND82A), RKQQC (isolate 99KS76A), and TPMKC (isolate 74MN1409). The spreader rows were sprayed with a bulked mixture of Pgt races suspended in a light mineral oil suspension using an Ulva+ sprayer (Micron Sprayers Ltd., Bromyard, UK) after heading stage.

In South Africa, populations were planted in Cedara in KwaZulu-Natal Province in an augmented design with the check lines ‘Morocco’ and ‘Kariega’ planted every 50 entries. Seeds for each line were planted in hill plots with 30 cm between plots in a row and 60 cm between rows. To initiate the disease, spreader rows containing ‘Morocco’ and ‘McNair’ were inoculated with a mixture of the Pgt races TTKSF and PTKST of the Ug99 race group using an ultra-low volume sprayer twice in the season—once during the booting stage, and again at flowering stage.

Field reaction of the RILs to stem rust were recorded as disease severity on the 0 to 100 modified Cobb scale [38], and infection response, based on the size of pustules and amount of chlorosis and necrosis visible on the stem [39]. Disease phenotyping of the population segregating for resistance was carried out after the susceptible check varieties in each trial had attained maximum severity. Following Stubbs et al. [40], the severity response value was multiplied by the infection response to obtain coefficient of infection (CI) values. The number of days to heading (in StP12 and StP13) and plant growth stage (in Ken13) were also recorded. The number of days to heading was measured as the day after planting when half the spikes in the plot fully emerged above the flag leaf. Plant growth stages were determined mainly by assessing grain development stages such as watery, milky, soft dough, and hard dough; and also for stages of booting and flowering, as explained by Zadoks et al. [41]. Additionally, all 11 parental lines were inoculated with race TTKSK (isolate 04KEN156/04) at the seedling stage to postulate the presence/absence of major genes providing resistance to this race. Disease inoculation and phenotyping procedures were carried out as described by Rouse et al. [42].

Statistical Analysis

To improve normality of the phenotypic data, arcsine-square root, log, and square root transformation methods were tested. The CI values were passed through the formulae [arcsine (√ 1-(CI+1)/100)], [log(CI+1)], and [(√ (CI+1)/100)] for arcsine-square root, log, and square root transformation, respectively. Each transformed dataset was tested for improved normality using the Shapiro-Wilk normality test. While the test confirmed that the transformed data improved in normality (Fig 1), significant departures from normality were still observed (S1 File). Based on the normality test, the square root transformed datasets were used in the next steps. Using genotypes as fixed effects, mixed linear models were then fitted to assess family and line-within-family effects, and also the effects due to differences in heading date (St. Paul data) and plant growth stages (Kenya data). Factors explaining significant amounts of variation were retained in the model and the effect estimates for each line was done using lme4 package (version 1.0–6) in R 3.0.3 (R Development Core Team 2013, http://www.r-project.org/).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 1. Distributions for disease coefficient of infection (CI) and their respective transformed datasets for stem rust in each of the four environments.

The Pearson’s correlations are represented as follows: “r1” between the coefficient of infection (CI) values and square root transformed data; and “r2” between the square root transformed data and data adjusted for trial differences. “W” represents the Shapiro-Wilk test statistic between the coefficient of infection (CI) values and square root transformed data.

https://doi.org/10.1371/journal.pone.0155760.g001

The phenotypic scale used in South Africa was different than in other locations as the disease was evaluated on a scale of 0 to 10 with scores: 0–2 (highly resistant), 3 (resistant), 4–5 (moderate resistant), 6–7 (moderate susceptible), 8 (susceptible), and 9–10 (highly susceptible) [43]. Therefore, all data sets (fixed effects estimated using lme4 for each environment) were fitted into a mixed model in SAS 9.1 (SAS Institute Inc, Cary, NC, USA) with environments as fixed effects and lines as random effects. Thus estimated best linear unbiased predictors (BLUPs) were used in QTL mapping. Heritability across all environments was calculated on an entry mean basis [44].

Genotyping and SNP Discovery

All lines used in the study were genotyped using the genotyping by sequencing (GBS) approach [28]. DNA was extracted from ground leaf tissue of parents and F_6:7 RILs using the BioSprint 96 DNA Plant Kit (QIAGEN, Valencia, CA). Extracted DNA was quantified using the picogreen assay and diluted to 20 ng/μl. The restriction enzymes PstI and MspI were used to generate double-digested complexity-reduced libraries following Poland et al. [29] with few modifications: 1) each sample was ligated to two unique barcodes to minimize sequencing bias of reads with certain barcodes; and 2) the concentrations of both barcode and common adapters were increased to 0.1 μM and 50 μM, respectively, in an attempt to capture most of the digested DNA fragments. Ten 96-plex libraries were generated, with each parent repeated at least six times to obtain higher read coverage of parental alleles. Each library was sequenced in one lane of Illumina HiSeq 2000, generating 100 bp single-end sequences. Sequences of all 11 parents and RILs of all 10 populations have been uploaded to NCBI’s sequence read archive (SRA: SRP057693; BioProject: PRJNA281776).

For each population, sequences of each RIL and the two parents were aligned to the Chinese Spring genome assembly [45] ordered into chromosomal pseudomolecules using population sequencing (POPSEQ) data generated by Chapman et al. [46]. Sequence alignment was carried out using the aln method in ‘bwa’ using default parameters. The software ‘Samtools’ was used to process the aligned sequences (view method) followed by SNP calling (mpileup method). SNPs were accepted if the alignment read depth was ≥ 3 and the minimum mapping quality was ≥ 25. Heterozygote calls were converted to missing allele calls to allow for the imputation algorithm to discern either parental genotype. The workflow of genotyping, map construction, data analysis, and genetic mapping in this study is presented in S1 Fig.

Map Construction and Joint QTL Mapping

Genetic maps for individual populations as well as the combined panel were constructed using IciMapping 4.1 [47] with a minimum logarithm of odds (LOD) value of 5.0. Genetic distances between the markers were calculated based on the Kosambi mapping function [48]. For each population, an initial map was constructed comprising only of SNP markers with proportion of missing alleles at ≤ 20% and minor allele frequency of ≥ 30%. To this, markers with missing data of up to 50% were added and the linkage groups were reordered. To construct a joint map combining all populations, SNP markers present in at least two RIL populations were used. Construction of joint map was carried out as explained above. Imputation of the missing alleles was not carried out prior to map construction and the native imputation algorithm within IciMapping was implemented during QTL mapping to reconstruct the missing haplotypes. This algorithm assigns parental genotypes to the missing markers by considering linkage relationship between markers based on their map order on each chromosome (or linkage group) within each family [49]. Genotypic data on the joint panel—both non-imputed and imputed—has been provided as supplemental information (S2–S7 Files). A joint QTL mapping across all populations was done using the joint inclusive composite interval mapping (JICIM) method in IciMapping 4.1. JICIM, an efficient and specialty method for joint QTL mapping for the NAM design, offers higher mapping power and can therefore detect small-effect QTL better [23, 50]. Within JICIM, presence of QTL along the chromosomes was scanned at an interval of 5 cM and a QTL was declared significant if the threshold was greater than the 1,000 permutation of the trait data by resampling method [51], at type I error α = 0.05. The thresholds for StP12, StP13, SA12, and Ken13 were estimated at 4.3, 4.5, 5.8, and 4.8, respectively. Epistatic effects among the detected QTL were estimated using the joint connected model in MCQTL 5.2.6 [52] by providing the QTL information, yet a lower LOD threshold of 3.0 was chosen. QTL mapping in each population was also carried out separately using the inclusive composite interval mapping (ICIM) approach in IciMapping 4.1 using LOD threshold of 3.0. For each detected QTL, the percent of phenotypic variance explained (R²) and allelic effects were also estimated.

Genotyping

We generated ten 96-plex GBS libraries representing 852 RILs from ten biparental populations. The sequencing of each library on one lane of Illumina HiSeq 2000 generated a total of 1.5 billion 100 bp reads, with 154 million reads on average per lane. On average, 90% of the generated bases passed the Q30 filter with a median Q-score of 34.98. The reads were then filtered for having intact barcode sequences and a complete PstI overhang, which led to 73% of total reads assigned to each individual. This is comparable to recent GBS studies in wheat [30] and barley [31] that follow the same library construction and sequencing protocols. This resulted in a read distribution per individual from 175,443 to 31,381,071 with a median of 1,225,681 reads per RIL. Each parent line was sequenced at least six times to obtain sequences with higher read depth for confidence in SNP-calling and imputation. As a result, a total of 158,182,462 reads were obtained for the 11 parents used in the study, with a median value of 13,348,322 reads, and ranging from 7,951,649 (Kudu) to 31,381,071 of the common parent LMPG-6 which was replicated 10 times (S1 File).

SNP Discovery, Linkage Mapping, and Segregation Distortion

After calling SNPs within each population using the Samtools program, SNPs were renamed in the format ‘Chr_position’ to reflect the base position on the reference chromosome. SNPs within each population were filtered to obtain SNPs with: a) nucleotide base as an allele for both parents and b) two parents were polymorphic for the SNP call. The amount of missing data in the resulting files ranged from 0% to 100%, with SNP counts ranging from 29,651 (LMPG-6/Kudu) to 41,418 (LMPG-6/Popo) (S1 File). Common SNPs between at least two populations, identified as the ones with the same SNP name as described above, were selected to create the joint linkage map. This led to identification of 13,413 SNPs that were used for construction of linkage maps. Of these, 11,221 SNPs were assigned to 21 linkage groups with all 21 wheat chromosome represented. The number of markers per linkage group ranged from 138 (Chromosome 3D) to 1,218 (Chromosome 3B). The total genetic distance covered by these groups was 1,920 cM, with one SNP marker placed at every 0.17 cM on average. Most SNPs were assigned to the B genome with 5,259 SNPs, followed by the A genome (3,871) and the D genome (2,091) (S1 File). Though the D genome is usually considered to be the least polymorphic of all wheat genomes [53, 54], our approach of reference-based SNP-calling may have helped to greatly increase the marker coverage (1 marker every 0.18 cM, on average) for this genome. It is also worth noting that the GBS approach has been shown to provide a better coverage of the D genome by earlier studies [29, 55].

As we used markers present in two or more populations to create the combined map, an accurate estimation of segregation distortion was affected by a large amount of missing alleles. The same held true for individual population maps where markers with up to 50% missing data were retained. While imputation of the alleles partly aided in giving better distortion ratio, it was not able to solve the problem completely. Therefore, we report here the segregation distortions in each family based on the combined map when the missing data per family was ≤ 50%. While this still does not present a clear picture, all populations can be compared by the virtue of the same map. Several markers in most chromosomes in each RIL population deviated significantly (p <0.001), yet no strong distorted pattern was visible at any particular locus on any chromosome in any population (Fig 2). The population LMPG-6/Ada had the lowest amount of segregation distortion at 32.0%, whereas LMPG-6/Romany had the highest (58.5%). The average amount of distortion across all RIL populations was 44.7%. These numbers are slightly higher in comparison to the numbers reported a recently published study that used GBS to create a consensus map (albeit in different populations than ours with much less missing data) and study several disease traits in wheat [56]. We believe that using markers with less or no missing data can more accurately illustrate the deviations from the expected segregation ratio of 1:1 distortions, especially for markers at higher significance levels (p <<0.001).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 2. Segregation distortion of loci across each RIL mapping population.

Chromosome names and the—log(p) value for SNP markers in respective chromosomes for each population is shown.

https://doi.org/10.1371/journal.pone.0155760.g002

An advantage of using a joint linkage map that combines several populations is that more recombination events are obtained per chromosome. This was evident in our joint panel as higher recombination frequencies were observed on each chromosome, relative to that in the individual populations (Fig 3). While an increased capture of recombination events leads to higher resolution in QTL mapping, it can also potentially unlock hidden genomic regions, thereby allowing for detection of otherwise unmapped QTL.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 3. Number of recombination events per chromosome in the joint map (gray bar) relative to the average number of recombinations per chromosome in all 10 populations.

https://doi.org/10.1371/journal.pone.0155760.g003

Population Characteristics and Stem Rust Reaction

Population structure among the populations was determined by singular value decomposition of the genotype data using the princomp package in R 3.0.3 (R Development Core Team 2013, http://www.r-project.org/). No significant structuring among the lines was observed (S2 Fig), which is expected of the NAM population design [21]. The lack of structure among the lines can also be credited to the parents, which are diverse from each other (S3 Fig). The lack of population structure among the populations despite having diverse founder lines is an attractive aspect of the NAM design. As different founder lines are crossed to a common parent, shuffling of the parental genomes in the progeny normalizes allelic differences by virtue of the common parent. Therefore, population stratification is effectively controlled in the NAM design, which minimizes spurious associations that could arise from population structure. Because of the lack of structure, relationship matrices were not used as cofactors during joint mapping.

All ten RIL populations and the parent lines were phenotyped over two years in four environments where disease epidemics were artificially established by inoculating spreader lines with local Pgt races. Screening for resistance to N. American and African Pgt races showed that all male parents used in the cross were resistant to moderately resistant (Table 1). Some susceptible pustules were observed on the stems of ‘Ada’ in the Kenya 2013 nursery, albeit at low severity (average of 11%). The female parent ‘LMPG-6’ succumbed to high disease severity in all four environments and exhibited susceptible infection responses. Phenotyping of all populations was initiated once disease severity on ‘LMPG-6’ was at its maximum at each environment. The frequency of resistant lines was higher in the Kenya 2013 environment than any other (Fig 1, S4 Fig), most likely due to the unusually cooler day and night temperatures at this site during the 2013 main season, which led to slower onset of disease than usual. Adjustment of phenotypic data was able to resolve this, as trait means among all environments were similar post-adjustment. The broad sense heritability across four environments was estimated at 0.55.

Seedling screening of the parent lines with the race TTKSK (Ug99) showed that all parents were susceptible to this race (Table 1). Seedling susceptibility of the parents to TTKSK indicated that the lines lacked seedling resistance genes effective against this race. The parent lines were screened for the presence of Sr2 using available DNA markers [57], and also for the presence of another APR gene, Sr57, using marker csLV34 of the linked leaf rust APR gene Lr34 [58]. Results indicated that four of the ten founder lines—‘Gem’, ‘Kudu’, ‘Paka’, and ‘Romany’—could contain Sr2, whereas none of the parents appeared to have Sr57. Marker screening of other APR genes, namely Sr55, Sr56, and Sr58, was not possible due to the unavailability of diagnostic markers.

Stem Rust QTL Mapping

The joint inclusive composite interval mapping (JICIM) approach discovered 59 QTL on 20 of 21 wheat chromosomes (Fig 4, S1 File). Though, the ‘effective’ number of QTL could be 48 if adjacent QTL (highlighted in green in Fig 4) on a few chromosomes are considered as a single QTL because of their proximity to each other. Most QTL (26) were detected on B genome chromosomes. The number of QTL detected in A and D genomes were 21 (on all A genome chromosomes) and 12 (on all D genome chromosomes except 6D), respectively. The detection of fewer number of QTL in D genome was quite interesting, and may be attributed to a lower coverage of this genome relative to the other two genomes. Despite the D genome having a generous number of markers mapped to its chromosomes, the genomic distance covered by these markers was only 50% relative to that of the A and B genomes. Fifteen of the 59 QTL were detected in more than one environment: QSr.umn-1B.2 and QSr.umn-5A.2 were common between Ken13 and StP13 environments; QSr.umn-2A.4, QSr.umn-2B.3, QSr.umn-2B.8, QSr.umn-3B.1, QSr.umn-3B.3, QSr.umn-3B.4, and QSr.umn-5B between StP12 and StP13; QSr.umn-2B.2 between Ken13 and SA12; QSr.umn-3B.5 between SA12 and StP12; QSr.umn-6B.2 between StP12 and Ken13. The QTL QSr.umn-5D.1, QSr.umn-6B.4, and QSr.umn-7B were common in three environments: StP12, StP13, and Ken13; No QTL was common among all four environments.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 4. Chromosomes with APR QTL to stem rust detected by the joint inclusive composite interval mapping (JICIM) method.

Multiple QTL in green color on a given chromosome are hypothesized to be a single QTL.

https://doi.org/10.1371/journal.pone.0155760.g004

The detected QTL explained a wide range of phenotypic variation (R²), with the lowest R² of 1.4% and the highest of 20.4%. These values indicate the presence of small to medium-effect QTL contributing towards stem rust resistance. However, small population sizes can inflate phenotypic estimation variation [59]. As some RIL populations in our study are small (smallest population with n = 52), the R² values as well as the additive effect estimates could have greater error vs. populations of larger size. The additive effect estimates for each QTL were of small magnitude (S1 File), a function of the scale of transformed data. For each QTL, the parental alleles contributed different effects towards stem rust resistance, as exhibited by the distribution of total additive phenotypic effects of each parent (Fig 5). Overall, each parent contributed towards resistance in each environment. The positive values associated with the parent lines for specific SNP markers in certain environments do not imply that the parents are not important sources of resistance. As the effects of all parental alleles are considered simultaneously during joint mapping, variations in disease response in each family, contributed by each parent, are expected to occur. Moreover, the choice of a parent for introgression of alleles for APR is also a function of several important factors defining the target environment such as frequency, and types of Pgt races, and disease pressure. While the contribution to disease resistance from alleles with low additive values might be difficult to visually observe in the field, APR genes acting in an additive manner should elevate the resistance, thereby assisting in phenotypic selection. Likewise, combining QTL that were common in multiple environments is also likely to lead towards significant reduction in disease damage.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 5. Heat map of additive effect estimates of alleles contributed by the 10 founder lines at the QTL for resistance to Pgt races.

QTL (columns) are named according to McIntosh et al. [89] with their chromosomal positions after the underscore (_) symbol. The allelic effect estimates for each founder allele (rows) are color coded by increments in the allelic effect estimate (legend). Each block represents the environments where the QTL were detected, as labeled.

https://doi.org/10.1371/journal.pone.0155760.g005

No epistasis was detected among QTL within the populations, neither in a joint setting with all populations combined. We did however observe few non-QTL markers on a few chromosomes with LOD scores in the range of 2.5–2.8. Therefore, the main reason for no significant epistasis is likely due to the stringent threshold (LOD value ≥ 3.0) required for declaration of such interaction. Recent studies have mixed results for detection of epistatic interactions between the mapped QTL as some studies report significant levels of interaction [60–63] whereas other studies do not [64–66]. Taking together the larger number of detected QTL and the small magnitude of the marker effects, the nature of resistance to stem rust in these ten populations could be assumed to be polygenic, with several loci acting in additive fashion. A point to consider is that the lack of epistasis should not deter a breeder from adoption of a QTL as long as it contributes to elevated APR levels.

Comparison of Joint Mapping to Single-population QTL Mapping

The ICIM approach of QTL detection conducted on each RIL population detected a total of 59 QTL on 20 chromosomes (no QTL detected on 3D) (S1 File). The populations LMPG-6/Gem and LMPG-6/Romany had the most QTL (13), and the population LMPG-6/Ngiri had the least (1); no QTL was observed in LMPG-6/Popo. There were no common QTL among all four environments in any individual RIL population, although four populations shared at least one QTL in three of the four environments. This difference is likely due to the presence of different rust races in the four environments as well as other environmental factors, as evident by the observed genotype by environment (GxE) interactions (S5 Fig).

We were mainly interested in knowing if the QTL detected during joint mapping could also be detected by conducting QTL mapping of each population separately. While only three QTL that were detected by JICIM were detected by ICIM based on flanking markers, a closer look at the QTL positions revealed that as many as 14 QTL on 11 chromosomes could be common between the two methods (Table 2). Nine of these 14 QTL are located on the same position, albeit QTL mapping reported different flanking markers based on the difference in the confidence interval; remaining 4 are within a distance of ≤ 3 cM. Six of the 14 QTL were detected in the same environments between the two methods, and eight were detected in different environments. The difference of the environments where the QTL were detected suggests that these loci might be involved in providing resistance that is specific to each disease environment, and perhaps to specific Pgt races. It could also indicate differences in environmental conditions, inoculum load, and difference in genetic backgrounds (e.g. different alleles). S1 Text expounds on this comparison.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Table 2. Common quantitative trait loci (QTL) detected between the joint inclusive composite interval mapping (JICIM) in nested association mapping of 10 RIL populations, and inclusive composite interval mapping (ICIM) methods in individual populations for stem rust adult plant resistance in four environments.

https://doi.org/10.1371/journal.pone.0155760.t002

Comparison of Joint Mapping Results to Previously Reported Genes and QTL Conferring Resistance to Pgt

Kenyan lines are known to be the sources of many Sr genes discovered to date [3]. Examples include Sr6, which is quite common in Kenyan lines [67], and Sr9b, which was first observed in Kenyan lines [68]. The line ‘Frontana’, present in the pedigree of the founder parent ‘Gem’, is known as one of the sources of Sr9b [3]. Hence, the NAM panel presented here may possess these Sr genes, in addition to previously unidentified genes. In an effort to explore the relationship among the QTL detected in our study and previously reported stem rust genes and QTL, we BLASTn searched sequences containing the SNP markers associated with QTL in our study against the Infinium iSelect 9K [69], 90K [70] SNP sequences, and the sequences of GBS-SNPs from the consensus map comprising of GBS-SNP and SSR markers developed by Saintenac et al. [30], in an attempt to determine the position of QTL detected in our study relative to those reported in peer-reviewed studies. However, this yielded in zero high-quality matches among the sequences; although alignments with high number of mismatches and alignment lengths of less than 50% were observed. Additionally, the markers associated with QTL were searched against the putative wheat gene models (v2.2) made available on Ensembl; no matches were found. Though the SNP markers flank a given QTL, they (and the sequences housing the SNP) might not necessarily be of high importance as they are in essence flanking a possible genic region on the chromosome. Therefore, further investigation of the flanked chromosomal regions will likely lead to gene discoveries as well as a better understanding of the nature of disease resistance. While we are unable to draw direct comparisons with previously published QTL mapping studies, our methodology of reference based SNPs used in mapping, thus obtained mapping results, as well as comparative mapping of GBS SNPs with the annotated reference wheat genome in near future could be important steps in the direction of identifying regions with potential candidate genes.

Three QTL were detected on 1A in two environments: StP12 and StP13; two QTL on 1B in three environments: Ken13, StP12, and StP13, of which QSr.umn-1B.2 was common between StP13 and Ken13; and three QTL on 1D in two environments: Ken13 and StP12. Previous studies that have reported QTL associated with APR to stem rust include: Rouse et al. [63] in Thatcher wheat, Bhavani et al. [64] in the CIMMYT bi-parental population PBW343/Kingbird, Yu et al. [60] in CIMMYT’s winter wheat breeding germplasm, Pozniak et al. [71] in durum wheat (Triticum durum Desf.) association mapping panel, and Singh et al. [62] in the durum wheat population Sachem/Strongfield. The TTKSK-effective gene Sr1RS^Amigo [72], though located on 1A, is qualitative in nature and not expected to be present in our panel based on the results from TTKSK screening. Similarly, QTL and genetic markers on 1B that are significantly associated with stem rust resistance have been reported in a durum wheat GWAS study [71], and spring wheat RIL populations [64, 73]. The genomic regions spanned by the QTL in these three chromosomes ranged from 34 Mb (QSr.umn-1D.2) to 184 Mb (QSr.umn-1A.1).

Five, seven, and one QTL were detected on chromosomes 2A, 2B, and 2D, respectively. Three of the seven QTL on 2B lied within < 5 Mb physical distance, with one QTL spanning only about 600 Kb. Several QTL conferring APR to Pgt races, including the Ug99 lineage races, on chromosomes 2A, 2B, and 2D have been detected in spring wheat as well as durum wheat mapping populations and genome-wide association study (GWAS) panels [72, 74–76]. The genes Sr32 (introgressed to 2A, 2B, and 2D), Sr9h, Sr28, Sr36, Sr39, Sr40, Sr47 (located on 2B) are effective to races in the Ug99 race group, but are unlikely to exist in our population based on screening for seedling resistance to race TTKSK. Several other genes such as Sr9a, 9b, 9d, 9e (located on 2BL); Sr21 (located on 2AL), and Sr38 (located on 2AS) are ineffective against African races but are resistant to one or more N. American Pgt races [77, 78]. The QTL QSr.umn-2A.5 is located only 600 Kb from the distal end of 2AL, yet no known Sr genes or QTL are known to exist on this end of 2A. Sr36, also located on 2B, is absent in our panel based on marker screening (S1 File).

While 3A harbored four QTL in each of the four environments, 3B had six QTL of which five were common in at least two environments; two QTL were detected on 3D in Ken13 and StP13 environments. The APR gene Sr2, located on 3BS, is segregating in four populations: LMPG-6/Gem, LMPG-6/Kudu, LMPG-6/Paka, and LMPG-6/Romany. The presence of Sr2 in ‘Gem’, ‘Kudu’, ‘Paka’, and ‘Romany’ was confirmed by marker screening (S1 File) as well as expression of the trait pseudo-black chaff (PBC) in StP12 and Ken13 environments on plant internodes (not observed in Romany). PBC is conditioned by the expression of the partially dominant gene Pbc causing dark coloring of the glumes and inter-nodal regions in an adult wheat plants, and is considered to be associated with Sr2 [79, 80]. The physical length of 3B used to map our sequences is 774 Mb and QTL were detected as early as around 69 cM up to 113 cM on distal 3BS. The gene Sr12 [3], as well as several other QTL identified by mapping studies [64, 71, 81] that confer resistance to stem rust of wheat also exist on 3B. Hence, these QTL could potentially resemble Sr12 or Sr28, of which the latter was absent in our marker screening (S1 File) and the former does not have diagnostic markers available. Sr51, a new gene resistant to Ug99, is located on distal end of 3DS and QSr.umn-3D.2 is also located approximately 7 Mb from the distal end of 3DS. However, Sr51 has been introgressed into bread wheat from Aegilops searsii and is not currently used in agriculture [72, 82]. QSr.umn-3D.2 may represent an ortholog of Sr51 in bread wheat, or could be a new gene.

We detected three QTL on 4A in Ken13, StP12, and StP13; three on 4B; and one on 4D. Chromosome 4A contains the genes Sr7a and 7b, and chromosome 4B hosts the gene Sr37 introgressed from Triticum timopheevi. Sr37 is resistant to the Ug99 races, and therefore is unlikely to be present in our population. Both Sr7a and Sr7b are resistant to multiple N. American Pgt races [77, 83]. Sr7a was first identified in several Kenyan wheat lines, [67, 68]; and ‘Ngiri’ could have acquired this gene from the line ‘Manitou’ [3], which is present in its pedigree. Thus, this gene may have been bred into one or more of the Kenyan varieties used to create our NAM population. Further screening of the population, preferably in single rust race nurseries, is needed to confirm this postulation. Likewise, the stem rust APR gene Sr55 (Lr67/Yr46/Pm46) is located on 4DL. The only QTL detected on 4D (QSr.umn-4D) lies on the short arm of 4D (21 cM on linkage group; between 58–61 Mb on 4D which is 121 Mb long) and therefore may not represent Sr55. The use of markers cfd23 and cfd71 to screen for Lr67/Yr46/ Sr55/Pm46 was inconclusive.

Four QTL were detected on group 5 of chromosomes: two on 5A, one on 5B, and three on 5D of which QSr.umn-5D.1 was common in StP12, StP13, and Ken13 environments, and could be involved in broad-spectrum resistance. Several QTL on chromosome 5A that provide resistance against stem rust have been reported in biparental and association mapping studies [60, 64, 71, 73], although no known Sr gene has been mapped to either 5A. Chromosome 5B contains Sr56, another stem rust APR gene, on its long arm [17]. No diagnostic markers are available to test presence/absence of this gene. Chromosome 5D houses the Sr genes Sr30, which is known to exhibit resistance to N. American Pgt races [57, 72]; and Sr53, which is a recent translocation from Aegilops geniculata [84], and is not expected to be present in our NAM panel.

On chromosome 6A, three QTL were detected in three different environments; four on 6B (of which QSr.umn-6B.2 was common between Ken13 and StP12; and QSr.umn-6B.4 was common between Ken13, StP12, and StP13); chromosome 6D was the only chromosome with no QTL detected. Located on 6A are the genes Sr26 introgressed to common wheat [85, 86]; and Sr13 in durum wheat (Triticum durum Desf.) [71, 87]. Sr26 is effective against all races of the Ug99 race group [57]. Sr13 is not prevalent in common wheat [3, 72]. Chromosome 6A also contains the gene Sr52, a TTKSK resistant gene recently introgressed into hexaploid wheat from its diploid relative Dasypyrum villosum [88]. Hence, it is unlikely that any of these three genes is present in our population. The Sr gene Sr11 is located on 6B [3], in addition to several QTL effective to stem rust on 6A and 6B [60, 64, 71, 86].

QSr.umn-7A, the only QTL on 7A, is located at the very beginning (3 cM) of the linkage group 7A as well as on the physical map (ranging 2–6 Mb). This is strikingly similar to previous QTL reported on 7AS by Bajgain et al. [75] and Singh et al. [65], and could represent the same locus. Chromosome 7B also had only one QTL detected (QSr.umn-7B) which was common in all environments except SA12, and may be an important QTL because of its effectiveness in multiple environments. These QTL, located on the short arm, explained 8%—16% of phenotypic variance. The only gene mapped to 7B is Sr17, which is ineffective to races of the Ug99 lineage as well as to the N. American Pgt races TTTTF, QFCSC, and MCCFC used in this study [77]. Previously reported QTL on chromosome 7B that confer resistance to the Ug99 group of races include the QTL linked to the microsatellite markers wPt-0318 [60] and cfa-2040 [71]. The two QTL detected on 7D were mapped close to each other (QSr.umn-7D.1 at 25 cM and QSr.umn-7D.2 at 26 cM; 16 Mb apart on physical distance) albeit in two different environments. Sr25, Sr43, and Sr44 provide resistance to different races within the Ug99 group and are all located on 7DL. Marker screening for Sr25 showed that the gene is absent in our panel, whereas Sr43 and Sr44 are recently introgressed segments from Thinopyrum elongatum and Thinopyrum intermedium, respectively, and are not currently in use in wheat breeding [72].

Our study uses de novo SNP markers, compared to other studies that use microsatellite or chip-based SNP markers. Hence, a direct comparison for names and positions of markers linked to previously reported Sr genes and QTL is rather difficult to conduct. As indicated above, the lack of sequence alignment among our SNP markers and published SNP assays prevented from making position-based definitive comparisons with previously detected QTL. We expect that additional research, either jointly or in individual populations, will help in further understanding of the relationship between the discussed Sr genes and previously reported QTL with the QTL detected in our study. In particular, construction of a consensus map with previous marker types (SSR, DArT) anchored to the physical map is rather essential to establish a clear connection between sequence-based studies and older studies. Additionally, publication of the wheat genome [45, 46] followed by functional annotation and construction of high density genomic maps can be expected to aid in such comparison, and also in identification of regions with candidate-genes.