Integrating GWAS, linkage mapping and gene expression analyses reveals the genetic control of growth period traits in rapeseed (Brassica napus L.)

Background Brassica napus is one of the most important oilseed crops, and also an important biofuel plant due to its low air pollution and renewability. Growth period are important traits that affect yield and are crucial for its adaptation to different environments in B. napus. Results To elucidate the genetic basis of growth period traits, genome-wide association analysis (GWAS) and linkage mapping were employed to detect the quantitative trait loci (QTL) for days to initial flowering (DIF), days to final flowering (DFF), flowering period (FP), maturity time (MT), and whole growth period (GP). A total of 146 SNPs were identified by association mapping, and 83 QTLs were identified by linkage mapping using the RIL population. Among these QTLs, 19 were pleiotropic SNPs related to multiple traits, and six (q18DFF.A03-2, q18MT.A03-2, q17DFF.A05-1, q18FP.C04, q17DIF.C05 and q17GP.C09) were consistently detected using both mapping methods. Additionally, we performed RNA sequencing to analyze the differential expression of gene (DEG) transcripts between early- and late-flowering lines selected from the RIL population, and the DEGs were integrated with association mapping and linkage analysis to confirm their roles in the growth period. Consequently, 12 candidate genes associated with growth period traits were identified in B. napus. Among these genes, seven have polymorphic sites in the coding sequence and the upstream 2-kb sequence based on the resequencing data. The haplotype BnaSOC1.A05-Haplb and BnaLNK2.C06-Hapla showed more favorable phenotypic traits. Conclusions The candidate genes identified in this study will contribute to our genetic understanding of growth period traits and can be used as targets for target mutations or marker-assisted breeding for rapeseed adapted to different environments.


Introduction
Brassica napus (B. napus, genome AACC, 2n = 38) is one of the most important oilseed crops and an important source of protein-rich livestock feed in the world. At the same time, due to the energy crisis, low air pollution and renewability of rapeseed oil, an increasing number of people regard rapeseed oil as an ideal green energy source. This species is an allopolyploid (AACC) that evolved from an interspecies cross between Brassica rapa (B. rapa, genome AA, 2n = 20) and Brassica oleracea (B. oleracea, genome CC, 2n = 18) approximately 7500 years ago [1]. The species originated in Europe and spread worldwide (to East Asia, Australia, and North America). The growth period of rapeseed can be divided into vegetative growth period (the duration from seedling to flowering), and reproductive growth period (the duration from flowering to maturity). Flowering is a major trait in the plant's growth period, as it represents the transition from the vegetative stage to the reproductive stage. Flowering time plays an indispensable role in adaptation to specific environments. It has been reported that the growth stages show high correlations with plant height, pod number and seed yield [2]. Therefore, understanding the genetic network underlying growth period traits provides a theoretical and practical basis for developing new cultivars adapted to different geographical environments.
Linkage mapping based on bi-parental populations has been widely applied to detect QTLs for flowering time traits in B. napus [3][4][5][6][7][8][9]. Long et al. [4] identified 36 significant level (SL-QTL) and 6 micro real QTL (MR-QTL) in the doubled haploid (DH) population (Tapidor × Ningyou7) and its derived F2 population. Raman et al. [6] identified at least 20 flowering time loci localized on ten different chromosomes in the SASDH DH population (Skipton × Ag-Spectrum), which explained 2.4-28.6% of phenotypic variation. Raman et al. [7] further constructed an integrated genetic linkage map of the SASDH population and identified QTLs for flowering time accounted for up to 40.2% of genetic variation. Liu et al. [8] found 22 QTLs for days to flowering (DIF) in the spring-type DH population under multiple spring environments. Four major QTL were located on A7, C2, and C8, and explained 10.0-46.5% of phenotypic variation. Li et al. [10] used DH populations to localize QTL for flowering time in oilseed rape under winter, semi-winter, and spring ecological conditions, and identified 55 consistent QTL, including 12 environmentally stable QTL and 43 environmentally specific QTL.
Genome-wide association studies (GWAS) using natural populations provide higher mapping resolution than linkage mapping based on bi-parental segregating populations and have greater cost-effectiveness. Based on the GWAS website (https ://bigd.big.ac.cn/gwas/), approximately 8373 associations related to 19 growth period traits have been identified across six crops (cotton, maize, rapeseed, rice, sorghum, and soybean) using GWAS technology [11]. Among these associations, 198 associations for seven growth period traits were identified in rapeseed. Xu et al. [12] found 41 SNPs correlated with flowering time with a diversity panel comprising 523 B. napus cultivars, and 12 SNPs were consistent with previously identified QTLs. Niklas et al. [13] found six QTLs for beginning of flowering (BOF) and one QTL for end of flowering (EOF) in 405 B. napus inbred lines. Zhou et al. [14] identified 131 SNPs strongly linked to four earliness traits (initial flowering day, maturity time, final flowering day and flowering period), 40 of which fell into or were physically close to published flowering time SNPs. Wei et al. [15] detected 12 SNPs significantly associated with flowering time using 3,318,6 high-quality SNPs from 327 accessions. From the above results, it can be seen that these studies focused mainly on the initial flowering time. Other growth period traits, such as final flowering stage (DFF), flowering period (FP) and maturity time (MT), have rarely been reported by mapping methods.
RNA sequencing has been widely applied to examine differences in global gene expression with the benefit of high sensitivity and is more cost-effective than microarray analysis. However, transcriptome analysis often obtained a large number of differentially expressed genes (DEGs) between samples. Therefore, the integration of DEGs with GWAS or linkage mapping has been considered to be an effective way to identify candidate genes related to complex traits [16][17][18][19].
In this study, using a panel of 588 accessions, we carried out a GWAS for five growth period traits, the number of days to initial flowering (DIF), the number of days to final flowering (DFF), flowering period (FP), maturity time (MT), and growth period (GP). Furthermore, the SNPs associated with DIF, DFF, FP, MT and GP were detected in the recombinant inbred line (RIL) population. In addition, we performed transcriptomic analysis of the leaves of early-flowering and late-flowering accessions at different developmental stages (vegetative and reproductive). The main objectives of our study are: (1) to dissect the genetic mechanism of growth period traits by GWAS using a natural population and QTL mapping using an RIL population derived from GH06 and P174; (2) to dissect the molecular mechanism of floral development by characterizing the DEGs of leaves for early-flowering and late-flowering lines using RNA sequencing technology; and (3) integrating the QTLs and DEGs to identify candidate genes that control growth period traits for further verification.

Phenotypic variations, correlations, and ANOVA for growth period traits
Five growth period traits (DIF, DFF, FP, MT and GP) of 588 rapeseed lines were investigated in 3 years. Significant variation among the genotypes for five growth period traits was observed: for instance, DIF ranged from 138 to 207 days with an average of 157.98 days in 2017, from 149 to 174 days with an average of 157.34 days in 2018, and from 137 to 189 days with an average of 151.87 days in 2019 (Table 1). The coefficient variation (CV) of FP and MT was the highest (over 10%) followed by DIF and DFF, and GP was the lowest. In addition, the phenotypic frequency distributions of five traits in the 3 years all showed approximately continuous and normal distributions, which showed that the group was suitable for association analysis (Fig. 1).To test the effects of genotype (G), environment (E) and their interactions (G × E), we conducted analysis of variance (ANOVA) for each trait (Additional file 1: Table S1a). Significant variations were observed among environments and genotypes for all five traits. The DFF has the highest broad-sense heritability (h 2 ) of 92.66%, the FP has relatively low heritability  (74.11%). Overall, these findings indicated that all five traits were stably inherited.
To determine whether there are any relationships between five growth period traits in B. napus, the average value of 3 years for each trait was used for correlation analysis ( Table 2). DIF showed a significant positive correlation with DFF and GP, but was negatively correlated with FP and MT. DFF was significantly correlated with FP, MT, and GP with correlation coefficients of 0.148, − 0.664 and 0.760, respectively. It can be concluded that a shorter flowering time corresponds to a shorter growth period and a longer maturity time.

Genome-wide association analysis for growth period traits in 3 years
Through a set of processes of library construction, pairedend sequencing and SNP calling, a series of 3,856,91 highly consistent and locus-specific SNPs (minor allele frequency > 0.05 and call frequencies > 0.9) were retained for the following analysis. Additional file 2: Fig. S1 shows the density distribution of SNP markers on different chromosomes. The 3,856,91 SNPs were covered and unevenly distributed on all 19 B. napus chromosomes.
To avoid false-negative associations, three general linear models (GLM), naïve, PCA and Q models, and three mixed linear models (MLM), K, Q + K and PCA + K models were chosen to evaluate the effects of population structure (Q, PC) and relative kinship (K). According to the Q-Q plots of the six models, the MLM model can control false positives well for each trait in 3 years (Additional file 3: Fig. S2).
To minimize the effect of environmental variation, best linear unbiased predictor (BLUP) value for each line were also calculated for each trait using the R package lme4 [20]. Therefore, the MLM models were conducted for GWAS of five growth period traits with 3 years' data and BLUP values. We mainly focused on the QTLs detected in at least two environments. In total, 146 SNP loci significantly associated (− log10(P) ≥ 5.58) with five growth-related traits were identified, including 13 related to DIF, 23 to DFF, 6 to FP, 19 to MT and 85 to GP for the whole panel of accessions (Additional file 3: Fig. 2 and Additional file 4: Table S2).

Identification of candidate genes for growth period traits in B. napus
Among the 146 SNP loci significantly associated with growth period traits, 60 SNPs were divided into 19 genomic regions using a haplotype block estimation, with sizes ranging from 7 bp to 270.31 kb, while the remaining 86 SNP loci were not present in the LD blocks (Additional file 5: Table S3). Then, we obtained the genes within the same LD block or within 300 kb to either side of the significant SNPs using B. napus 'Darmor v4.1' as the reference genome. Finally, a total of 101 candidate genes were identified as orthologous to Arabidopsis flowering genes reported in the Flowering Interactive Database (FLOR-ID), most of which were involved in six flowering pathways of aging, autonomous pathway, vernalization, photoperiod, GA, and circadian clock; other genes functioned as flower development and meristem identity and flowering time integrator (Additional file 6: Table S4). The flowering time integrators FLOWERING LOCUS C (FLC) control the transition from vegetative to reproductive meristem by integrating the signals from six pathways and then precisely regulating the expression of specific flower meristem identity genes APETALA2 (AP2) and FRUITFULL (FUL) [21].
By comparing the SNP regions, 28 SNP loci (within 300 kb) associated with growth period traits, were located in or near the QTL regions identified in previous studies: 23 were reported by Raman et al., seven were reported by Wang et al. and four were reported by Zhou et al. [3,14,22,23] (Additional file 5: Table S3 and Fig. 3). Among these loci, S1_5075025, S2_5719334 and S15_5922896 were reported in two studies, and S3_13708544 and S7_10803897 were detected simultaneously in three studies. At the same time, the candidate genes we identified in this study were also reported in previous studies [14,[22][23][24][25][26][27][28]. The flowering time integrator (BnaFLC.C02, BnaFUL. A03 and BnaSVP.A09) were reported in three or more previous studies (Additional file 5: Table S3). Overall, the

QTL linkage mapping of growth period traits in 2 years
The RIL population (GH06 × P174) grown in 2017 and 2018 was evaluated for DIF, DFF, FP, MT and GP phenotypes and then used for QTL analysis. The flowering time of GH06 was later than that of P174 in the two investigated environments. A wide range of variation (Table 3) and the normal phenotypic distribution for the growth period traits (Additional file 7: Fig. S3) were observed in the RIL population, indicating the quantitative inheritance suitable for QTL mapping. ANOVA was performed on the phenotypic data from 2 years, and h 2 ranged from 49% (DFF) to 53% (GP) (Additional file 1: Table S1b). The LOD score plots across 19 linkage groups are shown in Additional file 8: Fig. S4. A total of 17, 25, 7, 21 and 13 QTLs for DIF, DFF, FP, MT and GP were detected in 2 years and were located on all B. napus chromosomes Among the 83 QTLs, 14 QTLs were associated with at least two traits. The QTL q18DIF.A01-2, q18DFF.A01 and q18MT.A01-1 have the same QTL region and explained the highest total phenotypic variation of each trait. However, the additive values of DIF and DFF were positive, and the additive value of MT was negative, which indicates that the QTL from GH06 increased DIF and DFF, while the QTL from P174 increased MT. The QTL q17DIF.A06-1, q17DFF.A06-2, q17MT.A06-3 and the QTL q17DIF.A07-2, q17DFF.A07-2, q17FP.A07, q17MT. A07-3 also has the same region and opposite additive effect. These results indicate that the same QTL may have opposite effects on the different traits. Based on the rapeseed genome annotation and the physical locations of these QTL regions, a total of 74 Arabidopsis flowering homolog genes were identified (Additional file 10: Table S6).

Identification of DEGs for floral transition and flower development by RNA-sequencing
To identify candidate DEG transcripts controlling flowering time, we sequenced four RNA samples from leaves of early-flowering cultivar 18Z134 and late-flowering cultivar 18Z88 sampled at vegetative and reproductive development stages (EV, ER, LV and LR). Every sample had two replicates. The flowering time differences between 18Z134 and 18Z88 are shown in Additional file 11: Two types of comparisons were performed: (1) to identify DEGs between the two extreme lines at each development stage: LR vs. ER and LV vs. EV; (2) to identify expression changes between different stages in each line: EV vs. ER and LV vs. LR. The false discovery rate (FDR) ≤ 0.05 and absolute value of |log2 (fold change) | ≥ 1 were used as thresholds to judge the DEGs between the two groups. In total, 3727 DEGs were identified in LR vs. ER (2421 upregulated, 1306 downregulated), 2327 DEGs were identified in LV vs. EV (1557 upregulated, 770 downregulated), 3750 DEGs were identified in EV vs. ER (2135 upregulated, 1615 downregulated), and 2038 DEGs were identified in LV vs. LR (1281 upregulated, 757 downregulated). In addition, 1662 and 683 DEGs were common to LR vs. ER/LV vs. EV and EV vs. ER/LV vs. LR, respectively (Additional file 13: Fig. S5 and Additional file 14: Table S9).

GO and KEGG analysis of differentially expressed genes
As the first criterion, we analyzed DEGs between vegetative and reproductive stages (EV vs. ER and LV vs. LR) to identify the key phase transition-associated genes. The 683 common DEGs between EV vs. ER and LV vs. LR were subjected to an enrichment analysis for GO annotation terms. The top 20 significantly enriched GO terms are shown in Additional file 15: Fig. S6a and Additional  [65] file 16: Table S10a. Among these terms, the salicylic acid biosynthetic process (GO:0009697), response to cold (GO:0009409), negative regulation of floral organ abscission (GO:0060862), intracellular auxin transport (GO:0080162) and long-day photoperiodism, flowering (GO:0048574) were involved in the phase transition process. To further determine the metabolic pathways in the phase transition process, we performed KEGG enrichment analysis (Additional file 15: Fig. S6b and Additional file 16: Table S10b). Circadian rhythm-plant (ko04712) and plant hormone signal transduction (ko04075) were significantly enriched and participate in plant development and flowering processes.
As the second criterion, we searched for DEGs between early-flowering and late-flowering lines that result in flowering time variation. The 2327 DEGs between LV vs. EV and 3727 DEGs between LR vs. ER were subjected to GO and KEGG enrichment analysis (Additional file 17: Fig. S7). DEGs in the circadian rhythm were enriched in both GO and KEGG analysis between LR vs. ER, suggesting that the different expression of genes involved in circadian rhythm could result in flowering time difference.

Identification of transcription factors and hormone-related flowering genes
In this study, 50 transcription factor (TF)-encoding DEGs, including the basic helix-loop-helix (bHLH; five members), ERF (eight members), MIKC_MADS (eight members), NAC (three members) and WRKY (four members), were identified in the early-and late-flowering lines at two developmental stages. Additional file 18: Fig.  S8a shows the overall expression trend in the four samples, and most TF families were significantly upregulated in the reproductive stages.
In addition, 20 genes in the ABA signaling pathway, 23 genes in the auxin signaling pathway, five genes in the GA signaling pathway, five genes in the cytokinin signaling pathway, four genes in the JA signaling pathway and five genes in the SA signaling pathway were identified in our transcriptome data. The genes involved in the SA signaling were upregulated in the reproductive stage. While the genes from other hormone signaling pathways have no obvious expression trend (Additional file 18: Fig. S8b).

Identification of floral transition-and flower development-related genes
According to the annotation of unigenes in Arabidopsis, a total of 125 DEGs related to flowering time were identified in the comparisons LR vs. ER, LV vs. EV, EV vs. ER and LV vs. LR (Additional file 19: Table S11). These genes mainly included the photoperiod, circadian rhythms, vernalization, GA signaling, aging, flowering time integrator and flower meristem identity genes, and the expression value of these genes in the four samples is shown in Fig. 5.
In the photoperiod pathway, 36 unigenes homologous to EARLY FLOWERING 3 (ELF3), EARLY FLOWER-ING 4 (ELF4), CONSTANS-like 9 (COL9), cycling DOF factor 1 (CDF1) and LATE ELONGATED HYPOCO-TYL (LHY) were identified. The circadian rhythm, as an internal timekeeper, controls daily and seasonal changes, is an important part of the photoperiod pathway and plays a key role in controlling plant flowering [29]. Circadian clock associated 1 (CCA1) is a key component of the Arabidopsis circadian oscillator, and it interacts with LATE ELONGATED HYPOCOTYL (LHY) and TIMING OF CAB EXPRESSION 1 (TOC1) to inhibit transcription of the Evening Complex (EC) proteins ELF4 and ELF3. In the circadian clock pathway, two CCA1 genes, three TOC1 genes, two pseudo response regulator 7 (PRR7) genes and three pseudo response regulator 9 (PRR9) genes were included.

Identification of candidates for growth period traits by integrating QTLs with DEGs
To further understand the roles of these DEGs in regulating floral transition and flower development, the DEGs were integrated with the significant QTLs identified in either association analysis or linkage mapping. Therefore, the DEGs were considered candidate genes if they were located within the confidence interval (CI) of the QTLs identified by GWAS, linkage mapping, or both. According to the above criteria, a total of 12 DEGs located in the  (Table 4 and Fig. 6). BnaAGL6. A05, which encodes a MADS-box transcription factor, negatively regulates the FLC/MAF clade genes and positively regulates FT in Arabidopsis [30]. Among them, BnaC06g15270D (BnaLNK2.C06), BnaC09g48370D (BnaGA20OX3.C09) and BnaA03g39820D (BnaFUL. A03) were also identified as EDGs between BBCH20 and BBCH50 (vernalized and nonvernalized) associated with flowering time and yield QTLs [27]. In addition, an ortholog of the circadian clock pathway, BnaTOC1. A03 and floral meristem identity gene BnaFUL.A03 were identified in both GWAS and linkage analysis [31]. Thus, we considered these five genes to be our most promising candidate genes for future prospects.
To analyze the polymorphism and explore their relationship with five growth traits, the sequence variation was analyzed in the genomic sequence including the complete coding sequence and 2000 bp upstream of the ATG translational start codon in the 12 candidate genes in the 558 lines based on our resequencing data. Among 12 genes, seven genes had polymorphic sites: 1, 15, 4, 5, 1, 3 and 10 SNPs were identified in the BnaFUL.A03, BnaELF4.A04, BnaSOC1.A05, BnaELF4.A05, BnaLNK2. C06, BnaTOC1.C09 and BnaGA20OX3.C09, respectively, and formed different types of haplotypes. Detailed haplotype information is listed in Additional file 20: Table S12.
The association of each haplotype with five growth periods was then analyzed in the GWAS population. We compared the phenotypic variations of different haplotypes for the above agronomic traits (Additional file 20: Table S12, Additional file 21: Fig. S9 and Additional file 22: Fig. S10)

Confirmation of candidate gene expression using qRT-PCR
To verify the accuracy and reproducibility of the transcriptome analysis, nine candidate DEGS listed in Table 4 were selected for qRT-PCR analysis. As shown in Fig. 7, the expression of nine genes were consistent with the RNA-Seq results in the four samples. These results demonstrated the reliability of the RNA-sequencing results.

Discussion
Growth period are important traits that can influence crop yield and quality. Identifying the genetic loci for growth period traits may help to elucidate the genetic basis underlying the growth period and be valuable for developing cultivars adapted to different geographical regions. Linkage analysis and association mapping are the two most common methods used for mapping complex traits. Linkage mapping can exclude false positives of associated loci caused by high linkage disequilibrium, and association mapping can be used to narrow the confidence interval of linkage analysis when its QTL regions are large. The combination of two mapping methods further improves mapping efficiency and accuracy [32][33][34].
In the current study, we detected the SNP loci associated with five growth period traits (DIF, DFF, FP, MT, GP) via GWAS in the natural population in three-year environments and detected QTLs that influence growth in a RIL population over 2 years in B. napus. A total of 146 SNP loci were found to be associated with growth period traits (13 related to DIF, 23 to DFF, 6 to FP, 19 to MT and 85 to GP), which were detected in at least two environments (including BLUP). Among these SNP loci, 28 were consistent with at least one QTL identified in one or more previous studies. In addition, 101 candidate genes were identified as flowering orthologous genes based on the LD block or 300 kb to either side of the significant SNPs, most of which were also reported in previous studies [14,[22][23][24][25][26][27][28]. In particular, the flowering time integrator and floral meristem identify genes (BnaFLC.C2, BnaFUL.A3, and BnaSVP.A9) were reported in at least three studies. These results indicated that the results in the present study are notably reliable. A total of 83 QTLs were observed relating to DIF, DFF, FP, MT and GP by linkage mapping in 2 years. These QTLs could explain over 10% of the phenotypic variance for DIF, DFF, FP, MT and GP. However, most QTLs were only detected in a single-year environment, which indicates that the rapeseed growth period was largely affected by environmental variation. Finally, six QTLs were associated with growth period traits using both association and linkage mapping (Additional file 23: Table S13). S19_17595451 and q17GP.C09 were associated with GP. However, the other five QTLs with overlapping confidence intervals were related to different traits in linkage analysis and association analysis. The results were also reported in Zuo et al. [35].

SNP markers with pleiotropic effects
Some GWAS and linkage mapping studies have indicated that one locus can control multiple highly correlated traits [33,36]. We detected significant correlations between the five growth period traits (Table 2), therefore, it is possible that some SNPs can affect these traits simultaneously. Five QTLs found in the GWAS and 15 QTLs identified in the linkage mapping were co-associated with at least two growth period traits, which coincided with significant phenotypic correlations among the traits (Additional file 24: Table S14a, b). It has been reported that co-mapping of QTLs for correlated traits may result from either tight linkage of multiple genes [37] or single gene pleiotropy [38]. In our research, we found that a single gene can control multiple traits (Additional file 25: Fig. S11). This finding further confirmed that co-mapping of QTLs caused by single gene pleiotropy. The comapping of QTLs can help breeders identifying favorable alleles for multiple traits simultaneously in markerassisted breeding [39].

Identification of flowering DEGs
A GO and KEGG pathway enrichment analysis indicated that the common DEGs between EV vs. ER and LV vs. LR were enriched for circadian rhythm-plant and salicylic acid biosynthetic process (Additional file 15: Fig. S6b), indicating that circadian rhythm and salicylic acid play an important role in the B. napus floral transition, which was also reported in the Annona squamosa and Rosa chinensis [40,41]. We also observed significant enrichment for genes involved in the circadian rhythm between LR vs. ER, suggesting that the flowering time difference between early-and late-flowering lines was primarily due to the expression difference of circadian rhythm genes (Additional file 17: Fig. S7b, d).
Transcription factors have been reported to play crucial roles in the reproductive development of flowering plants [42]. The MADS-box transcription factor is a major group of regulators controlling floral transition, floral organ specification and floral development in flowering plants [43].In the present study, the bHLH, ERF, MIKC_MADS, NAC and WRKY transcription factors were identified. Most IF genes were upregulated during the floral transition of B. napus (Additional file 18: Fig.  S8a), suggesting their pivotal roles in the induction of floral transition. Plant hormones, such as auxin, cytokinin (CK), abscisic acid (ABA), gibberellic acid (GA), salicylic acid (SA), and jasmonic acid (JA), are involved in the floral transition process and plant flowering [44][45][46]. The effects of GA on flowering have been extensively studied in Arabidopsis [47,48]. 62 DEGs involved in six hormone pathways were identified in our transcriptome data, which further confirmed the importance of hormones in the flowering and flower development process (Additional file 18: Fig. S8b).
In this study, we also detected 125 flowering and flower development-related homologous genes based on BLAST analysis using sequences of 306 known flowering genes in Arabidopsis. Some of these genes are related to five major flowering pathways, including circadian rhythms/photoperiod, vernalization, GA signaling and aging, while others encode regulators functioning as flowering integrators or floral meristem identification ( Fig. 5 and Additional file 19: Table S11) [49].

Candidate genes for growth period traits
The integration of QTL mapping and transcriptome sequencing is a highly useful and popular strategy for discovering candidate genes of complex traits [18]. By integrating DEGs from RNA-seq with GWAS and linkage mapping, 12 genes mainly related to circadian clock, hormone, flowering time integrator and flower meristem identity were identified as candidate genes regulating growth period traits in B. napus.
Of the 12 candidate genes detected in this study, five genes are involved in the circadian clock/photoperiod pathways: BnaA03g40160D and BnaC09g05250D are homologous to TOC1 and interact with LHY, and CCA1 forms the negative feedback loop of the circadian clock [29]; BnaA04g22640D and BnaA05g05560D are homologous to EARLY FLOWERING 4 (ELF4), which is a part of a corepressor complex consisting of ELF4, ELF3, and LUX involved in the transcriptional regulation of APRR9 and elf4 mutations result in early flowering in non-inductive photoperiods [50]; BnaC06g15270D is homologous to NIGHT LIGHT-INDUCIBLE AND CLOCK-REGU-LATED 2 (LNK2), LNK1 and LNK2 integrate early light signals with core oscillator components (TOC1, PRR) to keep track of seasonal changes in day length [51]. The up-regulation of TOC1 and ELF4 in the late-flowering lines (LR) was in accordance with the mutant phenotype in Arabidopsis; BnaC02g36310D homologous to ICE1 is involved in vernalization pathways. Cold-activated ICE1 directly induces the expression of FLC, which represses SOC1 expression, resulting in delayed flowering [52]; BnaC09g48370D is homologous to GA20OX3, which encodes a gibberellin 20-oxidase involved in the GA pathway. Both GA20ox and GA 2-oxidase genes are regulated by several transcription factors, including SHORT VEGETATIVE PHASE (SVP) and SOC1 [53][54][55]; BnaC04g15640D encodes a gene homologous to TOE1, which is a type of AP2 transcription factor and the target gene of miR172. TOE1 and its microRNA regulator miR172 play an important role in plant growth and phase change in different crops [56][57][58]; BnaA05g05000D is homologous to AGL6 (encoding a MADS-box transcription factor). AGL6 regulates flowering time through control the transcription of two key regulators of flowering time: FLOWERING LOCUS C (FLC) and FT. The agl6-1D mutant, in which AGL6 was activated by the 35S enhancer, showed early flowering under both longday (LD) and short-day (SD) conditions [59]; the floral integrator genes play important roles in activating floral meristem formation genes to induce flowering. Bna-A05g05010D was identified as a SOC1 paralog, which is a floral activator and integrates signals from the photoperiod, vernalization, and autonomous pathways in Arabidopsis [60,61]. Overexpression of SOC1 could suppress the late flowering and delayed phase transitions during the vegetative stages [62]; BnaA03g39820D was homologous to FRUITFULL (FUL) of Arabidopsis, a MADS-box transcription factor negatively regulated by APETALA1(AP1), which mediates the vegetative and meristem identity transitions by forming FUL-SVP and FUL-SOC1 heterodimers. FUL has been reported to connect several flowering pathways as a downstream target gene of age, photoperiod, and ambient temperature pathways [31]; BnaC09g23670D is homologous to AGL14, AGL14 promote flowering transition and participate in flower meristem maintenance and determinacy by positively regulating TERMINAL FLOWER 1 (TFL1) expression [63,64].

Conclusions
In this context, we conducted systematic research (including GWAS, linkage mapping and gene expression analysis) to identify candidate genes regulating growth period traits in B. napus. In summary, 146 SNPs and 83 QTLs associated with the five growth period traits (DIF, DFF, FP, MT and GP) were identified by association and linkage mapping techniques. Six QTLs were associated with growth period traits using both mapping methods. Five QTLs found in the GWAS, and 14 QTLs identified in the linkage mapping were identified as pleiotropic SNPs. RNA-Seq analysis demonstrated that the genes mediated by photoperiod/circadian rhythms, vernalization, GA signaling, aging, flowering time integrator and flower meristem identity may regulate floral transition and flower development in B. napus. An integrated analysis of QTLs and DEGs identified 12 candidate genes for growth period traits. BnaTOC1.A03 and BnaFUL. A03 were identified in both GWAS and linkage analysis; BnaLNK2.C06, BnaGA20OX3.C09 and BnaFUL. A03 have been previously reported. Among 12 genes, 7 genes have polymorphic sites in the coding sequence and 2000 bp upstream of the TSS. According to phenotype analysis, the haplotype BnaSOC1.A05-Haplb and BnaLNK2.C06-Hapla showed more favorable phenotypic traits. Therefore, our data suggest that these genes may play important roles as regulators of flowering time and growth period in rapeseed.

Plant material and phenotypic evaluation
A diversity panel consisting of 588 rapeseed inbred lines (74 winter types, 428 semi-winter types, and 86 spring types) was used for the GWAS analysis in the present study [24]. These lines were grown at the Southwest University of Beibei (29°45′N latitude, 106°22′E longitude, and an altitude of 238.57 m), Chongqing, China, in 2017, 2018 and 2019. Each variety was planted in a plot with two rows (30-cm line width and 20-cm plant distance), and each row had 10 plants. A randomized complete block design with two replications was employed.
The RIL populations of 172 RILs derived from GH06 × P174 were previously described in Liu et al. [65]. The genetic linkage map contains 9164 SNP markers covering 1832.9 cM, and 2795 SNPs were applied for QTL mapping. The population was planted in the experimental field of Southwest University of Beibei (29°45′N latitude, 106°22′E longitude, and an altitude of 238.57 m) Chongqing, China, in 2017 and 2018. The planting method is consistent with the GWAS population.
Five growth period traits, DIF (the number of days from sowing to the date when 25% plants opened first flower), DFF (the number of days from sowing to the date when 75% plants had stopped blooming), FP (the difference between DFF and DIF), MT (the number of days from DFF to the date when 75% pods were yellow) and GP (the number of days from sowing to the pod ripen) were recorded in the field trials, respectively [14]. Flowering was scored in winter, semi-winter and spring types using the same standard [15,24,28]. Analysis of variance (ANOVA) was performed using the GLM procedure of SAS.

Genome-wide association analysis
GWAS for the growth period traits were performed using Tassel 5.0 software by general linear models (GLM) and mixed linear models (MLM) methods. Population structure, relative kinship, and LD analysis were completed in previous studies [24]. A negative log (1/n) was used as a threshold for significant association SNPs with traits, where n represents the SNP number used in the GWAS [66,67]. Significant markers in the same LD block were viewed as one QTL region. The QQ plot and Manhattan plot were displayed using qqman [68] and CMplot software.

QTL analysis
Windows QTL Cartographer Version 2.5 (WinQTL-cart2.5) software was used in QTL mapping via the composite interval mapping method [69]. The logarithm of the odds (LOD) threshold for a significant QTL was calculated with 1000 permutations at a significance level of p = 0.05. A confidence interval for each QTL was defined by LOD change from the peak position. The contribution rate (R 2 ) and additive effect of a putative QTL were also calculated from WinQTLcart2.5. If the phenotypic variance explained was larger than 10%, we classified the QTL as a major QTL. QTLs were assigned names by adding the prefix of trait abbreviation, and their locations, including the chromosome number, were determined according to McCouch et al. [70]. If more than one QTL was detected on the same chromosome for a trait, QTLs were serially numbered. MapChart software was used to draw the QTL position [71].

Gene expression analysis
Early-flowering material (18Z134) and late-flowering material (18Z88) were selected from the above RIL population based on 6 years of field data. The juvenile leaves at the vegetative stage (30 days after seeding when six leaves fully expanded) and reproductive stage (20 days after flowering when 50% of the flowers are open) were collected, respectively, and immediately placed in liquid nitrogen and stored at − 80 °C. Total RNA was extracted using RNAprep Pure Tissue Kit (Tiangen, China) according to the manufacturer's instructions. RNA was pooled to yield four RNA samples, EV (leaves of early-flowering material at vegetative stage), LV (leaves of late-flowering material at vegetative stage), ER (leaves of early-flowering material at reproductive stage) and LR (leaves of late-flowering material at reproductive stage). Each sample has two biological replications. The sequencing data have been deposited in NCBI Sequence Read Archive (SRA, https ://www.ncbi.nlm. nih.gov/sra) with the accession number PRJNA540020.
DEGs between the two samples were identified based on the criteria false discovery rate (FDR) < 0.05 and |log2 (Fold change) | > 1. Gene ontology (GO) enrichment analysis of the DEGs was performed using Blast2GO, and the significantly enriched GO terms (p < 0.05) were displayed using the online tool WEGO (http://wego.genom ics.org.cn). Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis of the DEGs was performed using the KOBAS2.0 website (http://kobas .cbi.pku.edu.cn/home.do). The transcript abundance was calculated based on the fragments per kilobase of exon model per million mapped reads (FPKM). The log2transformed (FPKM + 1) values were used to study gene expression.

Candidate gene predication
According to the decay of LD, genes within 300 kb up and downstream of the significantly associated SNP loci were mined [12,72]. The QTL intervals were aligned to the B. napus reference genomes (http://www.genos cope.cns.fr/ blat-serve r/cgi-bin/colza /webBl at). Based on the physical positions of the flanking marker, genes located in the QTL confidence interval were extracted. The prediction of candidate genes referred to the two conditions: genes of known function in B. napus or genes with functionknown orthologs in Arabidopsis. Genes within the confidence interval and meeting one of the above conditions were used for gene expression analysis.

Validation of candidate genes by qRT-PCR
For qRT-PCR, the samples used for RNA-seq were used to synthesize cDNA using the QuantiTect ® Reverse Transcription Kit (Qiagen), and qRT-PCR experiments were performed with SYBR qPCR Mix (Bio-Rad) according to the manufacturer's specification. At least three biological replicates along with two technical replicates were analyzed for expression levels. BnActin7 (5′-GGA GCT GAG AGA TTC CGT TG-3′ and 5′-GAA CCA CCA CTG AGG ACG AT-3′) was used as an internal control. The 2 −ΔΔCt method was used to calculate the normalized expression of target genes [73]. The candidate gene-specific primers are shown in Additional file 26: Table S15.