Association Mapping Analysis of Fatty Acid Content in Different Ecotypic Rapeseed Using mrMLM

Brassica napus L. is a widely cultivated oil crop and provides important resources of edible vegetable oil, and its quality is determined by fatty acid composition and content. To explain the genetic basis and identify more minor loci for fatty acid content, the multi-locus random-SNP-effect mixed linear model (mrMLM) was used to identify genomic regions associated with fatty acid content in a genetically diverse population of 435 rapeseed accessions, including 77 winter-type, 55 spring-type, and 303 semi-winter-type accessions grown in different environments. A total of 149 quantitative trait nucleotides (QTNs) were found to be associated with fatty acid content and composition, including 34 QTNs that overlapped with the previously reported loci, and 115 novel QTNs. Of these, 35 novel QTNs, located on chromosome A01, A02, A03, A05, A06, A09, A10, and C02, respectively, were repeatedly detected across different environments. Subsequently, we annotated 95 putative candidate genes by BlastP analysis using sequences from Arabidopsis thaliana homologs of the identified regions. The candidate genes included 34 environmentally-insensitive genes (e.g., CER4, DGK2, KCS17, KCS18, MYB4, and TT16) and 61 environment-sensitive genes (e.g., FAB1, FAD6, FAD7, KCR1, KCS9, KCS12, and TT1) as well as genes invloved in the fatty acid biosynthesis. Among these, BnaA08g08280D and BnaC03g60080D differed in genomic sequence between the high- and low-oleic acid lines, and might thus be the novel alleles regulating oleic acid content. Furthermore, RT-qPCR analysis of these genes showed differential expression levels during seed development. Our results highlight the practical and scientific value of mrMLM or QTN detection and the accuracy of linking specific QTNs to fatty acid content, and suggest a useful strategy to improve the fatty acid content of B. napus seeds by molecular marker-assisted breeding.


INTRODUCTION
Rapeseed (Brassica napus L.) is one of the most important oil crops in the world, providing not only edible vegetable oil but also its potential use in lubricants and biofuels (Saeidnia and Gohari, 2012). However, the physical, chemical, and nutritional qualities of rapeseed oil depend mainly on its fatty acid composition, which consists approximately of 60% oleic acid (C18:1), 4% palmitic acid (16:0), and 2% stearic acid (18:0) (Bauer et al., 2015;Wen et al., 2015). Rapeseed oil is considered by many nutritionists to be ideal for human nutrition and superior to many other plant oils (Zhao et al., 2008;Qu et al., 2017), as it can be heated to high temperatures without smoking (Miller et al., 1987), and reduces levels of undesirable low-density lipoprotein cholesterol in the blood plasma, reducing the risk of arteriosclerosis (Chang and Huang, 1998). Optimizing the fatty acid composition is an important breeding objective for rapeseed cultivar development.
In B. napus, fatty acid metabolism is influenced by both genotype and environmental factors. Efforts to improve the oil quality have yielded many high oleic acid Brassica lines, including B. rapa (Tanhuanpää et al., 1996), B. carinata (Velasco et al., 1997), and B. napus (Pleines and Friedt, 1989;Fei et al., 2012). Further, oleic acid concentrations >70% have already been achieved in rapeseed through hybrid breeding methods (Zhang et al., 2009). Fatty acid content is a typical quantitative trait controlled by multiple genes that regulate its desaturation (Wang et al., 2015;Chen et al., 2018), and numerous quantitative trait loci (QTLs) for fatty acids have been mapped to all 19 chromosomes of B. napus, with most being found on chromosomes A01, A02, A03, A08, A10, C03, A04, A07, A09, C01, C06, and C08 (Burns et al., 2003;Zhao et al., 2008;Liu and Li, 2014;Bauer et al., 2015;Lee et al., 2015;Teh, 2015;Wen et al., 2015;Javed et al., 2016). With the increasing availability of whole-genome-sequences and SNP array development, association mapping represents a powerful approach for dissecting the genetic basis of complex quantitative traits at high resolution, which could significantly increase the precision of estimating QTL locations (Meuwissen and Goddard, 2000). Recently, genome-wide association studies (GWAS) have been performed to detect the genetic variation associated with important agronomic traits in rapeseed using the Illumina Infinium Brassica 60K SNP array (Delourme et al., 2013;Lu et al., 2014;Hatzig et al., 2015;Luo et al., 2015), including seed weight and quality , seed oil content in a panel of 521 rapeseed accessions (Liu et al., 2016), and the composition of seven fatty acids (Qu et al., 2017). Although these studies have revealed loci for associated with fatty acid traits, no beneficial alleles have been detected within the B. napus accessions.
Numerous studies showed that FATTY ACID DESATURASE 2 (FAD2) is the major gene responsible for the desaturation of oleic acid to linolenic acid (Hu et al., 2006;Peng et al., 2010;Yang et al., 2012), and four paralogs of FAD2 were previously identified in B. napus (Scheffler et al., 1997;Yang et al., 2012). These paralogs were mainly expressed in the developing seeds, suggesting possible roles in controlling oleic acid content in B. napus . In addition, KCS18, is known to play a crucial role in regulating erucic acid biosynthesis in B. napus (Wang et al., 2008;Wu et al., 2008;. However, the identified QTL were not cloned and undertaken for contributing to the minor fatty acids. Furthermore, the genetic basis of fatty acid synthesis is still unclear. The multi-locus random-SNP-effect mixed linear model (mrMLM) is emerged as a powerful tool for quantitative trait nucleotide (QTN) detection and QTN effect estimation for complex traits (Wang et al., 2016;Li et al., 2017;Chang et al., 2018;Peng et al., 2018). For example, Li et al. (2017) detected 38 significantly-associated loci and identified numerous highly-promising candidate genes (e.g., TAC1, SGR1, SGR3, and SGR5), for branch angle across 472 rapeseed accessions. Zhang et al. (2018) identified 127 significant QTNs for stalk lodging resistance-related traits using mrMLM in a population of 257 maize inbred lines. As reported by Ma et al. (2018), 127 significant QTNs with maize embryonic callus regenerative capacity were identified in a population of 144 maize inbred lines, and many candidate genes were reported to relate with auxin transport, cell fate, seed germination, or embryo development, respectively. In the present study, we analyzed the fatty acid composition in 77 winter varieties, 55 spring varieties, and 303 semi-winter varieties of rapeseed grown in three environments, and genotyped all of the accessions using the high-through Brassica 60K SNP array (Clarke et al., 2016). Then, 32,543 SNPs from the 60K SNP array were used for genome-wide association analysis usingmrMLM. In total, 149 QTNs were identified using mrMLM, suggesting that this is an effective model for identifying candidate genes underlying complex traits. Subsequently, 95 candidate genes were annotated using BlastP against A. thaliana homologs, providing insight into the genetic control of fatty acid content in B. napus. Furthermore, novel fatty acid content-associated SNPs identified here may be useful for marker-based breeding programs aimed at improving the fatty acid content of B. napus seeds.

Plant Materials
A diversity panel consisting of 55 spring, 77 winter, and 303 semiwinter rapeseed accessions (B. napus; Supplementary Table S1) was used for the association analysis. These accessions were grown in three growing seasons (2015-2016, 2016-2017, and 2017-2018) in Beibei (106.38 • E, 29.84 • N), Chongqing, China. Three rows of 10-12 plants per accession were established in the experimental fields with a randomized complete block design and three replications. Self-pollinated seeds were harvested from plants at complete physiological maturity and used for the fatty acid analysis.

Fatty Acid Measurement and Statistical Analysis
Seeds (200 mg) were homogenized with a pestle and extracted in 2 mL petroleum ether:ether (1:1) for 40 min, and methylated with 1 mL KOH/methanol (0.4 mol L −1 ). The supernatants separated by adding distilled water were identified by gas-liquid chromatography on a Model GC-2010 (Shimadzu, Kyoto, Japan). Chromatographic analysis was carried out using a fused silica capillary column DB-WAX (30 m × 0.246 mm × 0.25 um) with default parameters (Qu et al., 2017). Fatty acid profiles were calculated as a percentage of total fatty acids in the seeds, and optimized with an R script of the best linear unbiased prediction (BLUP) (Merk et al., 2012). The resulting values for each accession were used in the association analysis. All experiments were performed in triplicate, and the mean, standard deviation (SD), coefficient of variation (CV), and minimum (Min) and maximum (Max) values of the oleic acid content were calculated using SPSS 15.0 (IBM Corp, Armonk, NJ, USA).

SNP Genotyping Data Acquisition and Analysis
The methods used for SNP genotyping and mapping were described in previous reports of  and Liu et al. (2016). Using the Brassica 60K Illumina R Infinium SNP array (Clarke et al., 2016), the genotype of each accession was generated at the National Subcenter of Rapeseed Improvement in Wuhan (Huazhong Agricultural University, Wuhan, China). The low quality SNPs (call frequency <0.9 and a minor allele frequency ≤0.05) were filtered in all accessions. In addition, SNPs not accurately mapped to the B. napus genome were excluded. The probe sequences of 52,157 SNPs were used to perform a local BlastN search against the B. napus "Darmor-bzh" reference genome (version 4.1, http://www.genoscope.cns.fr/ brassicanapus/data/; Chalhoub et al., 2014) using our previously published method (Wei et al., 2015). In total, 32,543 SNPs were analyzed further.
The Q matrix of population structure was calculated by a Bayesian model-based analysis in STRUCTURE 2.1 (Pritchard et al., 2000) with published parameters of Falush et al. (2003) and Qu et al. (2017). The optimal number of K values (K = 2; Supplementary Figure S1) was determined using the Evanno method (Evanno et al., 2005). The Q matrix was selected as the fixed covariate in the subsequent association analysis (Gajardo et al., 2015). To visualize genetic relatedness among all genotypes, the principal component analysis (PCA) was constructed using the GCTA tool (Yang et al., 2011). The relative kinship matrix for each association panel was calculated using SPAGeDi (Hardy and Vekemans, 2002), and the negative values were defined as zero between two individuals, following the method of Yu et al. (2006).

Genome-Wide Association Analysis
The mrMLM significantly improved the power and precision of the GWAS, which was previously used in B. napus . Therefore, the multi-locus GWAS method (mrMLM, https:// cran.r-project.org/web/packages/mrMLM.GUI/index.html) was performed to evaluate the trait-SNP association analysis in this study (Wang et al., 2016). Moreover, the phenotypic and genotypic datasets, kinship (K), and population structure (Q) were imported into the R package mrMLM, and significantly associated SNPs were identified by mrMLM with the critical log of odds (LOD) score of 3.0 (p = 0.0002) (Wang et al., 2016). The QTNs were named using the nomenclature described by McCouch et al. (1997). For example, q-C16:0-A03-1 indicated the first locus located on chromosome A03 associated with palmitic acid.

Candidate Gene Prediction
Candidate genes were identified using significant SNP markers, which were detected using mrMLM (Wang et al., 2016). The association regions and 100-kb region upstream or downstream of peak SNPs associated with fatty acid content were identified based on the physical distance of chromosomes of significant associated-trait SNPs in the B. napus "Darmorbzh" genome (version 4.1; Chalhoub et al., 2014). Subsequently, putative candidate genes were predicted according to the annotation of the SNP-tagged genome regions and confirmed by BlastP searches against the Arabidopsis genome with an E-value ≤1E-10. The function of these candidate genes was further annotated using the Kyoto Encyclopedia of Genes and Genomes (KEGG) database (https://www.genome.jp/kegg/ pathway.html). Highly-orthologous genes involved in fatty acid biosynthesis were analyzed further, which were defined as the environment-insensitive and -sensitive genes according to the frequency detected between the different ecological genotypes and environments. To identify the directed functional genes for fatty acid, sequences of these candidate genes, isolated from plants with higher-and lower fatty acid levels, were aligned using ClustalW (Thompson et al., 1994) implemented in Geneious 4.8.5 software (Biomatters, Auckland, New Zealand).

Analysis of the Expression Profiles of Candidate Fatty Acid-Associated Genes
Total RNA was extracted from the seeds of B. napus cultivar Zhongshuang No. 11 (ZS11) at 15 developmental stages (3-49 days after pollination) using the RNAprep Pure Plant Kit (Tiangen Biotech, Beijing, China), following the manufacturer's instructions. The cDNA library construction and RNA sequencing were performed using Novogene Bioinformatics Technology (Beijing, China). Transcriptome sequencing datasets were deposited in the BioProject database (BioProject ID PRJNA358784). The data were analyzed as previously described (Qu et al., 2015), and the expression profiles of the candidate genes were quantified in terms of fragments per kilobase of exon per million mapped fragments (FPKM), using Cufflinks with default parameters (Trapnell et al., 2012). These transcriptome datasets were previously deemed suitable for selecting candidate genes (Zhou et al., 2017). Candidate genes were validated using RT-qPCR analysis. Three biological replicates and three technical replicates were performed on a CFX96 Real-Time PCR system (Bio-Rad, Laboratories, Hercules, CA, USA), and the expression levels of candidate genes were calculated using the 2 − Ct method (Zhou et al., 2017). Hence, the expression values of the 106 candidate genes were normalized by Log 2 (expression values). Heatmaps of the candidate genes were drafted using HemI 1.0 (http://hemi.biocuckoo.org/). The specific primer sequences used in this study were obtained from the qPCR Primer Database  and are listed in Supplementary Table S3.

Phenotypic Variation and Correlation Among Different Rapeseed Genotypes
Extensive variation in fatty acid content was observed between rapeseed plants of different genotypes grown in over 3 years ( Table 1). The content of palmitic acid, stearic acid and linolenic acid varied slightly among the different ecotypic rapeseed accession at different years. For example, the mean palmitic acid  content varied from 2.53 to 5.49% in spring rapeseed, 2.66 to 5.78% in winter rapeseed, and 2.69 to 6.10% in semi-winter rapeseed. The stearic acid content varied from 0.24 to 2.94% in spring rapeseed, 0.05 to 3.98% in winter rapeseed, and 0.09 to 2.89% in semi-winter rapeseed. The linolenic acid content varied from 4.90 to 14.61% in spring rapeseed, 2.15 to 11.41% in winter rapeseed, and 6.54 to 12.31% in semi-winter rapeseed. However, considerable quantitative variation was found for the content of oleic acid, linoleic acid, eicosenoic acid, and erucic acid. For instance, the mean oleic acid content ranged from 14.49 to 72.21% in spring rapeseed, 10.65 to 76.19% in winter rapeseed, and 7.91 to 83.00% in semi-winter rapeseed; the linoleic acid content ranged from 10.98 to 28.32%, 10.07 to 27.61%, and 10.41 to 28.07% in spring, winter, and semi-winter rapeseed, respectively, the eicosenoic acid content were ranged from 1.92 to 17.16%, 1.12 to 18.05%, and 0.22 to 22.34% in spring, winter and semi-winter rapeseed, respectively, and the erucic acid content ranged from 0 to 53.41 µmol g −1 , 0 to 52.83 µmol g −1 , and 0 to 48.58 µmol g −1 in spring, winter and semi-winter rapeseed, respectively (Table 1). Moreover, the largest CV (coefficient of variation) was found among the oleic acid, eicosenoic acid, and erucic acid content in different ecotypic rapeseed at different environments, ranging from 25.33 to 40.62, 50.00 to 183.00%, and 54.28 to 199.83%, respectively (Table 1), indicating that extensive variation was widely detected in the panel of accessions.
In addition, the phenotypic values of fatty acid content were displayed among the ecotypic rapeseed accessions at different years (Figures 1A-U, Table 1). Of these, the palmitic acid, stearic acid, linoleic acid, and linolenic acid content were normally distributed, but the eicosenoic, oleic, and erucic acids content were skewed for three genotypic populations in different years ( Figures 1A-U). Plants with higher oleic acid content were more common than those with lower content for each ecological type of rapeseed (Figures 1G-I). Analysis of variance (ANOVA) was performed among the spring, winter, and semi-winter rapeseed ecological types in different years, and showed that genotype and environment have significant effects on the fatty acid content of rapeseed (Table 1).
For linolenic acid (C18:3) content, a total 48 QTNs were found and covered almost the whole B. napus genome ( Table 2). Among them, seven highly identical QTNs distributed on chromosome A01, A02, A05, A06, A08, A09, and C02, along with other minor loci could be identified in different ecotypic rapeseed accessions in at least 1 year of growth ( Table 2).
In all, 149 QTNs associated with fatty acid content were detected using mrMLM ( Table 2, Supplementary Table S2), while only 62 associated regions were detected using the PCA + K model in TASSEL 5.2.1 (Qu et al., 2017). Among these, 34 QTNs were overlapped, including the association regions on A08 and C03, which had been widely reported previously (Wang et al., 2015;Liu et al., 2016;Qu et al., 2017), indicating that these results were credible and reproducible. In addition, 115 novel loci were identified for fatty acids via mrMLM compared with MLM (PCA+K; Table 2), indicating that a multi-locus random effect MLM method was better able to detect QTNs of complex quantitative traits. Furthermore, of these, 30.43% novel QTNs (35/115) were repeatedly detected for fatty acid content among different ecotypic rapeseed accessions and environments, located on chromosome A01, A02, A03, A05, A06, A09, A10, and C02, respectively (bold in Table 2), and other novel QTNs (80/115) were found for fatty acid content in a single environment. Among these QTNs, 29 were simultaneously detected in three ecotypic rapeseed, with greater QTN variation in spring and semi-winter rapeseed (Figure 2, Supplementary Table S2). Our results provideinsight into the mechanism underlying fatty acid composition, and lay the foundation for marker assisted selection in breeding projects for improved rapeseed genotypes with high oil quality and an ideal fatty acid composition.

Identification of Candidate Genes
To predict candidate genes for loci significantly associated with fatty acid content, the reported and repeadly detected novel QTNs were used to confirm the genomic regions in the B. napus "Darmor v4.1" reference genome (Chalhoub et al., 2014). We identified five environment-insensitive and fifteen environmentsensitive association regions and screened for candidate genes within these regions. Subsequently, we extracted gene sequences within the GWAS peaks in candidate association regions, and identified 95 putative genes that possibly influence fatty acid content ( Table 3). Of these, 63.16% candidate genes (60/95) were screened in the overlapping and repeatedly detected association regions, while the remaining genes (35/95) were detected on the single QTN regions ( Table 3). Using peak SNPs on A08 and C03 (Bn-A08-p13066424 and Bn-scaff_15794_3-p89999, respectively), 9 and 4 candidate genes were selected in the association regions on each chromosome, respectively (Table 3). BnaA08g11130D and BnaC03g65980D are putative paralogs of 3-ketoacyl-CoA synthase 18 (KCS18), while BnaA08g11140D and BnaC03g66040D are putative paralogs of KCS17 (Table 3), based on comparisons of the physical positions of genes associated with erucic acid traits in a GWAS (Wu et al., 2008;Lee et al., 2015;Xu et al., 2015), indicating that there is a strong correlation between GWAS peak regions and candidate genes. In addition, the putative candidate genes were characterized and annotated, such as 3-methylcrotonyl-CoA carboxylase (MCCB, BnaA08g11650D), TRANSPARENT TESTA16 (TT16, BnaA09g05410D, and BnaC02g42240D), and MYB4 (BnaC03g60080D), which might be environment-insensitive genes located on chromosome A06, A09, and C02 respectively ( Table 3). Of these, 17 genes had been identified that might be involved in the fatty acid pathway, and 12 members were annotated for fatty acid metabolism in KEGG analysis ( Table 3).
The expression of candidate genes in low-frequency association regions identified in this study were influenced     Qu et al. (2017), and QTNs with bold font were identified in at least two ecotypic rapeseed and/or environments.
Oleic acid is a monounsaturated fat beneficial for human health that contributes to the nutritional and economic value of rapeseed oil. To provide insight into the genetic control of oleic acid content, we therefore aligned 95 gene sequences from the different rapeseed accessions, and identified nucleotides in the intronic regions of BnaA08g08280D and BnaC03g60080D that show significant differences between the high-and low-oleic acid lines (Figure 3).

Expression Patterns of Candidate Genes
We assessed the relative expression levels of the candidate genes during seed development of B. napus variety ZS11 (Figures 4, 5), which had a high oleic acid content and low erucic acid. The expression levels of the environment-insensitive genes showed no obvious variation during seed development, but KCS18 (BnaA08g11130D and BnaC03g65980D), and BnaC02g42910D showed higher expression levels during the middle stages of seed development (Figure 4), indicating that they might contribute to the accumulation of oleic acid during the middle stages of seed development. In addition, TT16 (BnaA09g05410D and BnaC02g42240D) and BnaC02g42240D were expressed at high levels in the early stages of seed development, while other genes showed low expression levels throughout seed development (Figure 4).
However, we found that the expression levels of the environment-sensitive genes varied throughout seed development (Figure 5). For example, OLEO2 (BnaC04g32530D) was highly expressed in the middle and late stages of seed development, while the expression of KCS9 (BnaC07g05570D) peaked in the early and middle stages (Figure 5). In addition, BnaA01g29500D and TT4 (BnaA02g30320D) were mainly expressed in the middle stages of seed development, but KCR1 (BnaA02g13310D) and TTG1 (BnaC07g29950D) showed high expression levels throughout seed development (Figure 5). Furthermore, other genes displayed different patterns of expression throughout seed development.

DISCUSSION
In B. napus, seeds fatty acids are mainly composed of palmitic, stearic, oleic, linoleic, linolenic eicosenoic, and erucic acids, which determine the rapeseed oil quality. Enhancing the oleic acid content and quality of rapeseed through modifying its fatty acid composition has become an important breeding goal. However, previous studies identified the effect of putative fatty acid genes and the interaction of genotype and environment on the fatty acid content of rapeseed (Zhao, 2002;Zhao et al., 2005Zhao et al., , 2008Wen et al., 2015). Here, we report that fatty acid content also varies significantly among different rapeseed ecological types (spring, winter, and semi-winter rapeseed varieties) and environments (Figure 1, Table 1), indicating the complexity of the biosynthetic processes underlying fatty acid content in rapeseed. Interestingly, accessions with high oleic acid content were common amongst the different rapeseed ecological types (Figures 1G-I), possibly because these accessions are artificially selected in breeding projects aimed at producing rapeseed with high oleic acid content. Therefore, identifying the relationship between favorable alleles and environments will be beneficial for improving the fatty acid content of rapeseed.
With the development of genome sequencing and computational technologies, the Illumina Infinium Brassica 60K SNP array has been developed and widely used for the genome-wide association analysis of B. napus as well as the analysis of some trait-associated genomic regions and candidate genes Qian et al., 2014;Gajardo et al., 2015;Hatzig et al., 2015;Lee et al., 2015;Wei et al., 2015;Gacek et al., 2017;Qu et al., 2017). Furthermore, the MLM (Q+K and PCA+K) was found to be a powerful model for GWAS in these previous    The candidate genes for fatty acid content around the isolated and overlapped QTNs, respectively.
Frontiers in Plant Science | www.frontiersin.org studies (Yu et al., 2006;Zhao et al., 2011;Xu et al., 2015;Qu et al., 2017). In the present study, the mrMLM was employed for a GWAS, which is confirmed as a precise and powerful tool for analyzing phenotypic and genotypic information derived from numerous accessions and SNPs (Wang et al., 2016;Li et al., 2017). In the present study, 149 QTNs significantly associated with fatty acid content were identified using the mrMLM ( Table 2, and Supplementary Table S3).
Of these, 34 associated SNPs overlapped with those obtained using MLM (Qu et al., 2017), indicating that the association analysis was reliable; however, eight novel association regions containing 35 QTNs were simultaneously detected among different ecotypic rapeseed grown in different environments, strongly suggesting that mrMLM is more powerful to detect SNPs associated with complex traits than MLM in GWAS. In addition, 29 QTNs for fatty acids were simultaneously detected in spring, winter and semi-winter rapeseed (Figure 2), indicating that the orthologous genes for fatty acid might be better identified using these singinficant QTNs. Furthermore, more QTNs associated with fatty acids were identified from the sping and semi-winter type than in winter rapeseed (Figure 2), indicating that the fatty acids were associated with their genotype. These results might be helpful for elucidating the mechanism that determines fatty acid composition in B. napus. In B. napus, fatty acid variation is controlled by multiple genes (Zhao et al., 2008;Wen et al., 2015). In this study, we categorized the candidate genes as either environmentinsensitive or -sensitive genes, according to the published results and their detection frequency between the rapeseed genotypes grown in the different environments. A total of 95 candidate genes were identified with known functions in the fatty acid biosynthesis pathway. Of these, 34 were environment-insensitive genes ( Table 3), including KCS18, which is known to play a crucial role in regulating erucic acid biosynthesis in B. napus (Wang et al., 2008;Wu et al., 2008;, and the putative KCS18 paralogs BnaA08g11130D (Chromosome A08) and BnaC03g65980D (Chromosome C03), which are most highly expressed during the middle to late stages of seed development (Figure 4), suggesting that they are key genes regulating the accumulation of fatty acids in rapeseed oil. In addition, CER4 encodes an alcohol-forming fatty acyl-coenzyme A reductase (Rowland et al., 2006), and KCS17 is known to be involved in the biosynthesis of saturated fatty acids (Tresch et al., 2012) and would therefore be expected   Table S4). The expression values of the candidate genes were calculated using three biological replicates with three technical replicates and normalized by Log 2 (mean expression values). The "scale" function in R was used to normalize the relative expression levels (R = log 2 /mean expression values). The heatmap was generated using Heatmap Illustrator 1.0 (HemI 1.0). to be more highly expressed during seed development. We found that the putative orthologs of KCS17 (BnaA08g11140D and BnaC03g66040D) and CER4 (BnaC03g66380D) were associated with fatty acid content in B. napus, but the expression of BnaA08g11140D and BnaC03g66040D, putative KCS17 paralogs, was downregulated during seed development (Figure 4), suggesting that functional segregation exists among the different paralogs of the ancestral KCS17 gene, which has been reported previously in the B. napus genome (Chalhoub et al., 2014). In addition, the environment-insensitive genes, including DGK2 (BnaC02g42690D and BnaC02g42700D), HCS1 (BnaA09g06170D), MYB4 (BnaC03g60080D), and TT16 (BnaA09g05410D and BnaC02g42240D), might also been involved in fatty acid biosynthesis (Tasseva et al., 2004;Deng et al., 2012;Yang et al., 2012;Chen et al., 2013). Most of the environment-insensitive genes were steadily expressed throughout seed development (Figure 4), and these might be the major factors regulating oleic acid accumulation in B. napus. Furthermore, 12 environment-insensitive genes were enriched in the fatty acid biosynthesis pathway according to KEGG pathway analysis (Table 3). Importantly, the nucleotide sequences of BnaA08g08280D and BnaC03g60080D differed between the high-and low-oleic acid lines (Figure 3), but these differences were all located in the intronic regions. Furthermore, 61 environment-sensitive genes showed wide variation in expression during seed development (Table 3, Figure 5), and these could be divided into early, middle, and late expression genes, respectively. For example, BnaA06g31040D showed high expression levels in the early stages of seed development, KCS9 (BnaC07g05570D) peaked at the early and middle stages, OLEO2 (BnaC04g32530D) had high expression during the middle and late stages, and BnaA01g29500D, and TT4 (BnaA02g30320D) expression markedly increased in the middle stages (Table 3, Figure 5). Several other candidate genes, including putative homologs of FAD (3, 6, and 7), KCS (9, 12, and 21), KCR1, and LACS9, were found to be associated with fatty acid biosynthesis FIGURE 5 | Comparative expression analysis of environment-sensitive genes associated with the fatty acid content during seed development. The abbreviations above the heatmap indicate the different developmental stages of the seeds from B. napus ZS11 (defined in Supplementary Table S5). The expression values of the candidate genes were calculated using three biological replicates with three technical replicates and normalized by Log 2 (mean expression values). The "scale" function in R was used to normalize the relative expression levels (R = log 2 /mean expression values). The heatmap was generated using Heatmap Illustrator 1.0 (HemI 1.0). (Peng et al., 2010;Tresch et al., 2012;Yang et al., 2012;Lai et al., 2017;Shi et al., 2017), but their contribution remains to be confirmed by further studies (Table 3). In addition, 14 gene members were involved in the fatty acid metabolism confirmed by KEGG database analysis (Table 3). However, there is no clear evidence indicating that these genes control the fatty acid content in B. napus.
In summary, 149 QTNs for fatty acid content (including 34 reported and 115 novel loci) were detected, strongly demonstrating that mrMLM is a powerful and suitable tool for detecting QTNs for fatty acid content in rapeseed. Among these putative candidate genes, 63.16% (60/95) and 36.84% (35/95) were distributed on overlapping and isolated QTNs, respectively. Based on the pervious reports, 29 genes are involved in the fatty acid biosynthesis, and 26 gene members were enriched in the fatty acid pathway by the KEGG pathway database, indicating that mrMLM is an accurate tool to estimate the effect of QTNs on complex traits. Whether these genes exert significant regulatory effects on the fatty acid content of the seeds remains to be investigated. Hence, further studies are needed. Our findings provide useful candidate genes for the marker-assisted selection and breeding of rapeseed lines with increased oleic acid content in the seed.

AUTHOR CONTRIBUTIONS
CMQ and JL conceived and designed the experiments; MG and XH conducted the experiments; ZX, XX, LJ, and KL collected and analyzed the data; SW, and LJ made sequence alignment; YL, LW, and RW carried out the field experiments; MG, MZ, and CLQ wrote the manuscript; JL and CMQ reviewed the manuscript.

ACKNOWLEDGMENTS
We would also like to thank Kathy Farquharson for critical reading of this manuscript.