A high density GBS map of bread wheat and its application for dissecting complex disease resistance traits

Genotyping-by-sequencing (GBS) is a high-throughput genotyping approach that is starting to be used in several crop species, including bread wheat. Anchoring GBS tags on chromosomes is an important step towards utilizing them for wheat genetic improvement. Here we use genetic linkage mapping to construct a consensus map containing 28644 GBS markers. Three RIL populations, PBW343 × Kingbird, PBW343 × Kenya Swara and PBW343 × Muu, which share a common parent, were used to minimize the impact of potential structural genomic variation on consensus-map quality. The consensus map comprised 3757 unique positions, and the average marker distance was 0.88 cM, obtained by calculating the average distance between two adjacent unique positions. Significant variation of segregation distortion was observed across the three populations. The consensus map was validated by comparing positions of known rust resistance genes, and comparing them to wheat reference genome sequences recently published by the International Wheat Genome Sequencing Consortium, Rye and Ae. tauschii genomes. Three well-characterized rust resistance genes (Sr58/Lr46/Yr29, Sr2/Yr30/Lr27, and Sr57/Lr34/Yr18) and 15 published QTLs for wheat rusts were validated with high resolution. Fifty-two per cent of GBS tags on the consensus map were successfully aligned through BLAST to the right chromosomes on the wheat reference genome sequence. The consensus map should provide a useful basis for analyzing genome-wide variation of complex traits. The identified genes can then be explored as genetic markers to be used in genomic applications in wheat breeding.

Advances in next-generation technologies have driven the costs of DNA sequencing down to the point that genotyping based on sequence data is now feasible for high diversity, large genome species. Genotyping methods usually involve restriction enzyme digestion of target genomes to reduce the complexity at a reasonable cost [7,[15][16][17][18][19][20]. Davey et al. [17] grouped approaches that apply genome complexity reduction methods into the following classes: reduced-representation sequencing, restrictionsite-associated DNA sequencing (RAD-seq), low coverage genotyping including multiplexed shotgun genotyping (MSG) and genotyping by sequencing (GBS). All of these methods are more or less similar technically. The Diversity Arrays Technology (DArT), Canberra, Australia (http:// www.diversityarrays.com/), has developed a GBS platform known as DArT-seq, which provides an opportunity to select genome fractions corresponding predominantly to active genes. Restriction enzymes used in this method separate low copy sequences from the repetitive fraction of the genome. These low copy sequences are informative for marker discovery. Representative fragments are then sequenced on Next Generation Sequencing (NGS) platforms [21,22]. Using a combination of restriction enzymes, DArTseq GBS offers affordable genome profiling through generation of high-density SNPs as well as PAV (presence and absence variations) markers [23][24][25][26]. In a standard DArT assay, approximately 200,000 genomic fragments are sequenced 10 times, on average, with approximately 2,000,000 tags per sample. As a significant percentage of samples in each experiment are processed in duplicate (to enable stringent marker selection based on scoring reproducibility), all sequence variants that are not legitimate SNP markers are filtered out. Large numbers of additional metadata produced by the analytical pipeline (DArTsoft; DArT P/L, Australia http://www.diversityarrays.com/ software.html#dartsoft) make further marker selection and sorting easy, and enable users to choose specific groups of markers that are most useful for their applications [24].
The GBS platform has been used for genetic characterization of more than 40,000 wheat germplasm accessions held by CIMMYT as part of its Seeds of Discovery (SeeD) initiative. A genetic map of GBS markers would be an important prerequisite for trait-based genetic analysis of this large diversity panel. The validity and usefulness of a genetic map depends on its suitability for mapping genomic regions correctly and precisely. In wheat, the Ug99 group of stem rust fungus Puccinia graminis Pers. f. sp. tritici Eriks. & E. Henn. is described as a major forthcoming threat to global wheat production [9,27]. Mapping and deploying adult plant resistance (APR) genes in popular high yielding but susceptible varieties is needed by wheat breeders. To date, more than 50 APR genes have been identified through QTL analyses [28] and some of them (such as Sr2, Lr46 and Lr34) are well characterized and widely used in breeding. Genetic maps that can identify APR genes would be useful for marker-assisted selection (MAS), map-based cloning and detailed molecular characterization in wheat breeding. In this study, three recombinant inbred line (RIL) populations, PBW343 × Kingbird, PBW343 × Kenya Swara and PBW343 × Muu (designated as PB-KB, PB-KS, and PB-MU respectively), were used to anchor the GBS tags to wheat chromosomes. Therefore, the objectives of our investigation were to construct a consensus map of GBS markers and validate known APR genes/QTL against stem rust, yellow rust and leaf rust through mapping, and by comparing them with wheat reference genome sequences recently published by the International Wheat Genome Sequencing Consortium, Rye and Ae. tauschii genomes.

Results
Marker distribution in the consensus map and three individual maps Three RIL populations were analyzed using GBS tags. In total, 13123, 18612 and 6936 markers were mapped on PB-KB, PB-KS and PB-MU populations, respectively ( Table 1, Additional files 1, 2, 3 and 4). Using genotype data of the three populations, a consensus genetic map was constructed by assigning 28644 markers to wheat chromosomes with 3757 unique positions ( Table 1). The maximum ratio of unique positions on the consensus map was located on chromosome 4D (24.4%), whereas chromosome 2D harbored the minimum ratio of unique positions (4.3%). Of the 28644 markers of the consensus map, 32.9%, 56.3% and 10.8% were mapped on the A, B and D genomes, respectively (Table 1). On average, the A, B and D genomes covered distances of 1252.6 cM, 1635.2 cM, and 414.7 cM, respectively.
Total genetic length of the consensus map was 3302.5 cM (Table 1), and average marker distance was 0.88 cM, reached by calculating the average distance between two adjacent unique positions. The length of 3016 marker intervals, corresponding to 80.3% of total marker intervals by 3757 unique positions, ranged from 0 to 1 cM ( Figure 1A). In all three populations, 6.4% of markers on the consensus map were found to be polymorphic; 71.4% of markers were polymorphic in one of the three populations, and 22.3% of them were polymorphic in two of the three populations ( Figure 1B, and Additional file 5). Compared to the three individual maps, the number of markers in common between any two individual maps was roughly 43% of the number of markers in the map with the least markers of the two ( Figure 1C). For example, there were 3002 markers in common between PB-MU and PB-KB, which was 3002/6936 = 43.2% of the number of markers in the PB-MUU map. Also, there were 5729 markers in common between PB-KB and PB-KS, which was 5792/ 13123 = 43.6% of the number of markers in the PB-KB map. Similarly, there were 1189 markers in common among the three maps, which was 1189/6936 = 17.1% of the number of markers in the PB-MU map.

Segregation distortion across the three RIL populations
Segregation distortions were estimated in the three RIL populations. For PB-KB, PB-KS, and PB-MU, 4561 (34.8%), 3686 (19.8%), and 2263 (32.6%) markers, respectively showed evidence of segregation distortion at the 0.05 significance level. The most significant (i.e., -logP) segregation distortion regions (SDRs) were observed in PB-KB, LG α : Number of linkage groups in each chromosome. UP β : Number of unique positions.
followed by the PB-MU and PB-KS populations ( Figure 2). The highest significant SDR was in chromosome 6BL detected in PB-KB, where -logP reached 21.6, and selection favored alleles from PBW343. In both PB-KB and PB-MU, SDRs detected in chromosome 3BS are near the Sr2 gene, where selection favored alleles from PBW343. An SDR detected in PB-KS on chromosome 7DS was in the region around the Lr34 gene, designated by Dyck [29], for resistance to leaf rust. In this region, selection favored alleles from Kenya Swara. Among all chromosomes in the three RIL populations, high segregation distortion was observed on chromosome 1B and selection favoring PBW343.

Distribution of the recombinant inbred lines against wheat rusts
Large variations for stem rust were observed in the three RILs populations. Screening for yellow and leaf rusts was also carried out (in population PB-KS) to detect pleiotropic loci linked with yellow and leaf rust resistance. Phenotypic variations for leaf and yellow rusts were observed as well, but smaller than those for stem rust (Table 1 and Figure 3). In the MS2009 (main season 2009) screening, percent stem rust severity in the PB-KB, PB-KS and PB-MU populations ranged from 10% to 70%, 1% to 60%, and 0% to 80%, respectively. During MS2010 (main season 2010), it ranged from 0% to 60%, 1% to 80%, and 0 to 80% in the PB-KB, PB-KS and PB-MU populations, respectively; in MS2011 (main season 2011), percent severity varied from 1% to 100% in population PB-MU; and in OS2010 (off season 2010), it was 5% to 80%, 0% to 100%, and 5.6% to 76.3% in PB-KB, PB-KS and PB-MU populations, respectively. Population PB-KS was screened for yellow rust in T2010 (Toluca, Mexico 2010), with percent severity ranging from 0% to 100%, and for leaf rust in OB2010 (Obregon, Mexico 2010), with percent severity ranging from 1 to 100%.
Comparison of consensus map with wheat genome reference sequence A draft of the wheat genome sequence was published very recently [30]. To verify the consensus map in this study, we BLAST the sequences of 28644 GBS markers (Additional file 14) against the genome sequence of Chinese Spring [30], rye, and the D genome of Ae. tauschii with E-value < 1e-5. In general, sequences of 3619 GBS markers (that is, 34.9% of the total GBS markers) cannot hit to the genome sequence of Chinese Spring, and sequences of 12.6% of GBS markers cannot be mapped to the same chromosome in the genome sequence of Chinese Spring as that on the consensus map. Genome-wide sequences of 52.4% of GBS markers can be mapped to the same chromosome in the Chinese Spring genome sequence as that on the consensus map ( Figure 5). This ratio ranges from 31.4% to 67.5% across the 21 wheat chromosomes ( Figure 5). For chromosome 1B, there was a very clear enrichment for rye hits indicating the 1B/1R translocation ( Figure 6). There were increased hits to Ae. tauschii on the D genome (Additional file 15).

Discussion
GBS is a preferred high-throughput genotyping method involving targeted complexity reduction and multiplex sequencing to produce high-quality polymorphism data at a relatively low cost per sample. Three RIL populations sharing one common parent (PBW343) were genotyped using the GBS approach. A consensus genetic linkage map distributed by 28644 markers was developed with 3757 unique positions (13.1% of the total number of markers) covering a 3302.45 cM genetic distance (Table 1; Additional files 1, 2, 3 4 and 5). Recently, a wheat genetic map of 40,267 SNP markers was reported [10] where data were generated with the help of SNP iSelect array comprising~90,000 SNPs. On this map, 13.7% of SNPs were specified to unique positions. The percentage of unique positions in the two maps was comparable. Fewer numbers of unique markers on chromosomes were partially due to the lack of polymorphic markers evenly distributed on wheat genome, and the lack of recombination events captured by these populations.
Compared with the A and D genomes, in the B genome, the maximum percentage of the total number of markers (56.3%; Table 1), the maximum percentage of the total number of unique positions (56.5%; Table 1), the longest genetic length (1635.2 cM), and the maximum number of detected QTL regions can be observed (10 out of 18; Table 2 and Figure 4). These results indicate that most of the recombination events happened on the B genome, which was in accordance with previously reported genetic maps [8,10,23,28,31] and also with genome size, since the B genome is the largest, followed by genomes A and D. Variation in the D genome of bread wheat is consistently low [9][10][11]. In the present study, 10.8% of markers on the consensus map were located on the D genome (Table 1). Yet in some regions (chromosome 3D in the PB-KB and PB-KS populations, and chromosome 7D in three populations, etc.), a high number of markers with unique positions can be observed (Additional files 1, 2, 3 and 4). Four QTL regions were detected on the D genome (Table 2). This map can be a useful resource for finding more genes located on the D genome to dissect the traits of interest.
In terms of the marker distribution across populations, the highest number of polymorphic markers was available in the PB-KS population, followed by PB-KB and PB-MU (Table 1). When comparing the maps from different populations, the number of markers in common between any two maps was approximately 43% of the number of markers on the smaller map ( Figure 1C). The number of markers on all three maps was 17.1% of the number of markers on the smallest map. The reduced percentage of markers in common to all maps may be due to the known structural diversity among the parents and the varying recombination frequency patterns across the genome for the three crosses. These trends indicate PBW343 may have the largest genetic distance from Kenya Swara, compared with Kingbird and Muu.
Although we had 28644 markers on the consensus map, polymorphism markers are still lacking in some chromosome regions of the respective RIL population. Possible reasons are: (1) vast parts of the chromosomes of the Triticeae are recombination deserts (the so-called genetic centromeres) [30], so most meiotic recombination events occur in genomic regions that correspond to~20% of the chromosome length, while there is little recombination in 1/3 to 2/3 of the chromosome; in the region of recombination deserts, it is difficult to explore polymorphic markers; and (2) the different variations among four parental genomes, and different genetic distances between PBW343 and the other three parents (Table 1 and Figure 2). In our study, the numbers of linkage groups for the consensus map and each of the three mapping populations were 29, 28, 43, and 35, respectively (Table 1). If two markers are physically located on the same chromosome but very far away from each other, it is very likely that they will act as unlinked loci in a population of limited size. Due to the recombination desert between them, they will fall into different linkage groups.
A popular variety in South Asia, PBW343 is known for harboring the 1B/1R translocation. On the consensus map, there are 3640 markers located on chromosome 1B, with around 2251 markers on its short arm (1BS) and 1389 markers on its long arm (1BL). Only 12.0% of 3640 markers on chromosome 1B were located uniquely, which is a low proportion across 21 chromosomes. Looking at chromosome 1B specifically, 9.6% of markers are located uniquely on chromosome 1BS, while 16.3% of markers are located uniquely on chromosome 1BL. In other words, the vast majority of markers on chromosome 1BS are represented in the map as clusters of co-segregating markers. Interestingly, high segregation distortion was observed in the three populations on chromosome 1B as compared to the others ( Figure 2). Also, regions on chromosome 1B having low hits to the Chinese Spring genome clearly had high hits to the rye genome ( Figure 6). These phenomena could be assigned to the 1B/1R rye translocation in PBW343. In addition to the 1B/1R translocation, translocations 2D/2R and 7B/4R have also been reported in wheat. Rye has been used extensively in CIMMYT wheat breeding, and the three parents (PBW343, MUU and Kingbird) used in this study were derived from CIM-MYT germplasm. In our study, we could not find evidence of these parents carrying the 2D/2R and 7B/ 4R translocations.
The consensus map constructed in this study can be used to locate major genes controlling target traits. Phenotypic variation in the three RIL populations suggests polygenic inheritance for APR to stem rust race Ug99 (Table 1 and Figure 3). QTL analysis revealed a couple of APR QTLs against stem rust fungus Ug99, which were mapped as reported. To anchor and determine the relationship between the APR QTLs found in the present study and those found in previous reports (Table 1), we calculated the correlation coefficient of the genotypes of array-based DArT markers and DArT-seq markers on the consensus map. Array-based DArT markers wPt-744022 and wPt-5896 are reported to be linked with APR to stem rust [9]. QTL QSr.cim-2BS1 reported in this study was flanked by DArT-seq markers, i.e., 4989818 and 1088282. The correlation between 4989818 and wPt-744022 was 0.81 in population PB-MU, which is highly significant. DArT-seq markers 1298718 and 1025982 were the two markers flanking QSr.cim-5BL ( Table 2). The correlation between wPt-5896 and DArT-seq marker 1298718 was 0.91 in population PB-MU, which is highly significant as well (Additional files 1, 2, 3, 4 and 5). A QTL on chromosome 3BS, which is most likely the Sr2 gene (in Table 2 designated as Sr2) was flanked by DArT-seq markers 1106039 and 1140316 in population PB-MUU. According to the consensus map reported by Yu et al. [31], the Sr2 gene is 9.2 cM apart from the array-based DArT marker wPt-3761. Marker 1140316 was highly correlated with wPt-3761 in population PB-MUU (correlation coefficient: 0.75). A gene for leaf rust resistance on chromosome 7DS was mapped and flanked by DArT-seq markers 1128052 and 4991056 ( Table 2). The correlation between 4991056 and wPt-733087 was 0.69, which is highly significant. Array-based DArT marker wPt-733087 was associated with leaf rust resistance and was found to be 9 cM apart from Lr34 gene-based marker csLV34 in the PBW343 × Diniza population (Singh et al., CIMMYT, unpublished). In addition to marker co-location for the above mentioned QTLs/gene(s), a recently published consensus map for Ug99 stem rust resistance loci in wheat [31] was compared with QTLs mapped in this study. Since this map did not contain GBS markers, exact co-location could not be made. However, based on the location of APR QTLs on chromosome arms for some of the QTLs, i.e., QSr.cim-2BS1, QSr.cim-2BS2, QSr.cim-2BL, QSr.cim-2DS, QSr.cim-3AL, QSr.cim-3DS, QSr.cim-1BL1 and QSr.cim-6DL could be co-located with the ones projected on the consensus map of Ug99 stem rust resistance [31]. Singh et al. [9] also reported an APR QTL for Ug99 on chromosome 7AS, which is likely the same as the one on that chromosome in the present study (Table 2 and Figure 4). QSr.cim-3BS1, QSr.cim-3BS2, QSr.cim-4AS, and QSr/Yr.cim-6BL were new QTL regions associated with stem rust and yellow rust. The high density GBS consensus map increased the mapping resolution of linkage mapping. Identified genomic regions (i.e., the genetic region of each QTL's flanking-marker interval in each individual linkage map) for stem rust resistance ranged from 0.1 to 15.8 cM (Additional file 6). The interval size of 14 of the QTLs was < =1 cM and most of them were < 5 cM (Additional file 6). QTLs in genomic regions of this size are valuable for further understanding the molecular basis and developing perfectly linked markers. Co-location of the APR QTLs/genes in their respective chromosomal regions (Table 2 and Figure 4) and a high proportion of markers BLAST to the correct chromosomes of the genome sequence of Chinese Spring ( Figure 5) indicate the validity and utility of the consensus map.
CIMMYT's Seeds of Discovery project has characterized more than 40,000 wheat gene bank accessions through the DArTseq GBS platform. The high density GBS consensus map reported in this study is an essential prerequisite for analyzing the GBS data of a large diversity panel. It will facilitate the genetic dissection of important quantitative traits either by linkage mapping as we reported in this paper, or by genome-wide association mapping. GBS markers associated with important traits can be utilized by designing primers according to their sequence, for genomics applications in wheat breeding.

Conclusions
A high density map of 28644 GBS markers using genotypic data of the three RIL populations with a common parent, PBW343. Total genetic length of map was 3302.5 cM with 3757 unique positions, and the average marker distance was 0.88 cM by calculating the averaged distance between two adjacent unique positions. The length of marker intervals ranged from 0 to 28.3 cM. The number of markers in common between any two individual maps was roughly 43% of the number of markers in the map with the least markers of the two. Significant variation of segregation distortion was observed across three populations. Three genes (Sr58/Lr46/Yr29, Sr2/ Yr30/Lr27, and Sr57/Lr34/Yr18) and 15 published QTL were validated. The common parent PBW343 harbors the 1B/1R translocation, and there was a very clear enrichment for rye hits on chromosome 1B. Also, there were increased hits to Ae. Tauschii D genome. The high-density and better quality of genetic maps will advance the genetic studies of complex trait in wheat and facilitate genomics-assisted breeding.

Plant materials
Three RIL populations were used for consensus map construction and QTL analysis for rust resistance. Moderately susceptible bread wheat (Triticum aestivum) parent PBW343 was used as a common parent and crossed with three other bread wheat APR parents: Kingbird, Kenya Swara, and Muu. PBW343, a major variety in India, is a selection (GID2430154) from CIMMYT line Attila with the pedigree Nord Deprez/VG9144//Kalyansona/Bluebird/3/Yaco/4/ Veery#5 [9]. Parents Kingbird and Kenya Swara have maintained high levels of APR to stem rust. Kingbird has shown a high level of APR in field tests conducted at Njoro, Kenya, during the last six cropping seasons, including the 2008 main season, which was characterized by very high stem rust pressure. Muu (pedigree: Pfau/Weaver//Kiritati; GID5090613) was found to be susceptible to wheat stem rust at the seedling stage but adult plants showed low disease severity in response to stem rust race Ug99 during multiple years of field trials in Kenya [32].

DArTseq™ GBS markers
DArTseq is a GBS platform developed by DArT PL, Canberra, Australia. It is a combination of complexity reduction methods developed initially for array-based DArT and sequencing of resulting representations on next-generation sequencing platforms. For sequencing- based DArT genotyping, two complexity reduction methods optimized for several other plant species at DArT PL (i.e., PstI/HpaII and PstI/HhaI) were used to select a subset of PstI-HpaII and PstI-HhaI fragments, respectively [23]. DNA samples were genotyped twice using two different 4-bp cutters on one end of the RE fragments (HpaII and HhaI). Although the general concept of discovering markers through sequencing of genomic representations was presented over two decades ago [33,34], it was only very recently that the cost and throughput of sequencers reached a point where any GBStype approach can compete effectively with hybridizationbased arrays (DArT and/or SNP; [11]).
The DArTsoft marker extraction pipeline produced large numbers of markers in each of the three populations. Markers were filtered on the basis of reproducibility (that is, the percentage of technical replicate pairs scoring identically for a given marker), call rate (that is, the percentage of samples for which a given marker was scored), and the average read depth (that is, the average number of sequence 'tag' counts contributing to the genotype calls for a given marker). Approximately 60% of samples from each population were assayed twice to derive reproducibility scores. The minimum threshold value for reproducibility was 95%, with an average value of 98.5% for SNPs and 99% for silicoDArTs. The minimum threshold value for call rate was 85%, with an average value of 99% for SNPs and 95% for silicoDArTs. The minimum threshold value for Average Read Depth for SNPs was 7, with an average of 18; for silicoDArTs, it was 8, with an average of 17.2. Markers with identical genotypes were placed in redundant bins, and markers with unique genotypes (those that did not belong in a redundant bin) were excluded from the mapping process. This provided an additional quality control step for marker selection. The markers were selected to minimize the number of missing calls as per the selection criterion. The only filtering for 'false homozygote' calls was the masking of apparent double crossovers after ordering of the linkage groups. Genotyping errors will present themselves as apparent double crossovers in the ordered map data as 'singletons' (data points with genotype scores that differ from those of the immediately preceding and following markers) [35]. For SNP calling, the variants were called within the data only (clustering sequences by sequence similarity), and no external reference genome was used. The sequence defined as the 'Reference' for each SNP pair was either the most common sequence in the population or the sequence that had been previously recorded from DArT genotyping analyses of wheat.

High-density linkage map and consensus map construction
Our three individual maps were constructed using DArT PL's OCD MAPPING program [36], which implements a marker-ordering algorithm combined with a tunable double cross-over (DCO) masking algorithm. Markers were clustered into linkage groups according to the method described by Anderson et al. [37]. Markers with identical genotypes were placed in redundant bins, and the resulting markers/bins within each linkage group were ordered using the traveling salesman path solver program Concorde [38]. Since silicoDArT markers were dominant markers, while SNP markers were co-dominant markers, silicoDArT markers were separated into paternal and maternal phases. A map was produced for each parent by combining the relevant silicoDArT markers with all of the SNP markers. This resulted in two maps that were joined in a single population consensus map using the SNP markers in common to facilitate consensus map construction. Apparent double-crossovers were masked before reordering the linkage groups and calculating recombination fractions, with Kosambi function used to estimate genetic distances.
Construction of consensus maps presents a challenging problem in wheat due to its structural diversity, particularly the chromosomal structure differences between the parents of mapping populations. We applied DArT PL's OCD MAPPING program [36] to order DArTseq, arraybased DArTs, and SSR markers. We then applied DArT PL's consensus mapping software [24,36].

Evaluating stem, leaf and yellow rust severity
The parents, highly susceptible bread wheat check variety, Cacuke and the three RIL populations were evaluated for stem rust severity at the Kenya Agricultural Research Institute (KARI) in Njoro during four crop seasons: main season 2009, main and off-season 2010, and main season 2011, denoted as Sr-MS2009, Sr-MS2010, Sr-OS2010, and Sr-MS2011, respectively. All three populations were evaluated for APR to stem rust, and PBW343 x Kenya Swara was screened for yellow and leaf rusts in one season at one location (Additional file 12). The RILs and parents were sown using a completely randomized design with two replicates. Field plots consisted of two 1-m rows spaced 20 cm apart with a 0.5-m pathway. Approximately 60-70 seeds were sown in each plot. The experimental block was surrounded by a spreader row consisting of varieties differentially susceptible to the Sr24 virulent variant TTKST. Hill plots of spreaders were also planted in the middle of the pathway on one side of each plot to facilitate uniform disease build-up and spread. On at least two occasions just prior to booting, freshly collected urediniospores suspended in distilled water were injected into culms in the spreader plots (1-3 plants/m) using a hypodermic syringe. Disease response in the field was assessed twice, first, when the susceptible check variety Cacuke displayed 50-60% stem rust severity and subsequently at peak disease development, when Cacuke displayed 100% stem rust at the mid-dough stage of plant growth. Percent disease severity was scored using the modified Cobb Scale [39]. The second rating was considered as the phenotype in this study.
Parents and population lines were evaluated for leaf rust reaction in field nurseries operated by CIMMYT in Ciudad Obregon, Sonora, Mexico, in 2010, denoted as Lr-OB2010. Replicated trials with parents and population lines were grown in Obregon in 2010. Leaf rust severity in each plot was visually scored (near anthesis flowering time) using the modified Cobb scale [39]. For yellow rust screening, parents and population lines were evaluated under field conditions in rust nurseries operated by CIMMYT near Toluca, Edo. Mexico, Mexico, in 2010 and in Njoro, Kenya, in 2010 and 2011, which were denoted as Yr-T2010, Yr-K2010, and Yr-K2011, respectively. Two replicated rows of parents and population lines were assessed in each trial. Each row was visually scored around anthesis flowering time for yellow rust severity with the percentage of leaf covered with disease infection calculated as described for leaf rust. Considering the low phenotypic variance of yellow and leaf rust resistance in the PB-KB and PB-MU populations (data not shown), we did not use them to do QTL analysis in this study. Phenotypic distribution and correlation of rust resistance across the three RIL populations are shown in Figure 3 and Additional file 12.

QTL mapping of the individual RIL populations
Inclusive composite interval mapping (ICIM) [40][41][42] implemented by QTL IciMapping 3.2 (available at www. isbreeding.net) was used to map additive and epistatic QTLs controlling stem, leaf and yellow rust resistance. In ICIM for additive QTL mapping, marker selection was performed just once using stepwise regression and considering all marker information simultaneously; this is a key step in determining the scanning profile. Phenotypic values were then adjusted by all markers retained in the regression equation, except the two markers flanking the current mapping interval [40][41][42]. Permutation tests were conducted using stem rust in the three RIL populations to determine the criteria for model selection in the first step of ICIM. For all three RILs, the probability of a marker moving into the model corresponding to the overall type I error α = 0.05 was approximately 10 -5 . The probability of a marker moving out of the model was set at twice the probability of a marker moving into the model. The LOD threshold to declare the existence of a QTL was calculated by permutation tests as well. Permutation tests revealed LOD thresholds of 3.50, 3.53, and 3.51 for PB-KB, PB-KS, and PB-MU, respectively. Considering that thresholds retained from permutation tests are always conservative [43], an LOD threshold of 2.5 was used to report QTLs and determine common QTLs across seasons and populations.
For epistatic QTL mapping, we tested all possible pairs of scanning positions by ICIM [41]. That is to say, we can detect digenic interactions regardless of whether the two interacting QTLs have significant additive effects or not. Due to the large amount of variables in digenic QTL mapping, we used a much stricter probability (10 -6 ) of a marker moving into the model. The probability of a marker moving out of the model was set at twice the probability of a marker moving into the model. An empirical LOD threshold of 4.0 was used to declare the existence of epistatic QTLs.

Common QTL across the three RIL populations
Due to the differences in the three individual linkage maps, it was difficult to directly find common QTLs across the three RIL populations based on the QTL or marker position in each linkage map. Therefore, we projected each QTL's flanking markers to the consensus map. If the flanking markers of one QTL were 15 cM apart from the flanking markers of another QTL on both sides, the two QTLs were declared as common QTLs.