High-density SNP-based genetic map development and linkage disequilibrium assessment in Brassica napusL

Genetic linkage maps are highly valuable tools for comparative genome analyses and the identification of genomic regions carrying major genes and quantitative trait loci (QTL) controlling agronomical traits. They are a prerequisite for further map-based cloning or marker-assisted breeding programs. In recent years, the establishment of genetic maps have benefited from the development of new types of molecular markers which take advantage of automated sequencing and genotyping technologies. While the first marker-based genetic maps were built with restriction fragment length polymorphisms (RFLPs), random amplified polymorphic DNAs (RAPDs) and amplified fragment length polymorphisms (AFLPs), dense genetic maps now include simple sequence repeats (SSRs) and more recently single nucleotide polymorphisms (SNPs). Dense genetic maps based on sequence-derived markers allow finer comparative genome analyses to be performed based on comparisons with sequenced related genomes and to accelerate the process of map-based cloning of major genes and QTL. They are also very useful tools to assist sequence assembly in whole de novo genome sequencing projects [1–3]. Moreover, by integrating genetic map data with genotyping data generated from collections of accessions/varieties linkage disequilibrium (LD) pattern along the genome of a given species can be investigated, which is a prerequisite for precise genome wide association studies (GWAS). GWAS performed with a large number of SNPs have been reported in a number of crop species such as maize [4, 5], Arabidopsis [6], barley [7–11], and rice [12–15] The success of GWAS to locate genes responsible for complex traits depends on the extent of LD, the number, the distribution and the diversity of markers and the underlying structure in the studied collections. Since the diversity of markers and the extent of LD may vary depending on the history of the collections [7, 15], they should be investigated prior to GWAS design.

Oilseed rape (Brassica napus) is a prominent oilseed crop in most world continents including America, Europe, Australia and Asia and is cultivated for food (oil) and feed (meal) as well as for non-food uses such as biofuels or lubricants. It is the second world oil crop after soybean (http://faostat3.fao.org/home/index.html; May 2011) with a world production of more than 60 million tonnes per year. B. napus is an amphidiploid species (AC genome, n = 19) that arose from hybridization between B. rapa (A genome, n = 10) and B. oleracea (C genome, n = 9) diploid species [16] within the past 10,000 years [17]. B. napus includes spring and winter oilseed rape, rutabaga or swede, and some fodder crops. It likely originated from a few interspecific hybridization events [18] and has only a short domestication history of about 400–500 years [17, 19]. For these reasons, the genetic diversity within B. napus germplasm is rather low compared to that of its two progenitor species B. rapa and B. oleracea. Moreover, two bottlenecks have occurred during breeding of modern oilseed rape varieties through the selection for low erucic acid content in the oil and low glucosinolate content in the seeds, which reduced the genetic diversity in modern varieties [20].

Over the last 20 years, many genetic B. napus maps have been built, which have been progressively integrating various types of markers [21–27]. These maps have been used for genetic studies of various agronomical traits including development traits [28, 29], seed quality [30–34], yield components [35–38] and disease resistance [39–44] as well as for genetic study of chromosome pairing [45]. The establishment of genetic maps of diploid and amphidiploid Brassica species, and their comparison and alignment to Arabidopsis genome sequence provided insights into Brassica genome organization and evolution after the different rounds of polyploidization and diploidization occurring in these species history. Extensive collinearity observed between A. thaliana and B. napus led to the description of a genomic block system determined by Parkin et al. [21], who demonstrated that the structure of the Brassica A and C genomes could be described with approximately 21 conserved blocks. A framework built of 24 genomic blocks (A-X) within the ancestral karyotype was then proposed that represents an extension of the above mentioned study [46]. This conserved block structure was then further investigated in related species such as B. juncea[47] or B. oleracea[48] and led to a block arrangement comparison in the A, B and C genomes [49]. It was also recently confirmed in B. napus using dense genetic maps with SSR [27], and SNP [50] markers.

The availability of high numbers of markers now makes it possible to investigate more precisely genome wide diversity and the extent of LD in oilseed rape (OSR). To date, the published studies relied on either a low number of lines or low number of markers. Ecke et al. [51] used 845 AFLP markers to examine the extent of LD in a set of 85 winter oilseed rape lines. Bus et al. [20] investigated patterns of genetic diversity and the extent of LD in 509 inbred lines corresponding to different germplasms of oilseed rape with 89 SSR markers. Xiao et al. [52] assessed the genetic diversity and the extent of LD in a panel of 192 inbred lines from all over the world but with a great proportion originating from China using 451 SSR markers. Harper et al. [53] carried out associative transcriptomics on 53 B. napus lines using >60 K SNPs and confirmed the low overall level of LD in B. napus.

In this context, our objectives were to (i) obtain a dense integrated SNP genetic map of B. napus built from four segregating populations (ii) investigate the polymorphism of these SNPs within and among different germplasm types of a large B. napus collection and (iii) assess the extent and pattern of LD between densely mapped SNP markers.

In this study, we could built a high density SNP integrated B. napus map and depict polymorphism level and LD decay over the linkage groups across different B. napus collections by integrating genotyping data of a large set of SNPs in both segregating populations and diverse collections.

As reviewed by Kaur et al. [54], SNP discovery is challenging in allopolyploid species such as B. napus. SNPs may arise both between allelic (homologous) sequences within subgenomes and between homoeologous sequences among subgenomes but also from polymorphisms between paralogous duplicated sequences. SNP discovery has been based on B. napus ESTs sequence analysis [55–57] or on second generation high throughput sequencing [50, 58–60]. SNP genotyping using Illumina GoldenGate assays was shown to be possible in B. napus species after careful selection of the SNP [55]. Here, we report the first high throughput genotyping study in oilseed rape using Illumina Infinium and GoldenGate technology. From the 1536 SNPs tested on the two platforms, 1104 were validated on both, suggesting that it should be possible to use one platform or the other depending on the required number of SNPs.

In this study, we generated an integrated map with 7367 markers including 5764 SNPs, which corresponds to 3.3 markers every cM. This large number of mapped markers was obtained by integrating four mapping populations. Marker density actually increased from one marker (one SNP) every 0.54 (0.7-0.8) cM on the individual maps to one marker (SNP) every 0.3 (0.39) cM on the integrated map. This density was comparable to that obtained in Wang et al. [27] (one marker every 0.34 cM) but the number of markers in common between populations was much higher in our study (60% of the markers were common to at least two maps in our study compared to only 20% in Wang et al. [27]). The size of the population was also larger in our study, which led to a more accurate ordering of the markers on the integrated map. Bancroft et al. [50] built a B. napus map with 23027 SNPs by transcriptome sequencing of 37 TNDH lines. These SNP were distributed in 527 recombination bins (one bin corresponded to SNPs having the same scoring genotyping data). Our integrated maps exhibited 3177 bin loci (one bin corresponded to SNPs mapped within 0.1 cM), thus the mapping resolution increased considerably with the increase in population size. The collinearity blocks we identified were compared to those reported in Panjabi et al. [47], Parkin [49], Wang et al. [27] and Bancroft et al. [50]. Compared to Bancroft et al. [50], which is at present the most complete study, five and 10 small blocks were missing on the A and C genomes, respectively. Nevertheless, three new blocks were identified in our study: S on A8, V-W on C2 and B on C7. In Bancroft’s study, two markers corresponding to these blocks were identified but not declared as a block. Moreover, the V-W and B blocks were identified in the homoeologous regions on A2 and A7 linkage groups, which supports their occurrence in these genomic regions.

A higher number of markers was mapped on the A genome than on the C genome, as previously reported by Bancroft et al. [50]. As in their study, here, this difference between the A and C genomes was more pronounced on the three crosses involving Asiatic parental lines. The hypothesis is that Asiatic cultivars are partly derived from crosses involving its progenitor species B. rapa[61, 62], which increased the genetic diversity of the A genome. This introgression of B. rapa genetic information was recently shown for ‘Nignyou7’ (one of the parent of the TNDH map) by Bancroft et al. [50]. Indeed, there was less difference between the number of markers mapped on the A and C genomes in the AMDH population which is derived from a cross between two French winter OSR varieties. No such difference was observed in Wang et al. [27] but the markers mapped in their study were mainly contributed by two crosses involving resynthetised B. napus, which is enriched in polymorphisms in both genomes. We probably succeeded in capturing this high level of polymorphism on the A genome because the original sequencing to identify the SNP was performed on material that included Asiatic varieties.

We observed a high level of segregation distortion especially in the DYDH and TNDY populations. Such segregation distortions were reported in many B. napus maps (e.g. [27]). Segregation distortion and clustering of the skewed loci are common features of microspore-derived DH populations in various species ([63] for a review), including oilseed rape and may be related to differential responsiveness to microspore culture between the two parental lines, which leads to skewed loci in regions involved in the microspore culture responsiveness.

Very good collinearity was observed between all the individual maps, which made it easy to integrate the four individual maps accurately. Only three inversions were identified on the A2 and C8 linkage groups. These inversions could be due to mapping inaccuracies and need to be confirmed. Due to careful selection at the beginning of the study to target homologous SNPs, very few SNPs were not assigned a position on the same linkage groups among the different maps. Those SNPs that did map to different linkage groups were located in duplicated regions within or between the A and C genomes where there is a high level sequence similarity. The map derived from the AADH population was 75% bigger than those derived from the other crosses, which was expected from the way the DH population was obtained. DH lines were produced after intermating F2 plants, which increased the number of recombination events. This type of highly recombinant population is of interest for obtaining better mapping resolution [64, 65]. The map derived from the AMDH population was the least dense with some regions missing, due to the lack of polymorphism between the two parental lines. In many cases, the monomorphic regions corresponded to quite complete B. napus/Arabidopsis collinearity blocks. The lower marker density can be related to the lower level of polymorphism revealed within winter oilseed rape or “00” winter oilseed rape compared with spring oilseed rape.

A moderate level of differentiation was observed between the different B. napus types (WOSR, SOSR, fodder rape) as revealed by the estimated differentiation indices (0.18 between the three groups and 0.22 between SOSR and WOSR). A similar result was reported by Bus et al. [20]. They used 89 SSR primer combinations to assess the diversity in a set of 509 oilseed rape lines which included WOSR, SOSR, fodder and swede rape lines from diverse origins. Xiao et al. [52] found a lower level of differentiation between their groups defined from a collection of 192 oilseed rape lines genotyped with 451 SSR markers. However, each of their groups, as defined by STRUCTURE analysis, was constituted of lines from China, Europe, Canada and Australia. Within the SOSR and WOSR groups, a very low level of differentiation was observed between the three subgroups which corresponded to the three quality types (“++”, “0+” and “00”) with 91% and 83% of the variation present within the subgroups in WOSR and SOSR, respectively. Examination of the collection structure with PCA showed distinct clustering of WOSR, SOSR and fodder rape lines while the first two axes did not account for a large part of the variation (19.5%). This differentiation was previously reported by Diers and Osborn [66], Hasan et al. [67] and Bus et al. [20] and can be related to the relatively distinct breeding history between these pools and their adaptation to different environments or uses. However, our data allowed two groups to be differentiated within SOSR which mainly contained European and Canadian or Asian spring oilseed rape, thus corresponding to a differentiation due to geographic origin. Some SOSR lines were located at an intermediate position between these two groups such as ‘Grouse’ and ‘Marnoo’ (Australian cvs), ‘Chine Wuhan’ and ‘Yeong Dang’ (Asiatic cvs) or ‘Industry’ (a European cv with high erucic acid and low glucosinolate content). No such differentiation was observed within the WOSR lines which all originated from Europe, as previously reported by Ecke et al. [51] and Bus et al. [20]. Fodder rape lines were located at an intermediate position between SOSR and WOSR although closer to WOSR than to SOSR. The exception was ‘Liho’ which is a spring fodder rape and grouped with the SOSR. This result is consistent with the fact that fodder types and oilseed types were derived from the same ancestral spring and winter pools and were bred for different uses, as reported by Bus et al. [20]. The same structuration was obtained whether we consider SNPs from the A or the C genomes. The percentage of variation accounted for by the first two axes was higher with the SNPs from the A than from the C genome when we used the whole SNP set where the markers from the A genome were overrepresented. When we considered the 1507 chosen SNP markers, the percentage of variation accounted for by the first two axes was similar with the SNPs from both genomes, indicating that an overrepresentation of markers in some regions can biased the results. A selection of evenly spaced markers over all the linkage groups is thus recommended for genetic structuration. The optimal number has to be assessed by testing different sample sizes.

To assess the genetic diversity in the collection and within the different germplasm groups, we examined the mean PIC and its evolution along each linkage group. The PIC values ranged between 0.1 and 0.35 depending on the position on the linkage groups and depending on the collections used. Similar observations were made in barley and maize with the same type of markers [4, 10, 68]. Our set of SNPs was only derived from exonic sequences, which could have lowered the level of revealed diversity compared to intronic SNPs [4]. On average, PIC values were lower in WOSR than in SOSR. This difference was more important for SNPs derived from the A genome rather than the C genome. This might be because the SOSR lines have more diverse geographic origins than the WOSR lines and have undergone differential selection to adapt to the different continents. In WOSR, the mean PIC values were not very different between the three seed quality types. The mean PIC value was only slightly higher in “++” WOSR type but this result should be taken with caution due to the different sizes of the three subgroups. Bus et al. [20] reported a lower genetic diversity for “00” seed quality WOSR varieties and the difference was also quite low. The variation in PIC values along the different linkage groups between the three WOSR seed quality types was very contrasted, with more numerous regions showing a PIC decrease for the “0+” and/or “00” types. These variations could be related to potential selection signatures within and between these types. In barley, specific chromosomal regions exhibited contrasting levels of diversity in different germplasm subgroups. A region of reduced diversity in winter barley in the central part of chromosome 5H was attributed to the small number of founding genotypes that contributed to the winter seasonal growth habit locus Hrn-1. Similarly, an abrupt decrease in diversity on the short arm of chromosome 3H observed in all groups coincided with the locus of non-shattering of ears after ripening [68]. We therefore investigated whether the regions with reduced diversity in either the “0+” and /or the “00” subgroups could be related to the position of genes controlling erucic acid and glucosinolate content. Erucic acid genes are located on A8 and C3. Numerous genes are involved in the glucosinolate pathway [69, 70] but the major QTL controlling total glucosinolate content are located on A2, A9, C2, C7 and C9 [71, 72], Delourme, unpubl. data. Their positions are indicated in Additional file 2: Table 2. However, there was no decrease in diversity clearly surrounding erucic acid genes in the “0+” and the “00” subgroups or total glucosinolate QTL in the “00” subgroup. The only exception could be at the top of C9. It can be hypothesized that many recombination events have occurred around the selected genes in these regions during different rounds of intercrossing between the varieties since the original crosses were made with the genitors of low erucic and low glucosinolate content. Significant decreases in diversity were observed in other regions, which could be related to breeding for other agronomical traits but a more precise investigation of QTL located in these regions should be made before drawing any conclusions.

The mean pairwise r² values are close to previous estimates of 0.027 [51] or 0.0247 [53] and confirms the low overall level of linkage disequilibrium in B. napus. LD was observed to decay below a critical level (r² value 0.2) within a map distance between 0.6 and 1.2 cM among the subgroups. This value is in accordance with previous studies performed either on smaller oilseed rape collections or with a smaller number of markers [20, 52]. This level is lower than that detected in a collection of 85 WOSR lines genotyped with 845 AFLP markers, where the LD decayed within 2–3 cM at r² < 0.2 [51]. This could be due to the fact that this latter collection comprised only “00” seed quality types. In our study, the extent of LD was also higher in the “00” WOSR collection than in the whole WOSR or the SOSR collections. From the size of the genetic maps and of the genomes in Brassica species, on average 1 cM can be roughly estimated to correspond to ~500 kb [51, 73]. This means that on average for the whole genome, the extent of LD is between 300–1000 kb depending on the collections but great variations were observed across the linkage groups (from ~100 kb on A10 to ~1700 kb on C9 in SOSR). In many species, LD decay varied across the germplasms used. For example, LD extends less than 1 kb for maize landraces [74] and roughly 2 kb for diverse inbred maize lines [75] but can be as high as 100 kb for commercial elite inbred lines [76]. LD decay can also vary considerably from locus to locus 1–4 kb [77] up to 800 kb [78]. Similar differences were observed in rice, 50–500 kb, [12, 15, 79], Arabidopsis, 10–250 kb [80, 81] or barley, 90–210 kb [82] or 4–8 cM [83]. Generally, the extent of LD is related to the mating system of the species, the breeding history of the species (e.g. the occurrence of bottlenecks) and the genetic diversity of the different germplasms [84]. LD decay is more rapid in outcrossing species and in pools with higher genetic diversity. Oilseed rape is bred as a selfing species and has undergone two bottlenecks for seed quality improvement (to eliminate erucic acid and decrease glucosinolate content), which led to a LD decay similar to other selfing species and to a higher LD extent in the “00” WOSR collection.

The LD patterns varied greatly among the linkage groups and these variations were different between the SOSR or WOSR types. Similar results were reported in maize [4, 5], barley [85] or rice [13]. Different patterns of LD along the chromosomes in various pools can be related to variation in recombination rate and in the history of recombination for specific chromosome regions within these pools. The centromeres were localised approximately on the integrated map by mapping centromeric markers developed by Pouilly et al. [86] to the DYDH map. These were then included on the integrated map and on the LD heatmaps (Additional file 2: Table S2 and Additional file 9: Figure S5). On many linkage groups, LD seems to be extended across the centromeric regions as reported in barley [10] but other regions which do not correspond to centromeres also showed extended LD such as on A2, A8, C8 or C9. Differences in allele frequencies [87] could have also influenced the distribution of LD along the chromosomes. In barley, larger LD extent in some chromosomic regions was caused by markers with low allele frequencies [85]. In rice, LD decay rates in indica and japonica subspecies were only weakly correlated across the genome in relation to a relatively long history of partial reproductive isolation of these self-fertilized subspecies [13]. Since SOSR and WOSR have a relatively distinct breeding history, a similar hypothesis can also be proposed in our case.

With high throughput SNP genotyping on four segregating DH populations, we developed an integrated genetic map for oilseed rape that comprises 5764 SNP and 1603 PCR markers. This significantly improves the marker density or mapping accuracy compared to previously published genetic maps. The genotyping of these mapped SNP markers in collections allowed polymorphism level and linkage disequilibrium to be studied in oilseed rape. Both were shown to vary across the different collections (winter vs spring, seed quality types) and across the linkage groups. Taking into account the length of the genetic map (~2500 cM) and the mean LD extent (0.7 – 1.2 cM for r² > 0.2), a relatively low number of evenly spaced SNPs (few thousands) would be necessary to perform genome wide association studies in oilseed rape. However, this number should be adjusted to obtain a sufficient SNP density throughout the genome and to take into account the variation in LD along the linkage groups. A set of 1536 public SNPs was set up, of which 72% were validated on a GoldenGate platform. They provide evenly spaced SNPs showing a good level of polymorphism in oilseed rape. Information regarding the other SNPs can be requested from J Pauquet (Jerome.Pauquet@biogemma.com). Our study provides a valuable resource for further genetic studies through linkage or association mapping, marker assisted breeding and Brassica sequence assembly and comparative mapping.