Integrate QTL Mapping and Transcription Profiles Reveal Candidate Genes Regulating Flowering Time in Brassica napus

Flowering at the proper time is an important part of acclimation to the ambient environment and season and maximizes the plant yield. To reveal the genetic architecture and molecular regulation of flowering time in oilseed rape (Brassica napus), we performed an RNA-seq analysis of the two parents after vernalization at low temperature and combined this with quantitative trait loci (QTL) mapping in an F2 population. A genetic linkage map that included 1,017 markers merged into 268 bins and covered 793.53 cM was constructed. Two QTLs associated with flowering time were detected in the F2 population. qFTA06 was the major QTL in the 7.06 Mb interval on chromosome A06 and accounted for 19.3% of the phenotypic variation. qFTC08 was located on chromosome C06 and accounted for 8.6% of the phenotypic variation. RNA-seq analysis revealed 4,626 differentially expressed genes (DEGs) between two parents during vernalization. Integration between QTL mapping and RNA-seq analysis revealed six candidate genes involved in the regulation of flowering time through the circadian clock/photoperiod, auxin and ABA hormone signal, and cold signal transduction and vernalization pathways. These results provide insights into the molecular genetic architecture of flowering time in B. napus.


INTRODUCTION
Brassica napus is a major oil crop and the largest source of high-quality vegetable oil worldwide. There are three types (spring, winter, and semi-winter), categorized by the amount of lowtemperature exposure needed for the transition from vegetative growth to flowering. Semi-winter types require little exposure to the low temperatures and are planted primarily in the Yangtze River basin of southern China. Spring rapeseed is planted in the high latitude region of China, and winter rapeseed is planted in the low latitude region. Winter rapeseed has many advantages over spring rapeseed, including grain yield, seed quality, resistance to adversity, etc., however, winter rapeseed cannot survive winter in northern China due to extreme low temperatures (<−20 • C). Several varieties (e.g., NTS57, Ganyou 4, Ganyou 221) were bred via distant cross to introduce the cold-tolerance of winter turnip rape (B. rapa) to winter oilseed rape, enabling these varieties to overwinter in northern China, which is important to alleviate the shortage of vegetable oil and prevent sand from blowing from the fallow fields. These cold-tolerant varieties require prolonged exposure to low temperatures to initiate flowering. Vernalization is critical for flowering, and therefore also to grain yield, adaptation, and survival of winter crops. Insufficient exposure to low temperatures resulting in inadequate vernalization delays the flowering of winter-type crops. It is possible to evaluate the degree of vernalization and obtain phenotypic data that can then be used to detect QTLs related to the vernalization of winter-type crops. A recombinant inbred line (RIL) wheat population was used to identify 15 QTLs linked to vernalization on nine chromosomes, and these QTLs accounted for 17 to 46% of phenotypic variation (Li et al., 2020). In total, five QTLs in four linkage groups accounted for 5.4 to 28.0% of the phenotypic variation in vernalization in an F 2 population derived from two varieties of perennial ryegrass (Jensen et al., 2005). In the rapeseed, flowering time was used to evaluate vernalization in 127 RIL population lines and identify QTLs linked to vernalization. In total, two QTLs (dtf7.1 and dtf1.1) in two linkage groups accounted for 53.1 and 8.3% of the phenotypic variation, respectively (Lee et al., 2021). One double haploid (DH) population was used to identify five consensus QTLs in the three environments, including two major QTLs located on chromosomes A03 and C08, respectively (Lee et al., 2021). Raman et al. (2012) used the flowering time to detect > 20 QTLs, including seven related to vernalization that accounted for 59.4% of the phenotypic variation. Four indexes of phenotype in the DH population in six environments were used to identify QTLs, and two vernalization QTLs in chromosomes A02 and A07 were identified (Shen et al., 2018). F 2 progeny and an F 2:3 population were used to identify 10 QTLs related to flowering time in four environments, one of which may be a novel QTL on chromosome A09 (Xu et al., 2020).
Several cold-tolerant B. napus varieties that can survive below −25.6 • C were bred in our laboratory. These strong winter types require more exposure to low temperatures for vernalization than do semi-winter and winter types, and this suggests that novel QTLs are responsible for regulating vernalization. In addition, if flowering time can be efficiently regulated using artificial measures, heterosis between two ecotypes (spring and strong winter types) can be exploited, and this is crucial for yield breeding of B. napus. Therefore, F 2 mapping of a spring-ecotype population derived from NTS57 and CY12 was used to detect QTLs related to the flowering time and vernalization for efficient breeding applications.

MATERIALS AND METHODS
An F 2 population was generated from a cross between two B. napus varieties, NTS57 (late-flowering, strong winter-type) and CY12 (early-flowering, spring-type), consisting of 174 individuals for QTL mapping. F 2 seeds were first sowed into plugs with 45 holes filled with nutrient soil in a greenhouse at the Gansu Agriculture University in Lanzhou of China. One seed per hole was sowed to foster seedlings, and 4 weeks after sowing them, they were transplanted into 25 cm diameter pots. The temperature in the phytotron was maintained in the range from 18 to 22 • C. About 9 weeks after sowing, seedlings were transplanted into a phytotron with a constant temperature of 4 • C and an 8-h light period and placed for 40 days in the phytotron. After vernalization accomplishment, seedlings were returned to the greenhouse, and the growth situation of seedlings was observed every day to accurately record the date of the first flowering. Leaf of 9 weeks' seedlings per F 2 individual was harvested. Per sample was quick-frozen in liquid nitrogen, and stored at −80 • C to extract DNA for genotype detection and construction of a genetic map.
Seeds of two parents, i.e., 17NTS57 and CY12, were sown into pots in the greenhouse with temperature ranging from 18 to 22 • C. The 9 weeks seedlings of the two parents were harvested (T0), and after harvesting, seedlings were placed into the phytotron with a temperature as low as 4 • C for vernalization. Leaves of two parents were harvested on 20 days (T1) and 40 days (T2), respectively, after initiating vernalization. Leaves harvested were quick-frozen and stored at −80 • C to extract RNA for transcriptome and qRT-PCR analysis.

Calculation of the Number of Days of First Flowering
Flowering time (FT), i.e., the number of days from sowing to the first flowering is a phenotypic date using QTL analysis, which was calculated using the following formula: FT = the date of first flowering from the date of sowing.

DNA Extraction and Targeted Sequencing
Following the protocols of Doyle and Doyle (1990), genomic DNA from the young expanded leaves was extracted from each F 2 plant. After quality assessment, each DNA sample was used to construct a library and perform targeted sequencing as described by GENOSEQ company (Wuhan, China). As sequence regions with known polymorphisms are attempted as targets for amplicon fragments, we utilized a total of 1,010 existing primer pairs from the previous SNP genotyping assays to target specific loci. These primers were mixed together in equal proportions. Mixed primers were diluted to ca. 54 nM aliquots and were used for multiplex PCR to amplify the intended target loci of each genome from F 2 progeny. Products in the amplified pools were modified by adding sequencing primer tags through secondary PCR. Modified pools were used to construct libraries for each individual. The effective concentration of libraries was determined by qRT-PCR. Libraries with more than 2 nM of the effective concentration were selected for sequencing by Pairedend 150 bp Method in the Illumina HiSeq system. After filtering unclean reads from raw reads obtained by sequencing, clean reads were obtained, which were mapped to the reference genome (ZS11-v20200127) using the MEM algorithm in BWA software. The mapping rate, coverage scope, coverage depth, and coverage uniformity of the reference genome were assessed.

Linkage Map and QTL Analysis
Loci polymorphisms of the clean reads mapped to the reference were detected using GATK software (version 3.7), and singlenucleotide polymorphisms (SNPs) and insertions and deletions (InDels) were obtained to construct a linkage map of the F 2 population according to the methods of Xie et al. (2010). QTLs were detected by a composite interval mapping method using QTL Cartographer software (version 1.17j). The LOD thresholds of the QTLs were determined by a 1,000 permutation test at a 95% confidence level.

RNA-seq Analysis
RNA was extracted from the leaves of two parents (NTS57 and CY12) after treatment at low temperature for vernalization and used in the subsequent RNA-seq. Total RNA was extracted using TRIzol Reagent (Tiange Biotech, China) according to the manufacturer's instructions. The library construction and sequencing were performed by Gene Denovo Biotechnology Co. (Guangzhou, China) on an Illumina HiSeqTM2500 platform. The raw sequences generated from the Illumina sequencing were filtered to remove the adaptor sequences and low-quality sequence reads from the data sets. Clean reads were then mapped to the B. napus reference genome using TopHat2 (version 2.0.3.12) (Kim et al., 2013). The gene expression level was calculated by the reads per kb per million reads. Genes with a fold change ≥ 2 and a p-value < 0.01 were filtered as DEGs using the edgeR package 1 . Enrichment analysis of the DEGs was accomplished based on the gene ontology (GO) database and Kyoto Encyclopedia of Genes and Genomes (KEGGs) database. Both GO terms and KEGG pathways with a Q-value of ≤ 0.05 were significantly enriched by DEGs. The clustered profiles with a p-value of ≤ 0.05 were considered to be significant profiles. The RNA-seq raw data have been submitted to the SRA of NCBI with the accession number Sra552729.

Total RNA Extraction and Real-Time Quantitative PCR
Total RNA in leaves was extracted as per the manufacturer's instructions (TIANGEN Biotech (Beijing) Co., Ltd.), and RNA integrity was detected by the electrophoresis. The RNA was reverse-transcribed (PrimeScript TM RT Reagent Kit with gDNA Eraser, TaKaRa) to obtain single-stranded cDNA. In total, three replicates were performed for each sample. After measuring the concentration, the cDNA was stored at −20 • C. The primers used for real-time quantitative PCR (qRT-PCR) were listed in Supplementary Table 1. qRT-PCR amplification reactions were performed using a LightCycler R 96 RTPCR System (Roche, Basel, Switzerland) with SYBR qPCR Mix (Invitrogen, California, United States). Actin-1 was used as an internal reference, and 1 http://www.r-project.org/  the relative expression of each gene was analyzed using the 2 − Ct method.

Phenotypic Feature of Flowering Time
For the number of days from sowing to first flowering, i.e., flowering time, the phenotypic data showed significant differences between the two parents. The flowering time of the late-flowering parent NTS57 is 121 and 20 days later than the early-flowering parent CY12 ( The broad sense heritability of FT trait among F 2 progeny is 0.621 ( Table 1). The flowering time of three individuals in the F 2 population ranged from 2 to 8 days earlier than that of the earlyflowering parent CY12, and that of two individuals was later than that of the late-flowering parent NTS57. These earlier-or laterflowering individuals are of higher value in the rapeseed breeding compared with the two parents (Figure 1).

A Targeted Sequence of F 2 Population and Variation Analysis
The mean coverage of the genome percentage and mean genome coverage depth of the parents and F 2 families were 0.08% and 108.29×, respectively. After filtrating, low-quality reads from raw reads obtained by the targeted sequencing of the two parents and 178 F 2 plants, there was an average of 308.56 thousand clean reads per plant, including 88.62 million bases. A total of 1,034 variants were obtained through comparison of the genotype of the two parents, including 946 SNPs, 88 InDels, 563 transitions (i.e., 282 A/G and 281 C/T), and 378 transversions. Among InDels, oneand two-base deletions/insertions were the most common. These variants were distributed equally throughout the genome.

Construction of the Genetic Linkage Map
A total of 1,034 markers showing polymorphism between the two parents were used to construct a genetic linkage map, 1,017 of which were grouped into 19 linkage groups. The total length of the linkage map was 793.53 cM, and the average genetic distance between the two adjacent markers was 0.051 cM. Several markers lacking signs of recombination between adjacent markers were merged into one bin, and 268 bins were contained in this linkage map with an average genetic distance of 0.233 cM between the two adjacent bins. The linkage group with the longest genetic length was C06, which was also the linkage group with the largest number of bins (28). The linkage group with the lowest number of bins (2 bins containing 19 markers) was A10, and its genetic distance (6.43 cM) was the shortest ( Table 2).

QTL Mapping of Flowering Time
A total of two QTLs (qFTA06 and qFTC08) associated with the flowering time were detected on two chromosomes (A06 and C08) in the F 2 population, with LOD values of 2.39 and 2.85, respectively. qFTA06 accounted for 19.3% of the phenotypic variation and was located on an 8.8 cM interval from 47.47 to 56.27 cM, mapping to the physical interval from 41.5 to 48.5 Mb on chromosome A06 of reference genome ZS11. qFTC08 accounted for 4.60% of the phenotypic variation and was located at the interval from 0.01 to 10.95 cM of the C08 linkage group (Table 3 and Figure 2).

Illumina Sequencing and Mapped Clean Reads
To verify transcription-level changes, six cDNA libraries were constructed using total RNA extracted from the leaves of two winter rapeseed varieties, NTS57 and CY12, treated at low temperature for 20 (t1) and 40 (t2) days and continuous room temperature (t0). The de novo transcriptome assembly and annotation were implemented. A Phred quality score of > 30 (Q30) was > 92%, and the guanine-cytosine (GC) content was consistently 48% for the four samples, indicating that the sequencing results were satisfactory. Raw sequences from cDNA libraries were filtered to obtain high-quality clean reads. After trimming and removing low-quality reads and ribosomal RNAs, 845.5 million clean reads remained and were used for quantitative analysis of gene expression. These clean reads were mapped to the reference B. napus genome using Tophat 2 software (2.

Differences in Gene Expression Between Early-and Late-Flowering Varieties
To identify the important genes related to flowering time, the difference in gene expression between two varieties, earlyflowering CY12 and late-flowering NTS57, was analyzed. NTS57 seedlings were maintained at low temperature for 20 days and contained 19,456 significant DEGs, including     (Figure 3). NTS57 seedlings were maintained at normal temperature and contained 10,820 DEGs, including 10,366 upregulated and 7,874 downregulated DEGs (| log2(TP/SP)| > 1 and P ≤ 0.001).
In total, 11,594 DEGs were common to both 20-day-and 40day-exposed B. napus, including 6,654 upregulated and 4,940 downregulated DEGs. Of these, 4,626 DEGs (2,842 upregulated and 1,784 downregulated), were specific to low-temperature vernalization, and these were most likely the basis for variation in flowering time. Therefore, these 4,642 genes were used to further narrow the region containing QTLs related to the flowering time (Figure 4).

Analysis of the Crossover Between the Candidates Identified by QTL Mapping and RNA-seq Reduced the Number of Responsible Genes
RNA-seq analysis was used to narrow down the candidate genes identified via QTL mapping and identify specific DEGs associated with vernalization at low temperature in early-and late-flowering B. napus genotypes. A total of 4,642 specific genes were identified by RNA-seq., and we speculate that some important genes affecting variations in flowering time were included in these gene groups. Two QTLs related to the flowering time, qFTA06 and qFTC08, were found on the chromosomes A06 and C08, respectively. A total of 2,926 genes were distributed in two QTL regions, including 1,194 genes in the qFTA06 region and 1,732 genes in the qFTC08 region. Crossover analysis was used  to integrate QTL mapping and global transcriptomic data. In total, eighty-four DEGs in early-and late-flowering varieties at low temperature were found in two QTL regions related to the flowering time. Of these, 46 genes were located in the qFTA06 region, and 38 genes were located in the qFTC08 region. A total of fifty-six genes was upregulated and 28 genes were downregulated in the late-flowering genotype (Figure 4). In total, six genes were candidate genes affecting flowering and were involved in flower differentiation, flowering time, and floral organ formation. The IDs of the candidate genes and their annotations are listed in Table 5. BnaC08G0115300ZS, located in the qFTC08 region and encoding a rapid alkalinization factor (RALF32), belongs to a family of cysteine-rich plant peptide hormones that are involved in the multiple physiological and developmental processes such as growth and development of shoots, buds, anther, and pollen. RALF32 binds to CrRLK1L FERONIA (FER) to inhibit the development of the shoot apex (Abarca et al., 2021). BnaC08G0356200ZS encodes plantspecific Rop nucleotide exchange factor 6, which regulates many developmental processes such as the transition from vegetative growth to the reproduction of flowering plants through the auxin and ABA-signaling pathways (Zhang and McCormick, 2007;Kim et al., 2020). BnaA06T0437200ZS encodes a transcription activator GLK2-like (also called mitogen-activated protein kinase; MAPK) that belongs to the serine/threonine kinase family and is involved in the vernalization pathway in the cold signal perception and responsive networks (Gou et al., 2018;Xu and Chong, 2018;Chuang and Tan, 2019). BnaC08G0010400ZS The 'down' or 'up' indicate the downregulation or upregulation of gene expression at the low temperature relative to that at the normal temperature.
encodes Cryptochrome 2 (CRY2), a conserved photoreceptor that accelerates flowering time by regulating the circadian clock (El-Din et al., 2003;Nefissi et al., 2011;Wu et al., 2021). BnaA06T0363000ZS, an ortholog of At3G17800 that encodes a UV-B-induced protein with an unknown function, was significantly downregulated in two varieties of rapeseed during low-temperature vernalization. BnaC08G0010400ZS encodes a hypothetical protein that is a homolog of the floweringpromoting factor (FPF1) in Arabidopsis. FPF1 confers the promotion of flowering time by affecting auxin homeostasis.

Expression of Candidate Genes in the Early-and Late-Flowering Varieties
To investigate variation in candidate gene expression, total mRNAs were extracted from leaves of early-and late-flowering varieties treated at low temperature for 20 (t1) and 40 (t2) days for vernalization and from unvernalized samples as a control (t0). The mRNA levels of the candidate genes were detected using qRT-PCR. The expression of BnaA06TG0363000ZS, BnaA06T0437200ZS, and BnaC08T0115300ZS was significantly lower in the early-flowering variety than in the late-flowering variety. Conversely, the expression of BnaA06G0332400ZS, BnaC08T0356200ZS, and BnaC08G0010400ZS was significantly higher in the early-flowering variety than in the late-flowering variety ( Figure 5). This pattern was consistent at both the RNAseq and mRNA levels, indicating that the RNA-seq results were reliable in this study.

DISCUSSION
As a sessile organism, winter rapeseed is exposed to a wide range of ambient cues and must adjust its growth and development according to this environmental variation. For overwintering or biennial crops, flowering at the proper time and in season is vital to avoid ambient damage and improve crop yield (Xu and Chong, 2018). In the northern region of China, winter rapeseed is planted in early autumn and flowers the following spring. Winter is long and frosty; therefore, strong winter rapeseed genotypes are required to avoid premature flowering (Liu et al., 2021). Simultaneously, the winter low temperatures provide sufficient cold for vernalization to initiate floral transition the following spring. In China, winter rapeseed is primarily distributed in the Yangtze River valley. No existing winter variety can survive in this cold and arid region. Fortunately, several strong winter varieties were bred by our laboratory, and these varieties have been planted successfully in cold and arid regions. Flowering time is one of the most important adaptive traits in the winter rapeseed. In previous studies, populations developed by crossing early-and late-flowering B. napus genotypes (e.g., RILs, DH, F 2 , F 2:3 posterities) were used in QTL analyses of the flowering time. Several QTLs associated with flowering time have been reported, and all but three of these (Chr. A01, A09, and C01) are distributed on 16 chromosomes (Fletcher et al., 2015;Jian et al., 2019;Scheben et al., 2020;Xu et al., 2020;Song et al., 2021 and Supplementary Figure 1). Both minor QTLs with small phenotypic contributions and major QTLs accounting for > 10% of the phenotypic variation have been reported. For example, Li et al. (2018) identified 100 QTLs associated with flowering time using a DH population containing 348 lines. These included 25 major QTLs located on chromosomes A06, A02, A03, and C06 that explained > 10% of the phenotypic variation. Two QTLs associated with flowering time were located on chromosomes A10 and C06 in an F 2 population . In this study, two QTLs associated with flowering time were found in two regions on chromosomes A06 and C08 in F 2 progeny developed by hybridizing strong winter genotype (NTS57) and spring genotype (CY12), respectively. One of these, qFTC08, was found in the 28.1 Mb interval of chromosome C08 and accounted for 8.6% of the phenotypic variation. The other, qFTA06, is a major QTL found in the 7.06 Mb interval of chromosome A06 and accounts for 19.3% of the phenotypic variation. QTL analysis identified molecular markers correlating with flowering time and allowed for the selection of oilseed rape based on the molecular markers. RNA-seq was used to detect DEGs in strong winter and spring types during vernalization. A total of 4,626 DEGs, including 2,842 upregulated and 1,784 downregulated DEGs, were identified, and their expression differed only at low temperatures during vernalization and not at normal temperatures. Several important genes for flowering time may be among these, and their differential expression may be the molecular basis of variation in flowering time between two varieties of oilseed rape. In QTL regions, qFTA06 and qFTC08, 1,138 and 1,704 genes were included, respectively. To narrow down genes that are important for flowering time, QTL mapping and RNA-seq results were integrated. In total, 2,842 genes were identified by QTL mapping and 4,626 DEGs revealed by RNA-seq were crossed and 84 overlapping genes were identified. In a subsequent analysis of KO items and the KEGG pathway, six genes were identified as candidates involved in flowering.
The six genes included in the major flowering pathway regulate the circadian clock/photoperiod, vernalization, auxin and ABA signaling, and the MAPK cascade in the cold signal transduction pathway. BnaC08G0115300ZS encodes a plant peptide hormone involved in the differentiation of the shoot apex (Abarca et al., 2021). Shoot apices can differentiate into either flower buds and flowers or branch buds and branches. BnaC08G0010400ZS encodes a homolog of FPF1 that modulates flowering time (Wang et al., 2014), and overexpression of FPF1 causes early flowering (Ge et al., 2004). In rice, FPF1like protein 4 (OsFPFL4) is a homolog of AtFPF1 that modulates auxin homeostasis to control flowering time . BnaC08G0115300ZS expression is higher in the late-flowering B. napus genotype than in the early-flowering genotype, suggesting that this gene suppresses flowering. BnaC08G0356200ZS encodes plant-specific Rop nucleotide exchange factor 6, which is involved in the regulation of the auxin and ABA-signaling pathway. The ABA-signaling pathway is involved in the regulation of photoperiod and circadian rhythms and promotes floral transition (de Montaigu et al., 2010). In our study, BnaC08G0356200ZS was upregulated and was higher in the early-flowering genotype, which implies that the upregulation of these two genes promotes oilseed rape flowering. BnaA06G0437200ZS encodes an ortholog of a transcription activator MAPK involved in the perception and transduction of the cold signal required for vernalization (Gou et al., 2018;Xu and Chong, 2018;Chuang and Tan, 2019). BnaC08G0010400ZS and BnaA06T0363000ZS regulate photoperiod and influence flowering (El-Din et al., 2003;Nefissi et al., 2011;Wu et al., 2021). In this study, these candidate genes were differentially expressed between early-and late-flowering genotypes during vernalization, which suggests that these candidate genes are important in regulating the floral translation of oilseed rape.

CONCLUSION
In this study, a genetic linkage map was constructed based on 1,017 markers that merged into 268 bins covering 793.53 cM. Two QTLs were detected that were distributed on the two chromosomes, including one major QTL on chromosome A06 that accounted for 19.3% of the phenotypic variation. RNA-seq analysis identified 4,626 DEGs between early-and late-flowering genotypes during vernalization. Integrated analysis using QTL mapping and RNA-seq identified 84 DEGs in the QTL regions. Of these, six candidate genes were involved in the regulation of flowering time and played an important role in the floweringrelated pathways.

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding author/s.

AUTHOR CONTRIBUTIONS
ZL, XD, GZ, HL, XF, YW, and HT conceived the study. GZ, CX, XC, and HL developed the experimental populations. JW and JC completed the phenotyping, genotyping, and bioinformatic analyses. XD and GZ performed the linkage mapping. GZ completed the gene expression analyses. All authors contributed to the writing and/or editing of the manuscript and approved the final version of the manuscript.