DH mapping population (185 lines) of B. carinata was developed from an F₁ cross between the DH parental lines, Y-BcDH64 and W-BcDH76 (Guo et al., 2012). In addition, a set of 81 diverse accessions (Supplementary Table 1; Jiang et al., 2007) of B. carinata including the parents of the YW DH population, mainly of Ethiopian origin, collected from Centre for Genetic Resources, Wageningen, The Netherlands and Germany, were used for genetic diversity analyses. All accessions were selfed to avoid any cross pollination.

DNA from the B. carinata accessions was extracted from bulked young leaf tissue using the DNA extraction Kit (DNeasy Plant Mini Kit, QIAGEN) and sent to DArT P/L (www.diversityarrays.com) for DArTseq based genotyping as described previously (Raman et al., 2014). Both presence-absence markers and SNP markers, collectively called “DArTseq markers,” were scored by DArT P/L.

The 69-bp long sequence of the filtered DArTseq markers identified in 81 B. carinata accessions were aligned to the reference genome of B. nigra (Yang et al., 2016), B. napus (Chalhoub et al., 2014), and B. oleracea (Parkin et al., 2014) by Blast using an e-value threshold of e⁻¹⁰ (Altschul et al., 1990) with ≥60 bp match length. The top blast hits for the sequences were assigned as the physical position of the corresponding markers. The blast matches to multiple loci with the same top e-value were considered with multiple positions, and assigned to the class of unassigned markers without unique positions to the reference genome.

In order to compare the QTL detected in the C genome of YW DH population with the QTL or candidate gene identified in B. napus, the physical position of the markers mapping to the C genome of YW DH population were assigned according to the top blast hits for the sequences of each marker against the reference C genome of B. napus (Chalhoub et al., 2014) by Blast using an e-value threshold of e⁻¹⁰ (Altschul et al., 1990) with ≥62 bp match length. If the blast matches to multiple loci had the same top e-value, the locus which showed the same chromosomal localization of the marker as in the genetic and physical map, was considered as the right location of the marker.

The population structure was estimated with software Structure 2.3.4 using Bayesian clustering and the admixture model (Pritchard et al., 2000). The number of subgroups (K) was set from 1 to 10 with three independent runs. The optimum number of subpopulations was determined by log likelihood of the data [LnP(D)] and ΔK method described by Evanno et al. (2005). The Q matrix was assigned to a sub-population (Remington et al., 2001). Nei's genetic distance (Nei et al., 1983) was calculated and used for the unrooted phylogeny reconstruction using a neighbor-joining (NJ) method as implemented in Powermaker version 3.25 (Liu and Muse, 2005). Phylogenetic tree was viewed using MEGA 4.0 (Tamura et al., 2007). The markers with unique position were used to estimate the LD of B. carinata. The parameter r² (the squared Pearson correlation coefficient) between all pairs of SNP markers was used to estimate LD by the software package TASSEL 5.0 (Bradbury et al., 2007).

The genome-wide association analysis with four statistical models was performed by TASSEL version 5.0 (Bradbury et al., 2007). The four models are as follows: the general linear model (GLM) included a naïve model without controlling for population structure; the Q model which controlled for population structure; the mixed linear model (MLM) including kinship and the Q+K model which controlled for both population structure and kinship. The threshold for the significance of associations between markers and traits, P < 8.32 × 10⁻⁵ (P = 1/total markers used; −log₁₀ (1/12,030) = 4.08) as in Li et al. (2016) was used in this study.

The field experiments with 185 DH lines of the YW population and two parental lines were conducted in randomized complete block designs, with three replicates except for 2010 in Wuhan where only two replicates were followed. The trials were conducted in three contrasting environments: (i) spring [two experiments, one each at Hezheng (Gansu province, China) and Xining (Qinghai Province, China)], (ii) semi-winter at Wuhan, representing the major rapeseed growing zone in China, and (iii) spring at Wagga Wagga (New South Wales, Australia; Supplementary Table 2, Supplementary Figure 1). The 81 diverse accessions were grown in Wuhan in 2013 and 2014 in randomized complete block designs with three replicates (Supplementary Table 2, Supplementary Figure 1). For Chinese environments, each plot was 3.0 m² and 30 plants were planted in three rows with three replicates, with a distance of 40 cm between rows and 25 cm between individual plants. For Australian environment, trials were sown in pots under the birdcage at Wagga Wagga Agricultural Institute, NSW, Australia in 2013 and 2014. The trials consisted of two replications of four rows and 100 column pot array with each replication consisting of two rows and 100 column pot array. Five plants were sown per pot.

Seed yield, seven yield related traits (flowering time, pod width, pod length, seed number per pod, seed weight, pod number on main inflorescence, and length of main inflorescence) and six seed quality traits (protein content, oil content, erucic acid, linolenic acid, linoleic acid, and oleic acid) were investigated for QTL mapping in this study. The six seed quality traits (protein content, oil content, erucic acid, linolenic acid, linoleic acid, and oleic acid) were also measured on 81 diverse accessions of B. carinata that were used for GWAS. The seed yield and related traits were measured as described by Shi et al. (2009). Genetic variation in seed and oil quality traits was determined by Near-Infrared Reflectance Spectroscopy (NIR) method (Gan et al., 2003). Compared to other traits investigated in this study, only flowering time was determined in three different contrasting environments. We reanalyzed QTL for flowering time that were determined in a previous study (Zou et al., 2014), with the data collected from Australian Spring environment (WW13 and WW14).

The broad sense heritability was calculated as h² = σ_g²/(σ_g² + σ_ge²/n + σ_e²/nr), where σ_g² is the genetic variance, σ_ge² is the variance representing genotype by environment interactions and σ_e² is the error variance, n is the number of environments and r is the number of replications. Analysis of variance (ANOVA) was performed using σ_g², σ_ge², and σ_e² estimated by using a Proc general linear model (GLM) in SAS software (SAS Institute Inc., 2000). Pearson correlations were calculated between phenotypic traits of interest (Weaver and Wuensch, 2013).

A dense genetic linkage map of the YW DH population based on 4,031 presence-absence DArTseq markers corresponding to 1,366 unique loci (Zou et al., 2014) was utilized in this study. A total of 772 single nucleotide polymorphism (SNP) markers generated from the DArTseq were added to the genetic map by integrating in genetic bins (Supplementary Table 3). The present integrated map covered a genetic distance of 2,084 cM, with an average distance of 1.53 cM between loci and was subsequently used for QTL analysis.

QTL associated with various traits of interest were identified by the composite interval mapping model (Wang et al., 2012) using the software WinQTL cartographer version 2.5 (http://statgen.ncsu.edu/qtlcart/WQTLCart.htm). The criteria for identification of QTL were followed as described by Shi et al. (2009). The significant QTL (P = 0.05) with overlapping suggestive QTL (P = 0.5) named as “identified” QTL. The term “consensus” QTL was used for the same genomic interval identified for the same trait detected in different environments and was integrated by meta-analysis using the software BioMercator v4.2 (Goffinet and Gerber, 2000). We designated the “major” QTL which had R² ≥ 20% and minor QTL which had R² < 19.9%, as suggested by Collard et al. (2005). QTL that could be detected in more than one experiment were classified as “reproducible” QTL, while others were classified as “non-reproducible.”

We chose the reference genome of B. napus “Darmor-bzh” (Chalhoub et al., 2014), and the Tapidor/Ningyou (TN, B. napus) DH population as reference to compare QTL with those identified in YW DH population in this study. TN DH population has been extensively investigated previously and QTL associated with seed yield and its related traits and seed quality traits have been identified across different environments using genetic maps based on traditional makers (Shi et al., 2009) and a newly constructed SNP-based genetic bin map (Luo et al., under review). DArTseq markers underlying the identified QTL regions that were identified on the C genomes of YW DH population and TN DH population were aligned to the reference genome of B. napus “Darmor-bzh” (Chalhoub et al., 2014) by Blast using an e-value threshold of e⁻¹⁰ (Altschul et al., 1990) with ≥62 bp match length. The details of gene annotation for “Darmor-bzh” were cited from Körber et al. (2016). When the QTL regions of YW DH population could be aligned to the candidate gene for the same trait in the C genome of B. napus, we presumed there is a homologous candidate gene in the QTL region of B. carinata. When the QTL regions of YW DH population and TN DH population showed alignment to the same physical position on the reference genomes, we assumed the homologous QTL control genetic variation for trait of interest in YW and TN populations.

A total of 54,510 polymorphic DArT-seq markers (38,396 presence-absence, also called silico DArTs and 16,114 SNPs) were identified in 81 accessions of B. carinata. After filtering, we selected 33,924 high quality markers that had overall call rates >90% and reproducibility >0.9, missing data of <10%, and MAF (minor allele frequency) >0.05 in the population, for the analyses of population structure and genetic diversity. The LnP(D) value decreased continuously with the change of K from 1 to 10, and the most significant change was observed when K increased from 1 to 2, ΔK also showed a peak at K = 2 (Figure 1, Supplementary Table 4). Accordingly, the 81 accessions could be divided into two major sub-populations (Figure 2A). Both parental lines of the YW DH population were found to represent different lineage as revealed by both population structure and NJ phylogenetic tree (Figure 2).

When we aligned the marker sequences (33,924) of the diverse B. carinata lines to the published Cⁿ genome sequence of B. napus (Chalhoub et al., 2014) and the C° genome sequence of B. oleracea (Parkin et al., 2014), respectively, much more B. carinata markers sequence (1,458) could be uniquely aligned to B. oleracea (5,126) than that to B. napus (3,668) under the same parameters. Besides, the quality of the alignment to B. oleracea were much better, 87.87% of the alignments to the B. oleracea with over 65 bp match length, which is much higher than aligned to B. napus (0.65%). To avoid the potential subgenomic differentiation of the C genome between B. carinata and B. napus, we used the aligned results against its progenitor B. oleracea, as C genome reference. As a result, a total of 10,999 markers that uniquely aligned to the genomes of B. nigra (B genome) and B. oleracea (C genome) were then used to estimate the LD pattern of the diverse lines. Our result indicated that the C genome showed a longer LD decay (400 Kb) than the B genome (250 Kb) (Figure 3, Table 1). The LD decay distance varied across 17 chromosome, from 150 Kb (B4) to 650 Kb (C4, C6, and C9; Table 1). These results revealed significant differences in the level of LD between different chromosomes and subgenomes.

Estimated means of parental lines, ranges and h²-values of traits measured in DH lines are given in Table 3. Both parental lines Y-BcDH64 and W-BcDH76 differed significantly for most traits measured across phenotypic environments (Table 2). W-BcDH76 had higher seed yield and 1,000 seed weight and early flowering time than Y-BcDH64. Estimated means for different seed quality traits in Y-BcDH64 were higher than W-BcDH76 in most environments except for oil content and oleic acid. Erucic acid was significantly different in the contrasting environments. Y-BcDH64 had higher erucic acid in semi-winter environment, but had a lower content in spring environment (Table 2).

A wide range of phenotypic variation was observed for all 14 traits between the parents and among the DH lines under different environments. Moderate to high h² was observed in the YW population for 10 traits measured in more than one experiment, ranging from 15.59% for pod width (spring) to 90.79% for linolenic acid (spring) (Table 3). Frequency distributions for trait means clearly showed transgressive segregation for all the traits in the DH lines (Figure 4, Supplementary Figure 2). ANOVA on the investigated traits across different environments showed that “genotype (G),” “environment (E),” and “G × E” have significant effects for almost all of the traits except for the flowering time (at Wagga Wagga, Australia; Table 3). However, significant positive correlations for the flowering time were observed in different experiments, WH09 and WH10, HZ12, and XN12, WW13, and WW14 (Figure 5A). The correlation coefficient for the flowering time observed between semi-winter (China) and spring (Australia) experiments was higher than that between semi-winter (China) and spring (China), and spring (China) and spring (Australia).

Phenotypic correlations among seed quality traits were also calculated (Figure 5B). The correlation between oil content and oleic acid, and oil content and linoleic acid showed the significant negative correlation, but a significant positive correlation was observed between oil content and erucic acid, as well as oil content and linolenic acid. Oleic acid showed positive correlation with linoleic acid, but significant negative correlation with linolenic acid and erucic acid. These results showed that the different fatty acid components in the seeds correlate with high erucic acid.

A total of 150 identified QTL for seed yield, flowering time, pod width, pod length, seed number per pod, seed weight, pod number of main inflorescence, and length of main inflorescence were identified in the YW DH population (Supplementary Table 5). After binning QTL which were repeatedly detected from different experiments for the same traits, a total of 116 consensus QTL were identified for the eight agronomic traits. Of these, 76 consensus QTL were detected in B genome while 40 were detected in C genome. The most number of consensus QTL were identified in chromosomes B1 and B8 (13 each), while the least number of consensus QTL (1) were identified in C1 (Figure 6, Supplementary Figure 3).

Six consensus QTL for seed yield were identified (only measured in two experiments, HZ12 and XN12, because of the high pod shattering) on chromosomes B2, B4, B6, C4, and C9, explaining phenotypic variation ranging from 4.2 to 16.7%. The parent W-BcDH76 contributed favorable alleles of the QTL on B2, B4, C4, and C9, whereas parent Y-BcDH64 contributed favorable alleles of two consensus QTL for seed yield (qSY.B6-1, qSY.B6-2) located on B6 chromosome (Supplementary Table 5). None of the consensus QTL was consistently detected across experiments. But four of the six consensus QTL overlapped with the QTL accounting for other yield related traits. For instance, qSY.B4-1 and qSY.C9-1 for seed yield were overlapped with the QTL for seed number per pod (qSN.B4-1, qSN.C9-1), and qSY.B6-1 was overlapped with the QTL for seed weight (qSW.B6-1) (Supplementary Figure 3).

For flowering time, a total of 28 consensus QTL were identified on 11 chromosomes in three macro-environments (Australian spring, Chinese spring, and Chinese semi-winter; Figure 7). The identified QTL (qFT35, C6) integrated in the consensus QTL could explain up to 59.40% of the phenotypic variation in flowering time (Supplementary Table 5). A total of three consensus QTL, all of which located in C6, were identified as “major QTL” individually explaining phenotypic variation over 20% (Table 4). For instance, the consensus QTL, qFT.C6-3 on chromosome C6 within 1 cM interval (genetic position from 35.5 to 36.5 cM) explained the maximum phenotypic variation (from 28.0 to 59.4%). Of the three major QTL, qFT.C6-1 was identified in the semi-winter environment of China with high R² (27.13% in experiment WH09 and 42.36% in experiment WH10), but low R² (2.35% in experiment XN12) identified in the spring environment of China. Among the 28 total QTL for flowering time, seven of them (qFT.B4-1, qFT.B4-2, qFT.B4-3, qFT.C6-1, qFT.C6-3, qFT.C8-1, and qFT.C8-2) were identified in at least two macro-environments (Supplementary Table 5). Two of the seven QTL (qFT.B4-1 and qFT.B4-2) were identified in all the three environments, and both were located in B4 chromosome. Four of the seven consensus QTL, qFT.B4-3, qFT.C8-2, and qFT.C6-1, qFT.C8-1 were identified in the two spring macro-environments (Australian spring and Chinese spring) and two Chinese macro-environments (Chinese semi-winter and Chinese spring), respectively. While, the last one of the seven consensus QTL, qFT.C6-3 was identified in both of the Chinese semi-winter environment and Australian spring environment. These results indicated the differential genetic effects of the QTL in response to environment.

A total of 296 identified QTL were detected for protein content, oil content, erucic acid, linolenic acid, linoleic acid and oleic acid in the seeds of the YW DH population (Supplementary Table 5). These identified QTL could be integrated as 166 consensus QTL and localized across the 17 chromosomes, explaining phenotypic variation varying from 0.2 to 50.4%. Four of these QTL, located in C4, explained ≥20% of the phenotypic variation and therefore were designated as major QTL (Table 4). Consistent with the distribution of QTL for seed yield and yield related traits in subgenomes, more QTL were identified in the B genome (92) compared with that in the C genome (74) for seed quality traits (Figure 6). In chromosome B8, the most number of consensus QTL (24) were detected for seed quality traits, whereas the least (1) were identified in chromosome C7.

Among the six seed quality traits, the maximum number of consensus QTL were identified for oil content (34) and the minimum were identified for protein content (20). The QTL for oil content were detected on all chromosomes of B. carinata genome except B9 and C7 (Supplementary Table 5). The R² explained by these QTL ranged from 2.67% (qOC.C5-1 on C5) to 10.4–23.4% (qOC.C4-1 on C4). Nine consensus QTL for oil content were repeatedly detected in at least two experiments. The alleles from the parent W-BcDH76, at the QTL located in chromosome B5, C3, C4, C5, C6, C8, and C9, contributed to the increase of oil content in YWDH population.

A total of 30 consensus QTL for erucic acid were detected on all chromosomes of B. carinata except C2 and C7 (Supplementary Table 5). The QTL qERU.C4-1 on C4, accounted for the maximum phenotypic variation (5.7–32.6%) for erucic acid. The alleles from the parent Y-BcDH64 for the four QTL (qERU.B5-1, qERU.C3-2, qERU.C9-1, and qERU.C9-2) contributed negatively to increase the erucic acid.

Using the phenotypic data of six seed quality traits investigated for the 81 diverse lines in two experiments, we also identified the markers accounting for genetic variation in seed quality traits by GWAS and compared with the results revealed by linkage analyses. A total of 12,030 DArT-seq markers with known physical position were used for GWAS. According to the QQ plots, the suitable model for each trait was selected. Seven markers were found to be significantly associated with five traits (Table 5). Linoleic acid showed a very limited phenotypic variation among diverse lines. As a result, we could not detect any association for this trait. One of the seven markers (5121285) was significantly associated for three traits (PRO-WH13, OC-WH13, and ERU-WH13). Two markers 5868483-1S and 5121285 showed significant association with linolenic acid (LEN-WH14) and erucic acid (ERU-WH13), respectively. Both these markers also located in the QTL region for the same traits in the YW DH population (Table 5). The marker 5,121,285 was also aligned to the homologous QTL region for oil content of the TN DH population of B. napus (Luo et al., under review).

A majority (78%) of the consensus QTL accounting for different traits in the YW DH population investigated could be overlapped with each other (Supplementary Figure 3, Supplementary Table 5). For example, one of the QTL accounting for seed weight (qSW.B4-1) was overlapped with the one for flowering time (qFT.B4-1) in B4 chromosome. Whereas, the other QTL for seed weight (qSW.B5-2) was overlapped with the one for pod width (qPW.B5-1) in B5 chromosome. In C4 chromosome, there were five overlapping QTL, qPRO.C4-1, qLEI.C4-2, qOLE.C4-2, qLEN.C4-1, and qOC.C4-1 for protein content, linoleic acid, oleic acid, linolenic acid, and oil content, respectively (Supplementary Figure 3, Supplementary Table 5).

To determine whether the homologous QTL genomic regions of B. carinata controlling variation for yield related traits and seed quality traits present in the C genome could also modulate genetic variation for these traits in the other oilseed Brassica crop containing C genome, we compared and aligned the QTL regions of YW DH population to the physical position of those QTL regions of TN DH population and candidate genes predicted in B. napus (the reference genome Darmor-bzh). The QTL accounting for flowering time, erucic acid and oleic acid of B. carinata could be aligned to the candidate genes of the B. napus for the three traits (Figure 8). For example, the genetic intervals of three identified QTL for flowering time (qFT34, qFT35 and qFT36) were aligned to the physical positions 28.55–29.42 Mb on chromosome C6 of the reference genome of Darmor-bzh (Figure 8A). There were two genes controlling flowering time (BnC6.FT.a and BnC6.FT.b; Wang et al., 2009) located in the 28.55 Mb regions on C6 chromosome of Darmor-bzh. This suggested that both homoeologues of FLOWERING LOCUS T gene are likely candidates for flowering time variation in YW DH population of B. carinata. The QTL accounting for erucic acid (qERU37) in C3 chromosome of YW DH population was aligned to the gene BnaC.FAE1 of B. napus ((Körber et al., 2016); Figure 8B). The QTL accounting for oleic acid in C3 (qOLE25) and C8 (qOLE47, qOLE48, and qOLE49) of the YW DH population were aligned to the gene FAR7 gene in C3 and C8 of oleic acid in B. napus (Körber et al., 2016), respectively.

Additionally, a total of 17 QTL detected in the C genome of the YW DH population could also be aligned to the homologous QTL region for the same traits in the C genome of B. napus (Supplementary Table 5). For example, the identified QTL (qFT34, qFT35, and qFT36) for flowering time was aligned to the homologous QTL region in the TN DH population (Luo et al., under review). Three QTL (qSW26, qSW27, and qSW31) for seed weight and one QTL (qPRO28) for protein content were also aligned to the homologous QTL region of TN DH population. For the erucic acid, ten of the identified QTL (qERU37, qERU49 to qERU56, qERU59) detected in C3, C6, and C8 chromosomes could be aligned to the homologous QTL region in the C genome of B. napus.

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest. The reviewer EF and handling Editor declared their shared affiliation, and the handling Editor states that the process nevertheless met the standards of a fair and objective review.