Comparative Transcriptome Analysis of Gene Expression and Regulatory Characteristics Associated with Different Vernalization Periods in Brassica rapa

Brassica rapa is an important Chinese vegetable crop that is beneficial to human health. The primary factor affecting B. rapa yield is low temperature, which promotes bolting and flowering, thereby lowering its commercial value. However, quickened bolting and flowering can be used for rapid breeding. Therefore, studying the underlying molecular mechanism of vernalization in B. rapa is crucial for solving production-related problems. Here, the transcriptome of two B. rapa accessions were comprehensively analyzed during different vernalization periods. During vernalization, a total of 974,584,022 clean reads and 291.28 Gb of clean data were obtained. Compared to the reference genome of B. rapa, 44,799 known genes and 2280 new genes were identified. A self-organizing feature map analysis of 21,035 differentially expressed genes was screened in two B. rapa accessions, ‘Jin Wawa’ and ‘Xiao Baojian’. The analysis indicated that transcripts related to the plant hormone signal transduction, starch and sucrose metabolism, photoperiod and circadian clock, and vernalization pathways changed notably at different vernalization periods. Moreover, different expression patterns of TPS, UGP, CDF, VIN1, and seven hormone pathway genes were observed during vernalization between the two accessions. The transcriptome results of this study provide a new perspective on the changes that occur during B. rapa vernalization, as well as serve as an excellent reference for B. rapa breeding.


Introduction
The vegetative and reproductive growth stage are important stages in the life cycle of land plants. After a period of vegetative growth, land plants undergo reproductive growth, which is regulated by various environmental factors and endogenous signals that induce flowering. This transformation process is called floral induction, also known as flowering transition. A regulatory network ensures that plants accumulate sufficient vegetative growth over time, transforming them into reproductive growth [1]. With the development of molecular biology, five pathways that influence bolting and flowering in rice, wheat, and Arabidopsis have been identified, including the photoperiod, gibberellin (GA), vernalization, autonomous, and age pathways [2][3][4][5][6].
leaves from the center of the plants were sampled every 5 days from both treatments for a total of 50 days. Three biological replicates of the leaves were randomly collected from each individual seedling. Each treatment sample was collected separately, immediately frozen in liquid nitrogen, and stored at −80 • C for future transcriptome sequencing (i.e., RNA-Seq) and quantitative real-time polymerase chain reaction (qRT-PCR) analyses. After random sampling during each period, eight seedlings from the previous 4 • C treatment and grown in a nursery greenhouse (25 ± 2 • C under natural light). The flowering times were recorded, which began when the first petal was fully expanded.

RNA Extraction and Illumina Sequencing
Total RNA was extracted from leaf samples using an RNAprep Pure Plant kit (Tiangen, Beijing, China) following the manufacturer's instructions. RNA concentrations and quality were verified using an Agilent 2100 Bioanalyzer (Agilent Technologies, Santa Clara, CA, USA). A total amount of 1 µg RNA per sample was used as input material for the RNA sample preparations. Forty-two cDNA libraries (2 accessions ×7 periods ×3 replicates) were constructed and sequenced using NEBNext®Ultra™ RNA Library Prep Kit for Illumina (New England Biolabs, Ipswich, MA, USA) following manufacturer's recommendations. mRNA was purified from total RNA using poly-T oligo-attached magnetic beads. cDNA fragments of preferentially 240 bp in length were selected, the library fragments were purified with AMPure XP system (Beckman Coulter, Brea, CA, USA).Each library generated > 6 GB raw data. Raw reads were processed in fastq format. The adaptor sequences and low-quality sequence reads were removed from the data sets, including reads that removed more than 10% of N and more than 50% of reads with a base of Q ≤ 10. The raw data have been submitted under BioProject number PRJNA615255 to the Sequence Read Archive (SRA) database at NCBI.

Self-Organizing Feature MapAnalysis
A self-organizing feature map (SOM) is a data matrix and visualization method based on neural networks [23]. For the treatment methods of this experiment, SOM was used to design a formula containing the condition factor, time factor, and interaction between these factors. 'JWW' and 'XBJ' DEGs were processed, normalized, and used in SOM machine learning for centralizing and clustering. Through this analysis, 21,035 DEGs were identified between 'JWW' and 'XBJ'. 21305 DEGs (Table S1) was normalized by using R language _scale, which is beneficial for the next step of SOM cluster analysis. The SOM approach was used to find a set of central points and map each object in the dataset to the corresponding central point under the most similar principle [24].
In neural network terminology, each neuron corresponds to a central point. The parameters for the SOM analysis were set as follows: som_model< -som (data_train_matrix, grid = somgrid (xdim = 30, ydim = 30, topo = "hexagonal"), rlen = 100, alpha = c (0.05, 0.01), keep.data = TRUE). The SOMgrid defines the size of the network as 30 × 30 with six sides arranged between the squares. A matrix of the gene-to-sample expression values was used as input data. Genes with similar expression values were self-organized into the hexagon vicinity. The resulting map was visualized to show time-specific gene expression, in which a center point simplicity set defined multiple clusters that were grouped into the same cluster as the object closest to the center point. The expression model clustering was divided into each of the 20 'JWW' and 'XBJ' clusters.

Gene Annotations andqRT-PCR
Gene ontology (GO) enrichment analysis was conducted on the DEGs using GOseq R software [25]. KOBAS software was used to perform DEG enrichment statistics in the KEGG pathway [26]. Brassica Database v3.0 was used as reference lists [19]. Twelve DEGs were randomly selected to confirm the transcriptome data by qRT-PCR analysis. Gene-specific primers were designed using Primer v5.0 (Table S2). The normalized eEF1Bα2 gene(BraA10g020830.3C.gene) was used as an internal control of gene expression [27]. The quality-assured RNA was used as template, with the synthesis of cDNA processed by a Takara PrimeScript RT reagent Kit (Takara, Nojihigashi, Kusatsu, Japan). The Bio-Rad CFX96 RT-PCR Detection system (Bio-Rad, Hercules, CA, USA) and SYBR Green II PCR Master mix (Takara) were used for the qRT-PCR reactions. The gene expression data were analyzed using the 2 −∆∆Ct method [28]. SPSS v19.0 (SPSS, Chicago, IL, USA) was used to conduct the one-way analysis of variance (ANOVA) with Duncan's multiple range test using a significance threshold of p < 0.05. The results were visualized using SigmaPlot v10.0 (Systat Software Inc., San Jose, CA, USA).

Confirmation of Investigation Periods
From the beginning of vernalization to the end of vernalization, the number of leaves of 'JWW' is intuitively greater than that of 'XBJ', which may be related to the degree of vernalization resistance ( Figure 1A). At the same time, according to its own bolting characteristics, it can be seen that at least 25 days of vernalization can promote faster flowering of 'JWW', whereas 'XBJ' only needs 10 days of vernalization to promote flowering ( Figure 1B

DEG Analysis of 'JWW' and 'XBJ'
The two B. rapa accessions were analyzed across other periods using 0 DAT as the control for each period. Illumina RNA-Seq analysis of the 42 samples produced 1,949,168,044 total reads and 291.28 Gb clean data. The Q30 was ≥ 93.51% and the average GC content was 47.36% (Table S3A). After processing the reads, 974,584,022 clean reads were obtained and sequence aligned with the designated B. rapa reference genome with an alignment efficiency ranging from 84.62% to 93.33%. The unique map ratio ranged from 77.98% to 90.88%, and the multiple map ratio ranged from 2.07% to 6.65% (Table S3B). Comparison with the B. rapa reference genome, 44,799 known genes and 2,280 new genes were identified.
Multiple groups were differentially analyzed using 0 DAT (J1 and X1) as the controls to identify vernalization-related DEGs (  To determine the functional distribution of unigenes, GO annotation was used to classify the DEGs generated by 'JWW' and 'XBJ' during vernalization ( Figure 3). The two B. rapa accessions exhibited a high degree of similarity in terms of unigene distribution. In the biological process category, unigenes were concentrated in the "metabolic process", "cellular process", and "single-organism process". In the cell component category, most unigenes were related to "cell", "organelle", and "cell part". In the molecular function category, unigenes were highly related to "catalytic activity" and "binding". These highly enriched GO terms indicated the basic regulatory and metabolic functions of the two B. rapa accessions during vernalization. For the metabolic processes in the biological category, we selected the vernalization periods common to 'JWW' and 'XBJ' for comparison, which are 0 DAT (J1 and X1), 25 DAT (J2 and X4), and 50 DAT (J6 and X6). Compared to 0 DAT, we found that the change in 25 DAT and 50 DAT expression was larger in 'XBJ' than in 'JWW'. Similarly, the catalytic activity-related genes in molecular functions are more pronounced in 'XBJ'. J2 and J6 are enriched in this term by 2723 and 3303 DEGs, respectively, compared with J1.X4 and X6 are enriched in this term with 3461 and 3630 DEGs, respectively, compared with X1. Overall, in these three GO analysis categories, the expression of 'XBJ' in many terms has changed more than 'JWW'.

Analyses of Dynamic Gene Expression Patterns during Vernalization
To determine the main transcriptional dynamics related to vernalization of the two B. rapa accessions, SOM was used to identify the dynamic changes of the 21,035 DEGs in different vernalization periods, as well as a cluster analysis, which was performed on multiple groups of expression patterns. The expression of normalized DEGs in each period is represented by different colors of each hexagon in the square; the red to blue gradient indicated upregulation and downregulation, respectively ( Figure 4A,C). Each DEG that exhibited a similar expression level was clustered into agrid and grids with the same color from each cluster were grouped together ( Figure 4B,D). The characteristic values of the gene expression levels of each grid were extracted; 20 clusters were defined according to the number of grids (Figures S1 and S2). Then, GO classification analysis was performed on each cluster to analyze the regulated DEGs of the two B. rapa accessions during vernalization in the three GO categories (Tables S4 and S5). Clusters that were notably regulated during vernalization were selected to analyze the differential expression patterns of the two B. rapa accessions. In clusters 7, 10, 11, 13, and 15 of 'JWW', comparing 0 DAT and 35 DAT 25 • C, DEGs exhibited a state of continuous induction during vernalization. These genes were upregulated starting at 25 DAT and continued through 50 DAT. These continuously upregulated genes possessed biologically important functions and were used for studying the sensitivity of 'JWW' to vernalization. GO annotations were used for classifying the continuously increasing DEGs. Cellular components were enriched in genes related to extracellular space, lysosome, nucleolus, and endosome. In the biological processes category, DEGs were mainly enriched in regulation of gene expression and epigenetic. In the molecular functions category, these clusters had different DEGs that were enriched. Specifically, receptor activity-related genes were enriched in clusters 10, 13, and 15. Clusters 7 and 11 were enriched in genes related to enzyme regulator activity. The genes related to binding, including oxygen and lipid binding, were enriched in each of the five clusters. In clusters 3, 12, and 19 of 'JWW', the DEGs remained stable and exhibited low expression levels at 25 DAT. Cluster 12 reached its lowest level at 50 DAT. The GO analysis revealed that genes related to extracellular space, nucleolus, and endosome in the cellular components category were enriched in these three clusters. Cluster 3 was enriched in genes involved in the regulation of gene expression and epigenetic in the biological processes category. The molecular functions category was enriched in activity, including receptor, nuclease, and enzyme regulator activities, and binding, including oxygen and lipid binding, genes in these three clusters.
In 'XBJ' at 10 DAT, clusters 7 and 10 were compared between 0 DAT and 25 DAT 25 • C, revealing that the genes were continuously upregulated. The GO analysis revealed that these genes were enriched in peroxisome and nucleoplasm in the cellular components category, cellular homeostasis in the biological processes category, and lipid binding in the molecular functions category; related genes were all enriched in these two clusters. Additionally, signal transducer activity and carbohydrate binding-related genes in the molecular functions category were enriched in cluster 10 of 'XBJ'. During the vernalization of 'XBJ', clusters 1, 2, and 8 exhibited stable or low gene expression and were enriched in nucleoplasm of the cellular components category. During the vernalization of 'XBJ', clusters 1, 2, and 8 with stable or low gene expression were enriched in nucleoplasm of the cellular components category. Meanwhile, related genes of cluster 2 were enriched in golgi apparatus, and related genes of cluster 8 were enriched in nuclear envelope. Among the biological processes category, only cluster 8, which contained low-expressing genes, was abundant in cellular homeostasis. Moreover, low-expressed genes were enriched only in cellular homeostasis. The molecular functions category was enriched in low-expressing genes related to signal transducer activity and lipid binding in clusters 1, 2, and 8.
To further understand these continuously upregulated and downregulated genes in the enrichment pathways, these clusters were analyzed by a KEGG pathway analysis (Tables S6 and S7). In 'JWW', cluster 7 possessed continuously upregulated genes that were notably enriched in glycolysis/gluconeogenesis, starch and sucrose metabolism, and amino sugar and nucleotide sugar metabolism. In cluster 15, the DEGs of three aspects, including proteasome, plant hormone signal transduction, and glycolysis/gluconeogenesis, were notably expressed. Notably, during vernalization, these continuously upregulated clusters were considerably enriched in the plant hormone signal transduction and starch and sucrose metabolism pathways. Moreover, the plant hormone signal transduction and starch and sucrose metabolism pathways in clusters 3, 12, and 19 were continuously considerably downregulated. Analysis of the 'XBJ"s gene clusters that were continuously upregulated and downregulated during vernalization revealed that the plant hormone signal transduction and starch and sucrose metabolism pathways were also considerable. Therefore, a detailed analysis was performed on the related gene expression profiles of the plant hormone signal transduction and starch and sucrose metabolism pathways. Additionally, in order to fully explain the effect of vernalization on the two B. rapa accessions, the photoperiod, circadian clock, and vernalization pathways were incorporated in the analysis.

DEGs of the Plant Hormone-Related Signal Transduction Pathway
The plant hormone signal transduction pathways exhibited different DGEs levels during the two B. rapa accessions vernalization. To identify the hormone similarities and differences of the two different bolting B. rapa accessions, the expression levels of the hormone-related DEGs were further analyzed ( Figure 5).

DEGs of the Photoperiod and Circadian Clock Pathways
The photoperiod and circadian clock pathway-related genes are often referred to as flowering-related genes. The two pathways of the red and blue light absorbing plants were compared and analyzed to the 23 flowering-related genes identified in this study (Figure 7; Table S10). In 'JWW' from J1 to J4, 12 genes were considerably upregulated (log2(J4 FPKM /J1 FPKM ) > 1). In 'XBJ' from X1 to X4, only 7genes were considerably upregulated (log2(X4 FPKM /X1 FPKM ) > 1) and 6were considerably downregulated (log2(X4 FPKM /X1 FPKM ) < −1). Phytochrome B(PHYB), two-component response regulator (APRR), and protein LHY (LHY) genes in the red-light pathway that positively regulate flowering locus T (FT) between the B. rapa accessions was slightly different during vernalization. It was worth noting that the two phytochrome B (PHYB) genes upregulated in J4 versus J1 were downregulated in J4 versus JCK. Ten APRRs were up-and downregulated during vernalization; the expression levels of APRRs in 'JWW' were generally higher than 'XBJ'. The expression of LHYs, which play intermediate roles in the red-light pathway, were not notably different between the two accessions. The cryptochrome-2 (CRY), gigantea (GI), and cyclic dof factor (CDF) genes in the blue light pathway were negatively regulated. CDFs were also regulated by red light-related genes. In the comparison between J4 and JCK, the GI and CDF genes in 'JWW' were upregulated during vernalization. Comparing X4 to XCK, vernalization slightly upregulated the GIs of 'XBJ', but downregulated CDFs. 'XBJ' also upregulated flowering genes, including zinc finger protein constants (CO) and FTs, during vernalization at higher levels than 'JWW'. These results indicated that the photoperiod and circadian clock pathways may play important roles in regulating B. rapa flowering during vernalization.

Identification of VernalizationPathway DEGs
Studies on the molecular mechanism of vernalization in dicotyledons have focused on the expression changes of FLC and its related regulatory genes. In this study, FLCs and their related vernalization regulating genes were identified (Figure 8; Table S11). VIN3, VRN2, and VRN1 genes were found to be negative regulators of FLC. Moreover, vernalization promoted the upregulation of VIN3s. The expression patterns of VRN1 between the two B. rapa accessions exhibited contrasting trends. Moreover, VRN2 did not change much during vernalization. Vernalization reduced all four FLCs and the reductions in 'XBJ' were slightly greater than 'JWW'. FRI positively regulates FLC [10], but was not affected by vernalization. It can be observed that FRI was not very different between X4 and XCK, and there was not much difference in FRI expression between J4 and JCK. Although FRI expression in both accessions continuously increased as vernalization progressed, it could not be explained by the effects of vernalization. Downregulated FLCs negatively affected the upregulation of FTs, which promoted the upregulation of the flowering factor, suppressor of constans overexpression 1 (SOC1). SOC1 in 'XBJ' was upregulated more than 'JWW'.

Discussion
Studies on transcriptome regulation based on phenotypes is a powerful research method [29,30]. Using this analytical method, DEGs and related regulatory pathways produced by Asiatic lily and tropical lotus during vernalization have been investigated [31,32]. With regard to B. rapa, the reference genome has been updated for several generations, and various regulatory pathways have been continuously studied [33,34]. However, only a few studies have been conducted on the vernalization control networks in B. rapa at different periods. Thus, in order to clarify the associated regulatory mechanisms, two bolting accessions were selected for transcriptome analysis in this study. A total of 974,584,022 clean reads and 291.28 Gb clean data were obtained. After comparison with the B. rapa reference genome, 44,799 known genes and 2280 new genes were identified. A SOM analysis was performed and screened 21,036 DEGs in 'JWW' and 'XBJ' (Figure 4). GO and KEGG analyses were also performed on clusters that were continuously up-and downregulated during vernalization. DEGs that were enriched in the hormone signal transduction and sucrose metabolism pathways of these continuously up-and downregulated clusters were identified. Changes in the DEGs of these pathways during vernalization in the two B. rapa accessions were the primary focus of this study. Moreover, DEGs in the photoperiod and circadian clock, and vernalization pathways were analyzed.

Changes in the Plant Hormone Signal Transduction Pathway during Vernalization
The transition from vegetative to reproductive growth is a necessary process in plant flowering and reproduction. The shoot tip meristem (SAM) is an important part of the plant that helps complete the transition from vegetative to reproductive growth. During the vegetative growth stage, leaves formed around the SAM will generate lateral meristems. Then, a series of environmental influences and induced changes in endogenous signals will change the state of the SAM and induce flowering [35]. Plant hormones play a role in the flowering time syndrome of Arabidopsis [36]. However, in B. rapa, the regulation of hormones during vernalization remains unclear. Here, a layered heat map of the 7 plant hormone pathways of the two B. rapa accessions during vernalization was constructed to analyze this process in detail ( Figure 5, Table S8).
Auxin functions in the initiation of floral primordia and flower development [37]. Aux/IAAs are primarily responsive to auxin genes. Aux/IAA stability and activity are regulated by auxin [38]. The interaction between Aux/IAA and TIR1 plays a role in regulating auxin signal transduction [39]. Mutations that occur in Aux/IAA and TIR1affect the interaction, resulting in low-auxin responses [40]. Similar to carbohydrates, auxin concentrations affect flower induction in different ways. Specifically, low concentrations promote flowering, while high concentrations delay flowering [36]. In Arabidopsis, the local accumulation of auxin releases ARF5 due to Aux/IAA inhibition, which subsequently promotes a rise in CH3 levels, as well as high expression levels of SAUR family genes, and is the leading cause of flowering on the lateral side of the SAM [41]. In this study, 'JWW' and 'XBJ' upregulated ARFs, CH3s, and SAUR family genes with vernalization. These genes may be upregulated in the same way that auxin promotes flowering in Arabidopsis. This result indicates that, under the influence of vernalization, the auxin genes of the two B. rapa accessions may affect flower meristems, as well as their bolting and flowering abilities at different times.
GAs are a plant hormone involved in the floral transformation process of Arabidopsis and play a key role in controlling flowering time [2]. In Arabidopsis, GID1 is a positive regulator of GA signaling. The GID1−GA complex accelerates the degradation of the DELLA protein [42]. In this study, two GID1s were upregulated in the first four vernalization periods, but notably declined in the fifth period of 'JWW'. However, in 'XBJ', GID1 levels varied depending on the vernalization stage. The three DELLAs exhibited the same trend in both accessions. As vernalization continued, two DELLAs were downregulated and one was upregulated. These results indicated that the GA signaling pathway may play a role in the transition from the growth period during vernalization of these two accessions.
ET is a gas phytohormone that participates in plant stress responses and regulates a variety of developmental processes [43]. In Arabidopsis, the ET receptor possesses fivemembers and has been classified based on certain characteristics into ET receptor 1 (ETR1) and ET receptor 2 (ETR2) [44]. The serine/threonine-protein kinase CTR1 gene plays an important role in ET signal transduction [45]. There is a MPK module located between CTR1 and ethylene-insensitive protein 2 (EIN2) [46]. A previous study demonstrated that once MPK6 loses its function, the ET signaling pathway will return to normal, indicating that MPK6 may not participate in this pathway [47]. However, in the two B. rapa accessions of this study, 'XBJ' exhibited a large change in two MPK6s during vernalization compared to 'JWW'. Moreover, EIN2 is the positive regulator of the ET signaling pathway [48]. EIN3 is located downstream of EIN2. EIN2 and EIN3 work together and act on the ET-responsive transcription factor (ERF1/2) [49]. During vernalization, ERF1/2 genes were differentially expressed in the two B. rapa accessions of this study, which may explain the different conversion flowering periods.
ABA is a plant hormone that regulates plant development [50]. PYL family proteins are receptors of the ABA signaling pathway and inhibit PP2C in an ABA-dependent manner [51,52]. PP2Cs interact with SnRKs and inactivate them due to phosphorylation [53]. A previous study demonstrated that ABA played an inhibitory role in flowering time by modulating DELLA activities in the GA pathway [54]. A separate study found that ABA may play a role in hindering flowering and affect flowering by altering FLCs [55]. In this study, PYL genes were differentially expressed in both accessions and exhibited different trends. In the fourth period of vernalization, PYLs in 'JWW' were considerably downregulated compared to the control, while the PYLs in 'XBJ' were upregulated and downregulated compared to the control. Moreover, vernalization continuously upregulated PP2C genes, thereby blocking the generation of some SnRK2 genes; the difference was more noticeable in 'XBJ'. Therefore, the ABA signaling pathway may participate in regulating different growth period transformations in these two B. rapa accessions.
JA is a lipid plant hormone that plays multiple roles in plant growth and development [56]. CKs are a class of classical hormones that regulate both the division cycle and meristem homeostasis [57]. SAs are involved in plant defense functions and regulate leaf senescence [58]. In this study, the genes involved with these three plant hormones changed during vernalization. jasmonic acid-amido synthetase JAR1 (JAR1) is one of 19 closely related Arabidopsis genes that issimilar to the auxin-induced soybean GH3 gene [59]. Methyl JA (MeJA)-treated vernalization in insensitive spring wheat exhibited flowering and considerably downregulated TaVRN1 and TaFT1 genes, suggesting that MeJA may modulate vernalization and flowering time in wheat [60]. Previous studies have reported that CK promotes the formation of floral meristems and exogenously applied CK promotes Arabidopsis flowering [61,62]. The phytochrome signaling intersects with SA signaling [63], which correlates with flowering [64].
To summarize, different hormone biosynthesis and signaling pathways were differentially expressed during vernalization in the B. rapa accessions of this study. Clearly, this complex hormone signaling network accelerates the transition from vegetative to reproductive growth during vernalization in B. rapa.

Gene Expression Changes in the Starch and Sucrose Metabolism Pathways
The rapid propagation of plants is closely related to flowering time; carbohydrates are thought to play a crucial role in this process [65]. Sucrose and starch are the primary products of photosynthetic carbon assimilation. In lilies, sucrose treatment accelerates flowering [66]. Moreover, SUS is a key enzyme involved in the sucrose metabolism of plants [67], while SPS is a key enzyme that catalyzes the first step of the sucrose synthesis pathway [68]. 'XBJ' possessed a new SUS gene (newGene_3912), and the expression of SPSs in this accession were slightly higher than 'JWW' (Figure 6, Table S9).
Starch is the main carbohydrate stored, which allows them to meet their long-term carbon needs [69]. The main component of starch is synthesized by starch synthases (SS2), GBEs, and PYGs [70,71]. 'XBJ' highly expressed three PYGs and GBEs in each vernalization period. These differences indicated that the starch and sucrose metabolism pathways may play important roles in regulating B. rapa with different bolting times, but different roles in activating the SAM. Maltose is produced during the conversion of starch to sucrose in leaves under dark conditions [72]. BAM could potentially produce maltose as a key metabolite and osmolyte via starch hydrolysis [73]. Maltose is the primary source of carbon during cellular metabolism and growth at night [74]. Four BAMs were upregulated during vernalization in the two B. rapa accessions. However, when comparing the same vernalization period to the control, five BAMs were upregulated in 'JWW' and three in 'XBJ'. These results indicated that maltose may provide the two B. rapa accession with sufficient energy required for normal growth during vernalization at night.
In Arabidopsis, TPS plays a critical role in controlling the transition from vegetative to reproductive growth [75,76]. TPS catalyzes the formation of trehalose-6-phosphate (T6P) from glucose-6-phosphate and uridine diphosphate-glucose (UDP-glu) and is a signaling molecule that relays information regarding carbohydrate availability to other signaling pathways [77,78]. The loss of TPS resulted in extremely delayed flowering in Arabidopsis. The TPS pathway also acts as a signal that coordinates the induction of flowering by regulating the expression of key floral integrators in the leaves and SAM [79]. In this study, the TPSs of 'XBJ' expressed during vernalization were less than 'JWW'. It may be that the late-bolting 'JWW' accession requires higher TPS expression in order to promote faster flowering. The upregulation of pgms promotes the production of polysaccharides, which has been demonstrated in mushrooms [80]. In this study, pgms increased at first and then decreased during vernalization. Therefore, it is speculated that pgms rose rapidly at the beginning of vernalization in order to provide sufficient energy for normal growth, and then it gradually decreased after the vernalization transformation period. In young and mature leaves, UGP is primarily involved in the sucrose biosynthesis pathway, the major form of transport of carbon in plants, which provides energy to plants [81]. In this study, two UGP genes exhibited opposing states during vernalization, in which, one was upregulated and the other was downregulated. The underlying mechanisms for these contrasting trends, however, require further investigation.

Photoperiod and Circadian Clock Pathways
The photoperiod and circadian clock pathways function as a whole, as plants coordinate their metabolism according to changes in photoperiods, which thereby induces flowering. Expression changes in the CO/FT mode is the core link between the photoperiod induction pathway [3]. The purpose of this study was to analyze the influence of vernalization on the related genes in this pathway under certain photoperiods (Figure 7, Table S10). CO is a typical circadian clock control gene that acts on the FT promoter region and is the main positive regulator of FT [82]. The red-light pathway positively regulates CO and FT, which were continuously regulated by CDFs through PHYBs, LHYs, and APRRs. Among them, APRRs were considerably different between the two B. rapa accessions and may affect bolting. The two PHYBs upregulated in J4 versus J1 were downregulated in J4 versus JCK, suggesting that PHYBs may not be affected by vernalization in 'JWW'. CDFs are intermediate regulatory genes in this pathway that are regulated by both red and blue light. Specifically, vernalization downregulated the CDFs of 'XBJ', but upregulated the CDFs of 'JWW'. Thus, CDF may be a gene that does not affect late-bolting B. rapa during vernalization.

DEGs of the Vernalization Pathway
Land plant flower development undergoes floral induction, floral primordia, and floral organ development phases. In the flowering induction stage, although there are no morphological changes in the SAM, a series of related genes centered around FT undergo various changes, forming a complex regulatory network to achieve the transition from vegetative to reproductive growth [83]. In B. rapa, vernalization is the most important process that induces flowering. Currently, studies on the underlying molecular mechanism of the vernalization pathway in A. thaliana mainly focus on changes in the expression levels of FLCs. Among them, BrFLC1 and BrFLC2 of B. rapa respond to the regulation of vernalization and affect flowering time [84]. In this study, three of the four FLCs decreased considerably during vernalization. The decrease in 'XBJ' was larger than that in 'JWW', and it continually decreased in the first five vernalization periods ( Figure 8, Table S11). VRN1, VRN2, and VIN3 are negative regulators of FLC, and previous studies demonstrated that VRN1 and VRN2 in Arabidopsis were not induced by vernalization and only play a role in the persistently low expression of FLC at the end of vernalization [12][13][14]; thus, the results of this study were inconsistent, as two B. rapa VRN1s exhibited opposing modes during vernalization. However, these differences require further investigation. Changes in VRN2 were consistent with previous results [13]. VIN3s were upregulated in the first five vernalization periods and downregulated in the last vernalization period. Therefore, it is speculated that B. rapa transitioned to reproductive growth just before the last vernalization period and, hence, the downregulation of VIN3s.

Conclusions
To date, studies on the mechanisms related to vernalization of B. rapa have focused only on the comparison between the vernalized and unvernalized, and related research approaches have been relatively narrow. Vernalization promotes plant bolting and flowering as an important trait in productivity and reproduction; thus, it is necessary to clarify its effects on bolting and flowering in different vernalization periods. In this study, two B. rapa accessions with different bolting times were investigated over a series of vernalization periods. For the related flowering traits, the vernalization periods were selected for comparative transcriptome analysis to describe the effects of different vernalization periods on B. rapa. According to the results, the B. rapa accessions exhibited specific differences at each vernalization period, indicating that research on vernalization should be a process of vernalization, not a comparison between vernalized and unvernalized. In comparing the two accessions, key genes of seven hormone pathways of the plant hormone signal transduction pathway during vernalization were notably different, and the difference between late-bolting 'XBJ' was greater than early-bolting 'JWW'. This difference was also confirmed in the vernalization pathway. In the starch and sucrose metabolism and photoperiod and circadian clock pathways, TPS, UGP, and CDF exhibited different expression patterns, but these differences require further analysis. Notably, VIN1, which does not play a role during vernalization in other dicotyledons, exhibited opposing expression patterns in this study, which may be a new avenue for research on vernalization in B. rapa in the future. The transcriptional regulation of B. rapa vernalization was summarized in this study (Figure 10), which provides a foundation for future investigations on the molecular mechanisms of B. rapa with different abstractions during vernalization. Supplementary Materials: The following are available online at http://www.mdpi.com/2073-4425/11/4/392/s1, Figure S1. Clusters expressions associated with 'JWW' vernalization periods. The x-axis represents vernalization periods, and the y-axis represents expression. Figure S2. Clusters expressions associated with 'XBJ' vernalization periods. The x-axis represents vernalization periods, and the y-axis represents expression. Figure S3. Linear relationships between qRT-PCR data and RNA-Seq data of related genes. Table S1. The 21035DEGs in both B. rapa accessions at different periods (J for 'JWW'; X for 'XBJ'). Table S2. Specific primers for differential gene sequences for qRT-PCR. Table S3. (A) Statistics of Illumina sequencing data (J for 'JWW'; X for 'XBJ'); (B) The sequencing and mapped results obtained from each sample (J for 'JWW'; X for 'XBJ'). Table S4. GO classification of DEGs in each clusters for the 'JWW'. Red indicates the number of terms genes that continuously increase; blue indicates the number of terms genes that continuously decrease. Table S5. GO classification of DEGs in each clusters for the 'XBJ'. Red indicates the number of terms genes that continuously increase; blue indicates the number of terms genes that continuously decrease. Table S6. Analysis of pathways involving up-and downregulated genes for clusters in 'JWW'. Table S7. Anlysis of pathways involving up-and downregulated genes for clusters in 'XBJ'. Table  S8. Plant hormone related DEGs in two B. rapa accessions during the vernalization process. Table S9. Starch and sucrose related DEGs in two B. rapa accessions during the vernalization process. Table S10. Photoperiod and circadian clock related DEGs in two B. rapa accessions during the vernalization process. Table S11. Vernalization related DEGs in two B. rapa accessions during the vernalization process.