Genetic dissection of seed oil and protein content and identification of networks associated with oil content in Brassica napus

High-density linkage maps can improve the precision of QTL localization. A high-density SNP-based linkage map containing 3207 markers covering 3072.7 cM of the Brassica napus genome was constructed in the KenC-8 × N53-2 (KNDH) population. A total of 67 and 38 QTLs for seed oil and protein content were identified with an average confidence interval of 5.26 and 4.38 cM, which could explain up to 22.24% and 27.48% of the phenotypic variation, respectively. Thirty-eight associated genomic regions from BSA overlapped with and/or narrowed the SOC-QTLs, further confirming the QTL mapping results based on the high-density linkage map. Potential candidates related to acyl-lipid and seed storage underlying SOC and SPC, respectively, were identified and analyzed, among which six were checked and showed expression differences between the two parents during different embryonic developmental periods. A large primary carbohydrate pathway based on potential candidates underlying SOC- and SPC-QTLs, and interaction networks based on potential candidates underlying SOC-QTLs, was constructed to dissect the complex mechanism based on metabolic and gene regulatory features, respectively. Accurate QTL mapping and potential candidates identified based on high-density linkage map and BSA analyses provide new insights into the complex genetic mechanism of oil and protein accumulation in the seeds of rapeseed.

Scientific RepoRts | 7:46295 | DOI: 10.1038/srep46295 and Brassica [21][22][23] . QTLs with pleiotropic effects for SOC and SPC have also been observed in Brassica species 15,24 . Identifying the underlying candidate genes and recognizing the pleiotropic effect or correlation between these two traits would greatly increase the breeding efficiency.
Nevertheless, in most of these studies, the QTLs spanned large ranges of the genome due to low mapping resolutions, and generally there was a lack of marker sequence information, hindering the integration of these QTLs into the reference physical map and the identification of the underlying candidate genes. In addition, the genomic sequences of the parental cultivars were also unknown in most of the previous studies, which also complicated the identification and analysis of candidate genes. Recent innovations in genome sequencing technology and bioinformatics and single nucleotide polymorphism (SNP) markers, especially the 60 K SNP Infinium Array, have been developed and widely applied for QTL mapping and the genetic dissection of important agronomical traits in B. napus [25][26][27][28][29][30][31] . Simultaneously, the reference genome of B. napus has been released 32 , making it feasible to perform a comparative genome analysis between the genetic map and the physical map using SNPs with sequence information. Additionally, bulk segregant analysis (BSA) combined with next-genomic sequencing technologies has been used as a fast track approach to more rapidly identify candidate genomic regions. This approach involves the selection of 20-50 lines with extreme phenotypic values, with pooling at equivalent concentrations followed by sequencing of the pools for subsequent variance analysis 33 . BSA has been used successfully to map candidate genomic regions for early flowering in cucumber 34 , 100-seed weight in chickpea 35 , and fruit weight and locule number in tomato 36 .
QTL studies have led to the identification of numerous loci that are responsible for variation in SOC and SPC, but the corresponding genomic regions with their complex metabolism and regulatory mechanisms have not been studied with great accuracy. Therefore, candidate genes within QTL regions combined with the relevant metabolic pathway in B. napus have been not reported. There is extensive genome homology and microsynteny between Brassica species and Arabidopsis thaliana (A. thaliana) because of their close evolutionary relationship 37 . More attention has been focused on genes and regulatory factors involved in acyl-lipid metabolism, and more than 120 enzymatic reactions and at least 600 genes involved in acyl-lipid metabolism in A. thaliana have been well documented by Li-Beisson et al. 38 . A. thaliana acyl-lipid and storage protein-related orthologous genes mapped in a confidence interval (CI) of QTLs on the B. napus genome might be valuable to identify the candidate genes and dissect the complex metabolism and regulatory mechanism associated with SOC and SPC.
This present report consists of the following three parts: (1) the construction of a high-density genetic linkage map based upon SNP markers combined with non-SNP (SSR and STS) markers that have been used previously to generate a primary linkage map in the KenC-8 × N53-2 (KNDH) population 39 ; (2) precise QTL mapping of SOC and SPC across multiple environments based on a high-density genetic linkage map and rapid detection of associated genomic regions (AGRs) using BSA; (3) identification of candidate genes within QTL regions and provision of primary insights into the complex metabolism and regulatory mechanism of SOC and SPC in B. napus.

Phenotypic variation and correlation of SOC and SPC in the parents and the DH population.
Combined with the data that were previously reported by Wang et al. 39 , a total of 12 and 11 sets of SOC and SPC trial data were used to analyze the phenotypic variation (PV), respectively (Table S1). The SOC and SPC of the two parents (N53-2 and KenC-8) differed significantly in all investigated environments. The SOC of N53-2 was 10% higher than that of Ken-C8 in all microenvironments except 10HG. While the SPC of Ken-C8 was nearly 5% The blocks with green bands at the outermost circle represent the 19 genetic linkage groups and the marker distribution. The third circle (from the outermost one) represents the 19 chromosomes of B. napus. The lines connecting them represent their col-linearity relationship. The 12 inner circles with 3 colors represent 12 microenvironments in winter-type, spring-type and semi winter-type rapeseed growing areas, respectively. The short bars with a red and blue color within the 12 inner circles represent QTLs identified in different environments and linkage groups, and they are located between the outermost two circles representing consensus QTLs.
higher than that of N53-2 in all semi-winter environments, they also showed significant differences, though to a reduced extent for the winter environments. The KN population exhibited a normal or near-normal distribution and transgressive segregation with a small degree of both traits (Fig. 2). These results indicated that quantitative inheritance was suitable for QTL identification, with polygenic effects and favorable alleles mainly inherited from one of the two parents. Further analysis revealed that both parents and populations in different trials all showed significant differences (Table S1), suggesting that both traits were significantly influenced by the environment. In addition, SOC and SPC were significantly negatively correlated with average coefficients of − 0.65; the highest correlation up to − 0.77 was identified in the 11GS (Table S2), suggesting competition among sinks for assimilates.

Identification, integration and co-localization analysis of QTLs for SOC and SPC. A total of 164
identified QTLs were observed for SOC in 12 trials, with an average CI of 7.43 cM. These QTLs were distributed on 12 chromosomes, excluding A04, A05, A06, A07, C01, C04 and C07, and a single QTL could explain up to 22.24% of the PV. The largest number of QTLs (136) was distributed on A03, A08, A09, C03, C05 and C06, accounting for 82.93% of the total number ( Fig. 1 and Table S3). The 164 identified SOC-QTLs were then integrated into 67 consensus QTLs, and the average CI was reduced from 7.36 to 5.26 cM (Tables 2 and S3 and Fig. 1). Fifty-three consensus QTLs were distributed on A03, A08, A09, C03, C05 and C06, accounting for 79.10% of all consensus QTLs, and 35 of them (76.09%) showed environmental stability. There were 59, 30, and 16 QTLs detected in winter type, semi-winter type and spring type environments, respectively ( Figure S1a). Only one QTL (cqOC-A9-9) could be detected in all three macroenvironments. Eleven, two and one QTL could be expressed stably at least over 2 years in winter, semi-winter and spring type environments, respectively. Six major QTLs (cqOC-A2-3, cqOC-A9-9, cqOC-A9-10, cqOC-C5-3, cqOC-C5-4, and cqOC-C6-6) were repeatedly detected in at least two trials and displayed a large effect (PV > 10%), among which cqOC-A9-9 and cqOC-C6-6 were detected stably in six microenvironments, while cqOC-C6-6 explained the highest PV (22.24%). Although accounting for less than 10% of the PV value, cqOC-C3-4 could be detected simultaneously in eight microenvironments and was considered a major QTL that merits attention. A total of 68 identified QTLs for SPC were observed in 11 microenvironments, which were distributed on 12 chromosomes and accounted for 2.19-27.48% of PV with an average CI of 9.31 cM ( Fig. 1 and Table S3). The same results were obtained for SPC-QTLs, most of which (64.71%) were distributed on A03, A09, C03 and C05, with the difference that none of the SPC-QTLs were detected on A08. The identified SPC-QTLs were then integrated into 38 consensus QTLs with a narrower average CI of 4.38 cM (Table S3 and Fig. 1). Twenty-three QTLs could be detected in at least two microenvironments, and cqPC-C5-6 could be identified in four microenvironments, while five QTLs (cqPC-A3-4, cqPC-A9-2, cqPC-C1-2, cqPC-C5-5 and cqPC-C5-7) were identified in three microenvironments. None of the consensus QTLs could be detected in all three macroenvironments, and six and one QTL could be expressed stably at least over 2 years in winter and semi-winter type environments, respectively (Fig. S1b). Four QTLs, cqPC-A2 and cqPC-C1-2, which were repeatedly detected in at least two trials with a medial effect (20% > PV > 10%), and cqPC-A3-2 and cqPC-A3-5, which were detected in one trial but had a large effect (PV > 20%), were also treated as major QTLs.
In addition, as observed from the additive effect (Table 2), the N53-2 alleles in the majority of SOC-QTLs (81.25%) increased SOC, but most SPC-QTLs (81.58%) decreased SPC, consistent with the higher SOC and lower SPC of N53-2 compared with Ken-C8 (Table S1).
The close genetic relationship between SOC and SPC could be confirmed by mapping the co-localized loci. All 11 co-localized QTLs affected both traits and showed opposite additive effects (Table 3), and the favorable alleles of the 10 co-localized QTLs that were derived from the high oil content parent N53-2 could increase SOC and decrease SPC, with only uqC5-1 demonstrating the opposite. Close co-localization explained the high genetically negative correlation between SOC and SPC, and suggested a complexity of the competitive mechanism of these two traits.
Additionally, some QTL clusters associated with both traits distributed on chromosomes A09, C03 and C05 were also detected. For example, one of these clusters falling in interval from 287.50-322.64 cM on C03 influenced SOC and SPC with consistently negative effects ( Figure S2). These results indicated that tight linkage was also an important factor in addition to pleiotropy for their competitive relationship. In addition, there were SOC-QTL clusters, but none of the SPC-QTLs were found on A08, which suggested that more allelic variations control SOC compared with SPC on chromosome A08 between two parents.

BSA analysis based on extreme bulks and confirmation of QTL mapping. Twenty-four DH lines
with high SOC (47.2-52.7%) and 24 with low SOC (38.0-42.4%) were selected to prepare two extreme bulks. As a significant negative correlation of SOC and SPC, high and low SOC bulks corresponded to low (19.8-23.9%) and high (22.5-27.2%) SPC, respectively.
Illumina high-throughput sequencing was used to generate 274.08 million and 274.85 million reads from the low and the high SOC bulks, with a coverage of 84.30% (35.85-fold) and 90.74% (39.00-fold) of the reference genome of B. napus, respectively. In contrast, 97.47% (17.16-fold) and 96.56% (19.82-fold) genome coverage was obtained for the two parents, KenC-8 (low SOC and high SPC) and N53-2 (high SOC and low SPC), confirming the accuracy of the subsequent analysis (Table S4). The graph of Δ (SNP-index) was generated ( Figure S3) using the genome sequence of KenC-8 as a reference.
One hundred and thirty-three AGRs underlying SOC were identified on 17 chromosomes (apart from A01 and C07) (Table S5), and 38 AGRs overlapped and/or narrowed the SOC-QTLs from QTL mapping based on the present genetic linkage map (Table S6). The major QTL cqOC-A9-9 was validated by a significant region AGR_A9-8, and the major QTL cqOC-A9-10 was divided into three narrower AGRs (AGR_A9-9, AGR_A9-10 and AGR_A9-11). It is noteworthy that some AGRs with a significantly high Δ (SNP-index) and QTL hotspots were found to co-localize on chromosome A03, A09, C03, and C05 (Table S6). For example, the associated genomic region (AGR_C3-16) was identified to overlap with four SOC-QTLs (cqOC-C3-4, cqOC-C3-5, cqOC-C3-6 and cqOC-C3-7) that could be stably expressed in at least two microenvironments, among which cqOC-C3-4 could be assessed stably in eight microenvironments (Fig. 3). Identification and analysis of potential candidates related to acyl-lipid metabolism and seed storage underlying SOC and SPC, respectively. Accurate QTLs for SOC were identified using a high-density linkage map and confirmed by BSA. Through comparative mapping between the linkage map and the physical maps of B. napus based on SNP sequence information, a satisfactory colinearity relationship between the genetic linkage map and the reference genome of B. napus was validated ( Fig. 1), which could help us to accurately identify genomic regions and the corresponding candidate genes underlying these QTLs. For example, SOC-QTL cqOC-A9-9 could be stably identified in six microenvironments, explaining up to 14.94% of the PV with 0.98 cM of CI, which was mapped within a 0.07 Mb region in the reference genome (Table S7). Finally, 59 genomic regions of SOC-QTLs (88.06%) and 29 genomic regions of SPC-QTLs (76.32%) could be determined accurately (Table S7).
A total of 448 genes involved in 21 of 22 acyl-lipid pathways were identified within 41 SOC-QTL regions (Figures S4 and S5, and Table S8). Additionally, a whole genome search was implemented to determine the distribution of orthologues related to acyl-lipid metabolism, and a total of 3058 orthologues of A. thaliana acyl-related genes were identified in the whole genome ( Figures S4 and S5). Minor acyl-lipid related gene enrichment / prevalence in the CI of SOC-QTLs was found compared to the whole genome (Table S9). For example, six orthologues of AtLACS9 that encoded long-chain Acyl-CoA synthetase (plastidial) were distributed on the chromosome of the B. napus genome, four of which displayed SNP/InDel variation in the SOC-QTL regions. Moreover, four of eight orthologues of AtFAD6 that encoded oleate desaturase fell within SOC-QTL regions and exhibited SNP/ InDel variation. According to the re-sequencing results of the two parents from BSA, 355 of 449 acyl-lipid related orthologues had SNP/InDel variation between the two parents (Table S10). These genes were involved in 21 of 22 acyl-lipid metabolism pathways. The 25, 13, 24, 20 and 9 genes among them were involved in FA synthesis in plastids, FA elongation, the synthesis of triacylglycerol (TAG), TAG degradation and beta-oxidation, respectively, accounting for 24.20% of the total potential candidates associated with lipid metabolism within SOC-QTL regions ( Figure S6 and Table S10). Additionally, a large proportion of genes involved in lipid signaling and cuticular wax syntheses were mapped, accounting for 11.0% and 19.0% of the total acyl-lipid orthologous genes. This finding supported the complexity of the oil accumulation, which was also likely regulated by some genes related to lipid signaling, and the cuticular wax biosynthesis pathway probably also played a potential role in the accumulation of oil during embryonic development.
Two and eleven potential candidates related to oil body structural protein and seed storage protein were identified in 12 SPC-QTL regions (Table S11). BnaA09g02110D, annotated as an orthologue of oleosin 4, was located in an overlapping region of cqPC-A9-2 and cqOC-A9-2. In addition, seven seed storage protein-related orthologues (BnaA09g13220D, BnaA09g11520D, BnaC03g65080D, BnaC05g02160D, BnaC05g44560D, BnaC05g44510D and BnaC06g18820D) co-localized within the SOC and SPC-QTL interval. It is also noteworthy that some important potential candidates related to acyl-lipid metabolism were identified in overlapping regions of SOC-QTLs and AGRs of BSA (Table S10). BnaA09g48250D, a homologue of AtCAC2 involved in the key steps of de novo FA synthesis, was identified in AGR_A9-12 and the major QTL cqOC-A9-10, and it contained a nonsynonymous SNP in an exon. BnaA08g12780D, corresponding to AtFAD6 encoding oleate desaturase, was located in AGR_A8-4 and cqOC-A8-3. Twenty-two potential candidates related to acyl-lipid metabolism were identified within three major QTLs overlapping with AGRs: cqOC-A9-9 (1), cqOC-A9-10 (14) and cqOC-C5-4(7) (Table S10). Three genes (BnaA09g48250D, BnaC05g21140D and BnaC05g24510D) involved in plastidial FA synthesis, and two genes (BnaA09g51530D and BnaC05g21070D) involved in FA elongation were identified in the these regions.
Six potential candidates (BnaA08g12780D, BnaA09g48250D, BnaC03g65980D, BnaC03g67820D, BnaC05g43390D and BnaC05g44510D) that were reported to be involved in acyl-lipid metabolism in B. napus or A. thaliana were selected to measure the expression levels in the developing seed at 15, 25, and 35 days after flowering from both parents. All six genes showed significant expression difference between parents during different developmental periods (Fig. 4). The relative expression level of BnaA09g48250D, which was homologous to the AtCAC2 gene encoding BC (subunit of ACCase), from the overlapping region of AGR_A9-12 and the major QTL cqOC-A9-10 increased significantly in both parents and were clearly higher in N53-2 during the middle and late embryonic development periods. BnaA08g12780D and BnaC03g67820D, which are homologous to AtFAD6 (oleata desaturase), were determined from cqOC-A8-3 and cqOC-C3-7, their relative expression level was reduced but higher in N53-2 during the late developmental period. BnaC03g65980D encoded FAE1 and was involved in FA elongation, showing a rapid increase in expression throughout the process of embryonic development and a significantly higher level in N53-2 than in KenC-8. Of interest, the relative expression level of BnaC05g44510D located in cqOC-C5-7 and involved in cuticular wax synthesis was clearly higher in N53-2 than in KenC-8 but decreased rapidly and was lower in N53-2 during the late embryonic development period. This phenomenon might be associated with the usage of more precursors for the synthesis of TAG in N53-2, especially during late embryonic development.
Genetic dissection based on the metabolic pathway and interaction network using potential candidates. To investigate the genetic basis of oil and protein accumulation in the seeds of B. napus, a magnificent metabolic pathway was constructed with pyruvate as the centric junction based on potential candidates underlying SOC and SPC-QTLs (Fig. 5). The pathway was involved in almost all primary carbohydrate metabolism, including starch and sucrose metabolism, the pentose phosphate pathway, glycolysis/gluconeogenesis, the Calvin cycle, the TCA cycle, FA synthesis, acyl-lipid metabolism, cutin, suberin and wax biosynthesis and amino acid metabolism. There were 167 and 98 orthologues identified in the SOC-and SPC-QTL regions, respectively, which were involved in this enormous pathway, 58 of which could be identified in both traits. IPGAM2, which encodes the 2,3-bisphosphoglycerate-independent phosphoglycerate mutase that catalyzes the formation of glycerate-2P 40 , was identified in overlapping CIs of qOC-C5-11 for SOC and the QTL qPC-C5-5 for SPC. During glycolysis/gluconeogenesis, some orthologues of key rate-limiting enzymes, such as HKL, PFK, GAPCP and PGK, which are involved in a series of catalyzing the conversion of glucose into pyruvate, were candidates in the co-localized CIs or major QTLs. The Calvin cycle carries out carbon fixation and is an important process for photosynthesis. The TCA cycle is of central importance to many biochemical pathways, and it provides the precursors of certain amino acids and the NADH used in numerous other biochemical reactions through the oxidation of acetyl-CoA derived from carbohydrates, oil and proteins. More than half of the orthologues involved in these two key metabolic pathways could be detected in SOC-and SPC-QTL regions. These results suggested that primary carbohydrate metabolism indeed plays a basic role of providing the initial material for protein synthesis and oil accumulation during embryonic development. To further dissect the mechanism of oil accumulation in seeds, the interaction of potential candidates from SOC-QTL regions was analyzed based on A. thaliana orthologues. The primary results revealed that the whole network incorporated 5817 nodes and 38887 edges ( Figure S7). Annotation and clustering elicited an interaction network associated with lipid metabolism (Fig. 6). Some biologic processes, such as protein, secondary metabolism, photosystems, amino acid metabolism, cell size or cycle, and signaling pathways, participated in the regulation of lipid metabolism through some pivotal genes, such as CAC2, SNC4, ACX2 and PLDGAMMA1, which could be expressed during the embryo stage of the plant 41 .

Discussion
Most studies investigating important agronomic traits in B. napus were conducted with lower-density genetic maps [6][7][8][9][10][11] and some with high-density linkage maps 25,26,29 , but QTL analyses of SOC and SPC simultaneously based on high-density linkage maps have not been reported. In this study, we constructed a high-density integrated genetic linkage map with 3106 SNP and 101 non-SNP (SSR and STS) markers covering 3072.7 cM of the KNDH population. The release of the B. napus genome and the application of SNP markers with known sequence information facilitated alignment of the genetic linkage and physical maps of B. napus. Thus, a high-density   39 . Here, 67 SOC-QTLs integrated from 164 identified QTLs were mapped, and ten new SOC-QTLs were contributed due to 4 new microenvironments ( Table 2). The main cause of the increased QTL number detected in the present map was the higher map density rather than the addition of new environmental data (the number of consensus QTLs increased from 24 to 67, only 10 of which were contributed by new environmental data). Molecular-assisted breeding is the primary objective for detecting genetic variation in important agronomic traits. This will help us to identify stable QTLs for molecular-assisted selection (MAS) to compare previously published QTLs for oil content in other segregating populations (with different genetic backgrounds) in our current study. This method has been reported in many studies 42,43 . To validate the consistency and availability of the QTLs detected using the present high-density map for MAS, we compared SOC-QTLs with the published QTL map for oil content in other segregating populations with different genetic backgrounds. Twenty-six SOC-QTLs were found to overlap with loci checked by other populations (Table S12). The major QTL, cqOC-A9-9, was also detected and overlapped with qOC-A9-4-TN in the TN population, which was used to identify the seed-oil content and erucic acid content QTLs 44 , and the QTLs for many other important agronomic traits [45][46][47] . Additionally, qOC-A3-4, which could be stably expressed in 8 microenvironments, could also be identified to the co-localized locus in the TN population. These QTLs are very stable and useful for further MAS analyses. It is interesting that some QTLs detected in other populations were identified and divided into multiple QTLs. For example, qOC-A8-1-TN and qOC-A8-RNSL from TN and RNSL, respectively, were overlapping, and both of them were found to cover 5 QTLs in the present study. This result also suggested that the present genetic linkage map had a better resolution than the maps used in previous reports.  The recently proposed NGS-based BSA approach is a cost-effective and rapid method of trait mapping supported by the identification of the target genomic region in rice 33 . A large number of AGRs were identified rapidly, and significant AGRs with a high Δ (SNP-index) were consistent with the SOC-QTL hotspots, confirming the accuracy of the QTL mapping based on the present map. In addition, we found that highly significant AGRs were subdivided into multiple QTLs (Table S6) rather than being subdivided or narrowed QTLs as previously reported 34,35,48,49 . These results suggested that the present high-density linkage map had sufficiently high resolution to dissect genetic variation in primary mapping populations.
The application of a high-density linkage map and satisfying linear relationship between the genetic and physical map in our study could be helpful to more accurately identify candidate genes underlying SOC and SPC. A large number of potential candidates involved in acyl-lipid pathways and seed storage were identified within SOC and SPC-QTL regions, respectively (Tables S10 and S11). Additionally, for the most important candidates of SOC, a whole genome search was conducted, followed by a search for enrichment/prevalence in the QTL confidence intervals compared to the whole genome. Minor acyl-lipid related gene enrichment/prevalence was found in the CIs of SOC-QTLs, which indicated that the important genes clustered around QTLs and were consistent with the frequently identified QTL hotspots in the present study. Six potential candidates related to acyl-lipid metabolism were selected as representatives to validate their expression levels between the two parents. The relative differences in the expression of five genes showed a relationship with the SOC phenotype of the two parents in developing seed. The results provided clues for exploring the casual genes underlying SOC. Single gene mutations may affect the expression of some genes in the common metabolic network. Consequently, further identification of the network is necessary, and simultaneously the rapid capture of causal genes requires the use of transcription and fine mapping methods 50 .
In the present study, the simultaneous identification of QTLs for SPC together with SOC using the high-density linkage map revealed that some QTLs affected both SOC and SPC in an opposite manner (Table 3). SOC and SPC, which have a negative correlation, are controlled by many genes with additive and epistatic effects 19,23,51,52 . Oil and protein are synthesized in the endoplasmic reticulum (ER) during the same period of seed development 53 , which suggests that these pathways interfere and/or compete with one another for ER resources. TAG synthesis competes with seed storage protein synthesis due to sucrose utilization, and the reduction of one product would result in higher levels of the other. Pyruvate, which can be converted into carbohydrates via gluconeogenesis, to FAs or energy through acetyl-CoA and amino acids, is a key intersection in the network of metabolic pathways and unites several key metabolic processes. Here, we constructed a magnificent primary carbohydrate pathway in the center of pyruvate using potential candidates underlying SOC-and SPC-QTLs (Fig. 5). A. thaliana orthologues involved in this pathway provide more information to identify the interference and/or competition mechanism from a metabolic perspective. In addition, some genes associated with cuticular waxes, cutin and suberin have also been observed, such as CER4, CYP86, ACC1 and MAH1, and they probably affect seed storage triacylglycerides through competition with the common fatty acid pool 54,55 . Interestingly, most orthologous genes involved in long chain fatty acid synthesis, especially C20:X and C22:X, fell within SOC-QTL regions. These results were in agreement with the findings of Gu et al. 56 , in which long chain FAs had a significant and positive relationship with oil content.  Table 3. Unique QTLs controlling both SOC and SPC.
Oil and protein biosynthesis during embryo development is a complex and hierarchical biological process that is regulated by some transcription factors and complex gene interaction networks 57,58 . Therefore, it is not inadequate to dissect the mechanism responsible for oil and protein accumulation in seeds only through primary carbohydrate metabolic pathways. In A. thaliana, metabolic pathways associated with the biosynthesis and degradation of acyl-lipids have been summarized, and most of the genes have been described 38 . The construction of an interaction network provided new insights into the interaction among related genes (Fig. 6). It has been reported that the yellow seeded Brassica genotypes contained greater oil, higher protein and lower fiber contents than the seeds of the black/brown seeded genotypes 59 . In plants, flavonoids are the major secondary metabolites involved in the Brassica seed coat pigmentation process 60 . The orthologues TT4 and TT5, which were classified into secondary metabolism pathways, were involved in flavonoid biosynthetic processes and regulated lipid metabolic pathways by interacting with CAC2 (Fig. 6). In addition, the result verified by the cross-talk between the lipids and potential flavonoid metabolism were driven by the finding that acetyl-CoA acts as a substrate for both pathways 61 . In addition, 4CL family proteins associated with lignin biosynthesis are linked to lipid metabolism through many nodes, which could help to explain the association with negative additive effects of QTLs for seed oil and lignin content (unpublished data). These results not only enhance our understanding of the complex regulatory mechanisms of oil accumulation but also serve as a basis for follow-up to identify causal loci together with transcriptome and genome-wide association analyses [62][63][64] .

Materials and Methods
Plant material and trait measurement. The segregating double haploid (DH) population, named KN, was used in this experiment. The KNDH population was derived from a cross between the parental lines KenC-8 with low oil content and N53-2 with high oil content, which was first constructed by Wang et al. 39 . A total of 300 lines of the KNDH population were implemented in the genotype and phenotype analyses in this study.
The field experiments were conducted in 3 macroenvironments as described by Wang et al. 39 : Dali of Shaanxi Province (DL), Wuhan and Huanggang of Hubei Province (WH and HG) and Sunan of Gansu Province (GS), which were winter, semi winter and spring-type rapeseed planting areas. Year-location combinations were treated as microenvironments. The materials were planted in Dali, Wuhan, Sunan and Huanggang for, successively, six, three, two and one year, respectively. Four types of microenvironmental data were added in this experiment apart based on a report by Wang et al. 39 . In total, there were 12 microenvironments. The SOC data from 2012 and 2013 in DL and WH, excluding the report by Wang et al. 39 , were collected. The SPC data were collected in all microenvironments except 08DL. The trial locations of WH and HG were the experiment bases of Huazhong University of Science and Technology (Wuhan), and DL and GS were the experiment bases of the Hybrid Rapeseed Research Center of Shaanxi Province. No specific permissions were required for the field trials. The field experiment was the same as reported by Wang et al. 39 .
Open-pollinated seeds were collected for the SOC and SPC measurements. The SOC data in the first eight trials (08DL, 09DL, 10DL, 10HG, 10GS, 11DL, 11GS, 11WH) were measured by nuclear magnetic resonance (NMR) according to Burns et al. (2003) with modifications 6 . The SOC data for the other four trials and all SPC data were measured by near-infrared reflectance spectroscopy as described by Gan et al. 65 .

SNP marker analysis.
Whole-genomic DNA was extracted from young seedlings using the CTAB method 66 . A total of 300 DH lines with two parental lines were successfully genotyped using the Brassica 60 K Infinium HD Assay SNP arrays according to the manufacturer's protocol provided by Illumina Inc. (http://www. illumina.com/). Imaging of the arrays after BeadChip washing and coating was performed using an Illumina HiSCAN scanner. Allele calling data for each probe was exported using Genome Studio genotyping software v2011 (Illumina Inc.). SNP data filtering was accomplished according to Zhang et al. 26 . Furthermore, SNPs with identical genotypes across the KNDH population were classified into a bin using Perl language based on a close linkage according to Zhang et al. 26 , and the first SNP of each bin was selected as its name.
Linkage map construction, QTL analysis and integration. The construction of genetic linkage groups, QTL analysis and integration were performed as described by Wang et al. 39 . Non-SNP (SSR, STS, and SRAP) markers were used to construct the primary linkage map used by Wang et al. 39 and also applied to construction of the map, together with the SNP-bin markers. Carthagene software was used to construct the high-density linkage map 67 . Windows QTL Cartographer 2.5 software was used to detect putative QTLs with a CIM 68 . For the CIM model, the scan walking speed was 2 cM, and the window size was 10 cM with five background cofactors. The LOD threshold (3.2-3.5) for the detection of significant QTLs was calculated by a 1,000-permutation test based on a 5% experiment-wise error rate. To avoid missing QTLs with small genetic effects, QTLs that appeared repeatedly in at least two microenvironments, below the LOD threshold (3.2-3.5) but above 2.0, were considered micro-real QTLs according to the definition of Long et al. 46 . All these QTLs, including significant and micro-real ones, were called identified QTLs 69 . The identified QTLs were named by combining the microenvironment and the chromosome number, e.g., qOC-10DL13-1. The identified QTLs detected in different trials with overlapping CIs were then integrated into consensus QTLs using the BioMercator V4.2 program 70 . If an identified QTL had no overlapping CIs with others, then it was also regarded as a consensus QTL. Consensus QTLs were designated by the initial letters "cq", such as cqOC-A1-1. The consensus QTLs with overlapping CIs for SOC and SPC were integrated into unique QTLs, and unique QTLs were named by the initial letter "uq", for example, uq-A9-1. The method used for QTL nomenclature and integration were as described by McCouch et al. 71 and Wang et al. 72  Comparative analysis of the genetic map and physical map, and the identification of genomic regions corresponding to CIs of QTLs and potential candidates underlying QTLs. The method used to identify the homologous locus of SNP probes in B. napus reference genomes was described according to Cai et al. 73 . The sequences of the 50-mer probes provided by Illumina Inc. were used as queries to search for homologous loci using the NCBI-Blastn local program against the B. napus "Darmor-bzh" reference genomes (http://www.genoscope.cns.fr/brassicanapus/data), as described by Altschul et al. 74 . The E-value threshold was set to 1e-10. All SNPs assigned to the linkage map were used to validate the alignment relationship of the genetic map and physical map. If an SNP was mapped to multiple loci in the B. napus reference genome, only the location that corresponded to the particular linkage group of the locus was selected 30 . Additionally, if a SNP had multiple homologous loci on the same chromosome, an accurate position for a particular locus was determined manually by referring to the physical positions of its upstream and downstream SNPs. Genome regions corresponding to the CIs of consensus QTLs, herein denoted QTL regions, were defined through both closely linked SNPs within QTL CIs and co-linearity between the linkage map and the reference genome. Genes that fell within QTL regions were regarded as candidates. According to the BSA re-sequencing, alleles that were within QTL regions and had SNP or InDel variations in introns, exons or within 1 kb up-and downstream between two parents were considered potential candidates, and otherwise were excluded.
BSA sequencing and AGR analysis. DNA was extracted from individuals in the KNDH population with extremely high and low SOC (24 DNA samples for each extreme trait) and bulked equally to generate high and low extreme bulks. Four libraries, including two parents (KenC-8 and N53-2), were constructed and then sequenced using the Illumina HiSeq2500 platform at the Novogene Bioinformatics Institute, and 125-bp paired-end reads were generated with an insert size of approximately 350 bp. Short reads obtained from both parents and two DNA bulks were aligned against the B. napus "Darmor-bzh" reference genomes to obtain the consensus sequence using BWA software 82 . Reads of high and low SOC bulks were separately aligned to the KenC-8 (low SOC parent) consensus sequence reads to call the SNPs. Variant calling was performed for all samples using the Unified Genotyper function in the GATK software 83 . SNPs were used in the Variant Filtration parameter in GATK. Homozygous SNPs between two parents were extracted, and the read depth information for homozygous SNPs in the above bulks was acquired to calculate the SNP Index 33 . The Δ (SNP-index) was calculated based on subtraction of the low SOC bulk from the high SOC bulk SNP-index. The average SNP-index of SNPs in a certain genomic interval was calculated using a sliding window analysis with a 1-Mb window size and 10-kb step size, and the SNP-index graphs and corresponding Δ (SNP-index) graphs were plotted. The genomic regions with a significant Δ (SNP-index) (P < 0.05) were identified as AGRs.
Quantitative PCR (qPCR) analysis of 6 acyl-lipid-related potential candidates. To investigate the potential candidates, we used qRT-PCR for six potential candidates related to acyl-lipids from different associated loci. RNAs were extracted and quantified from developing B. napus seeds at 15, 25 and 35 days after flowering (DAF), following the user manual for the RNA prep Pure Plant Kit (TOYOBO, DP441). cDNA was synthesized from 2 mg of total RNA using a cDNA Synthesis Kit. Expression analysis was performed with the SYBR premix EX TaqTM kit (TaKaRa, Japan) using an ABI 7900HT Fast Real-Time PCR System (Applied Biosystems, Framingham, USA). For each reaction, three technical replicates were evaluated. The primers of the six target genes and reference gene Actin are listed in Table S13.