Regulatory Networks of Flowering Genes in Angelica sinensis during Vernalization

Angelica sinensis is a low-temperature and long-day perennial herb that has been widely used for cardio-cerebrovascular diseases in recent years. In commercial cultivation, up to 40% of flowering decreases the officinal yield of roots and accumulation of bioactive compounds. Although the regulatory mechanism of flowering genes during the photoperiod has been revealed, the networks during vernalization have not been mapped. Here, transcriptomics profiles of A. sinensis with uncompleted (T1), completed (T2) and avoided vernalization (T3) were performed using RNA-seq, and genes expression was validated with qRT-PCR. A total of 61,241 isoforms were annotated on KEGG, KOG, Nr and Swiss-Prot databases; 4212 and 5301 differentially expressed genes (DEGs) were observed; and 151 and 155 genes involved in flowering were dug out at T2 vs. T1 and T3 vs. T1, respectively. According to functional annotation, 104 co-expressed genes were classified into six categories: FLC expression (22; e.g., VILs, FCA and FLK), sucrose metabolism (12; e.g., TPSs, SUS3 and SPSs), hormone response (18; e.g., GID1B, RAP2s and IAAs), circadian clock (2; i.e., ELF3 and COR27), downstream floral integrators and meristem identity (15; e.g., SOC1, AGL65 and SPLs) and cold response (35; e.g., PYLs, ERFs and CORs). The expression levels of candidate genes were almost consistent with FPKM values and changes in sugar and hormone contents. Based on their functions, four pathways that regulate flowering during vernalization were mapped, including the vernalization pathway, the autonomic pathway, the age pathway and the GA (hormone) pathway. This transcriptomic analysis provides new insights into the gene-regulatory networks of flowering in A. sinensis.


Introduction
Angelica sinensis (Oliv.) Diels (Apiaceae, alt. Umbelliferae) is a perennial herbaceous species, originally native to China, with a dense population in Gansu province at an altitude of 2200-3000 m [1]. The roots have been widely used as a traditional Chinese herbal remedy in China, Korea and Japan for nourishing the blood, regulating menstrual disorders, relaxing bowels, etc., for over 2000 years [2,3]. Currently, over 50 compounds have been isolated from the roots, including phthalides, organic acids and polysaccharides, which confer pharmacological effects, including cardio-cerebrovascular, hepatoprotective, antioxidant, antispasmodic and immunomodulatory effects [4][5][6].
At present, over 43,500 ha of commercial of A. sinensis cultivation supplies the increasing market demand; however, up to 40% of flowering of plants in the field makes the roots lignified, decreasing the yield of roots and the accumulation bioactive compounds [7][8][9]. Previous studies have demonstrated that A. sinensis is a low-temperature and long-day plant that must experience vernalization (0 to 5 • C) and long-day conditions (>12 h daylight) for the transition from vegetative growth to flowering [10]. Thereafter, the flowering could be effectively avoided after the seedlings stored below freezing temperature (<0 • C) [11][12][13] or significantly inhibited when the plants grew under sunshade at photoperiod [14].
In the present study, the alterations of transcripts in A. sinensis with three different treatments, including uncompleted (T1), completed (T2) and avoided vernalization (T3) were analyzed by combining isoforms with Illumina sequencing, and 104 DEGs involved in flowering and cold response were observed. The changes in soluble sugar and hormones (GA 3 , IAA and ABA) were examined via a spectrophotometer and HPLC.

Isoform Sequencing and Analysis
The total RNA of SAM was extracted using a Trizol reagent, the integrity of the RNA was determined using an Agilent 2100 Bioanalyzer and agarose gel electrophoresis, and the purity and concentration of the RNA were determined using a NanoDrop micro-spectrophotometer. mRNA was enriched by Oligo (dT) magnetic beads and then reverse-transcribed into cDNA using a Clontech SMARTer cDNA Synthesis Kit. The cDNA was amplified by PCR for 13 cycles to generate large-scale, double-stranded cDNA and then purified using AMPure XP beads to select cDNAs > 4 kb in size. SMRTbell library was constructed after the cDNAs were damage-repaired and end-repaired, and the sequencing adapters were ligated. The SMRTbell template was annealed to bind the primer and polymerase.
The raw reads of the cDNA library were sequenced on a Pacific Biosciences Sequel platform by Gene Denovo Biotechnology Co., Ltd. (Guangzhou, China), and then the raw reads were analyzed using a SMRT Link (V8.0.0) [21]. Briefly, high-quality circular consensus sequences (CCS) were extracted from the raw reads, and full-length non-chimeric (FLNC) reads were obtained by removing the primers, barcodes, poly (A) tail trimmings and concatemers. Then, the FLNC reads were clustered to generate the entire isoforms [22]. Similar FLNC reads were used, the minimap2 was used to cluster hierarchically to obtain the consistency sequence (unpolished consensus isoforms), the quiver algorithm was used to further correct the consistency of the sequence and then high-quality isoforms were obtained with a prediction accuracy ≥0.99. Finally, the full-length isoforms were analyzed and annotated against the databases using BLAST: NCBI non-redundant protein (Nr) (http://www.ncbi.nlm.nih.gov, accessed on 1 December 2021), Kyoto Encyclopedia of Genes and Genomes (KEGG) (http://www.genome.jp/kegg, accessed on 1 December 2021), Eukaryotic Orthologous Groups of proteins (KOG), Swiss-Prot (https://www.uniprot.org, accessed on 1 December 2021), Gene Ontology (GO) (http://www.geneontology.org/, accessed on 1 December 2021) and COG/KOG (http://www.ncbi.nlm.nih.gov/COG, accessed on 1 December 2021) with the BLASTx program (http://www.ncbi.nlm.nih.gov/BLAST/, accessed on 1 December 2021) at an E-value threshold of 1 × 10 −5 to evaluate the sequence similarity with genes of other species.

Illumina Sequencing and DEGs Analysis
The procedures of total RNA extraction, integrity and purity determination, as well as cDNA library construction, for the samples T1, T2 and T3 were conducted according to the "Section 2.2" description. After the samples were sequenced, the raw reads were filtered using a FASTQ system to obtain high-quality clean reads [23] with the following parameters: removing reads containing adapters, removing reads containing more than 10% unknown nucleotides (N), and removing low-quality reads containing more than 50% low quality (Q-value ≤ 20) bases. The paired-end clean reads were mapped to the fulllength isoforms using HISAT2. 2.4 to obtain the unique-and multiple-mapped reads [24] with "-rna-strandness RF" and other parameters set as a default.
The expression level of each transcript was normalized to the values of the fragments per kilobase of the exon model per million mapped reads (FPKM). Differential expression analysis of transcripts was performed using DESeq2 software [25] between different groups. The differential expression levels at T2 vs. T1 and T3 vs. T1 were determined with the criteria of the false discovery rate (FDR) < 0.05 and |log 2 (fold-change)| > 1. To date, the genome of A. sinensis has not been reported to be sequenced and assembled. The function of DEGs was analyzed using BLAST against the full-length isoforms of A. sinensis based on the Nr, KEGG, KOG and Swiss-Prot databases.
For GO enrichment analysis, all DEGs were mapped to GO terms in the GO database, gene numbers were calculated (FDR ≤ 0.05) for every term, and GO terms that were significantly enriched in DEGs compared with the full-length isoforms of A. sinensis were defined using a hypergeometric test. For the KEGG pathway enrichment analysis, significantly enriched metabolic pathways or signal transduction pathways in DEGs compared with the full-length isoforms of A. sinensis were defined with FDR ≤ 0.05 as a threshold.

qRT-PCR Validation
The procedures of total RNA extraction, integrity and purity determination for the samples T1, T2 and T3 were conducted according to the description labeled "Section 2.2". Primer sequences of selected DEGs (Table S1) were designed using an NCBI primerblast tool and synthesized by Sangon Biotech (Shanghai, China). First-strand cDNA was synthesized using a FastKing RT kit (KR116; Tiangen, Beijing, China) with one cycle at 42 • C for 15 min and then at 95 • C for 3 min. qRT-PCR gene expression was carried out using a SuperReal PreMix Plus (SYBR Green) (FP205; Tiangen, Beijing, China) with one cycle at 95 • C for 15 min, followed by 40 cycles at 95 • C for 10 s, 60 • C for 20 s and 72 • C for 30 s. Melting curve analysis was performed after incubation at 95 • C for 15 s, 60 • C for 1 min and 95 • C for 1 s. The actin gene was used as a reference control gene [9,15,16,26]. Herein, the cycle threshold (Ct) values and standard curves of the ACT gene at different volumes (0.25, 0.5, 1.0, 1.5, 2.0 and 3.0 µL) were built to correct the gene expression level ( Figures S3 and S4). The relative expression level (REL) of genes was calculated using a 2 −∆∆Ct method (Ct, cycle threshold value of target gene) [27].

Measurement of Soluble Sugar and Endogenous Hormones Contents
Soluble sugar content was measured using phenol-sulfuric acid [28]. Endogenous hormones (GA 3 , IAA and ABA) contents were measured using a HPLC method [29].

Statistical Analysis
In order to obtain the precise estimation of PCR efficiency, each experiment for qRT-PCR validation was performed with three biological replicates, along with three technical replicates. Statistical analysis was performed via ANOVA and Duncan multiple comparison tests, and SPSS 22.0 was the software package used with p < 0.05 as the basis for statistical differences.

Full-Length Isoform Analysis
A total of 1,031,219 high-fidelity reads were extracted after 33 full passes of raw reads, 76,840 polished high-quality isoforms were obtained using a Quiver calculation and 64,117 full-length isoforms were generated after the FLNC reads were clustered ( Figure S5). A total of 61,241 isoforms were annotated on KEGG (60,792), KOG (42,645), Nr (61,161) and Swiss-Prot (51,070) databases ( Figure 1A), and the top 10 species in terms of distribution against Nr, were Daucus carota, Apium graveolens, Nyssa sinensis, Angelica sinensis, Petroselinum crispum, Camellia sinensis, Mikania micrantha, Vitis vinifera, Oliveria decumbens and Beta vulgaris ( Figure 1B). Since the genome of A. sinensis has not been sequenced, the isoforms are needed to compare with other species and the limited genes of A. sinensis in the NR database.

Illumina Sequencing of T1, T2 and T3
To reveal the molecular mechanisms responsible for regulating the flowering of A. sinensis during vernalization, a comparison of gene transcription in response to different temperatures and durations was performed. In this study, 9.14, 9.19 and 8.63 million raw reads were generated for T1, T2 and T3, respectively. After raw data filtering, 9.12, 9.16 and 8.61 million clean reads were collected. After mapping on the isoforms, 5.33, 5.50 and 5.23 million unique mapped reads, as well as 2.27, 2.31 and 2.05 million multiple mapped reads, were obtained from the T1, T2 and T3, respectively. Meanwhile, the exon rate of all the three treatments reached 100% (Table 1).  Figure S8). The cluster heat map of the DEGs at T2 vs. T1 and T3 vs. T1 was shown in Figure 3, which provides further analysis such as identification of the gene functions or gene response analysis.

GO and KEGG Enrichments of DEGs
Based on the GO system, the DEGs were classified into three ontologies including biological process (BP), cellular component (CC) and molecular function (MF) ( Figure S9). Based on the KEGG database, at T2 vs. T1, the 4212 DEGs were enriched into 127 pathways, such as plant hormone signal transduction, starch and sucrose metabolism, and phenylpropanoid biosynthesis; at T3 vs. T1, the 5301 DEGs were enriched into 128 pathways, such as biosynthesis of secondary metabolites, the MAPK signaling pathway, and plant and starch and sucrose metabolism ( Figure S10).

Genes Involved in Sucrose Metabolism
Twelve genes involved in sucrose metabolism were observed to be differentially expressed ( Figure 5B, Table S4), with genes encoding trehalose-phosphate synthase (TPS1, TPS5, TPS7 and TPS10), sucrose synthase (SUS2 and SUS3), sucrose-phosphate synthase (SPS and SPS1), phosphoglucomutase (PGMP and PGM1) and beta-amylase (BAM1 and BAM3). Three genes were confirmed using qRT-PCR: the RELs of these three genes (TPS1, SUS3 and SPS) were consistent with their FPKM values at T2 vs. T1 and T3 vs. T1 ( Figures 5B and 7A). The soluble sugar content exhibited a 0.77-fold decrease at T2 vs. T1 and a 1.13-fold increase at T3 vs. T1 ( Figure 7B). Here and below, "*" represents a significant difference at the p < 0.05 level between T2 vs. T1 and T3 vs. T1 for the same gene; the use of different letters represents a significant difference (p < 0.05) at different temperature treatments.

DEGs Involved in Circadian Clock
Two genes involved in circadian clock were observed to be differentially expressed ( Figure 5D, Table S6), including a transcription factor part of a circadian clock (ELF3) and cold-regulated protein (COR27). These two genes were confirmed using qRT-PCR, and the RELs were consistent with their FPKM values at T2 vs. T1 and T3 vs. T1 (Figures 5D and 9A).

DEGs Involved in Downstream Floral Integrators and Meristem Identity
Fifteen genes involved in downstream floral integrators and meristem identity were observed to be differentially expressed ( Figure 5E, Table S7); these included floral integrators (FTIP3, FTIP7 and SOC1) and floral meristem identity (AGL65, SPL1, SPL6, AP2, AP2-1, AP2-2, AP2-3, AIL5, ANT, At2g41710, TEM1 and SOK2). Six genes were confirmed by qRT-PCR, and the RELs of six of these genes (FTIP3, FTIP7, SOC1, AGL65, SPL1 and SOK2) were consistent with their FPKM values at T2 vs. T1 and T3 vs. T1 (Figures 5E and 9B). Here and below, "*" represents a significant difference at the p < 0.05 level between T2 vs. T1 and T3 vs. T1 for the same gene; the use of different letters represents a significant difference (p < 0.05) at different temperature treatments. Here and below, "*" represents a significant difference at the p < 0.05 level between T2 vs. T1 and T3 vs. T1 for the same gene; the use of different letters represents a significant difference (p < 0.05) at different temperature treatments.

Discussion
Vernalization is the process whereby flowering is promoted by prolonged exposure to a cold treatment given to a fully hydrated seed or a growing plant; without vernalization, plants show delayed flowering or remain vegetative. Meanwhile, vernalization can be lost at high temperatures and avoided below freezing temperatures [33]. Extensive studies have reported that vernalization suppresses the expression of genes encoding the repressors of flowering [34]; moreover, the gene network of flowering time has been mapped with known genetic and epigenetic regulators in the model plant Arabidopsis thaliana [31,32]. For A. sinensis, the effective temperature range for vernalization is from 0 to 5 • C with a duration 57 to 85 d, and the vernalization can be effectively avoided with exposure of seedlings to temperatures from −2 to −12 • C [10][11][12][13]. However, gene-regulatory networks involved in the flowering of A. sinensis during the vernalization have not been revealed. By storing seedlings at 0 • C for 14 d (uncompleted vernalization, T1) and 60 d (completed vernalization, T2), as well as at −3 • C for 125 d (avoided vernalization, T3), a total of 4212 and 5301 genes were differentially expressed, and 151 and 155 genes were involved in flowering at T2 vs. T1 and T3 vs. T1, respectively. Based on their biological functions, 104 coexpressed genes were classified into 6 categories: FLC expression, sucrose metabolism, hormone response, circadian clock, downstream floral integrators and meristem identity, and cold response.

Genes Involved in FLC Expression
FLC is a MADS-box transcriptional regulator that acts as a potent repressor of flowering [35]. Several genes inhibit FLC expression; among these, 11 genes encoding proteins are involved in vernalization response, including VIN3-like proteins (VIL1 and VIL2) that inhibit the expression of FLC and FLM, which are associated with an epigenetically silenced state and with acquisition of competence to flower [36][37][38][39], and VRN1-like (VRN1-L), which may act as transcriptional repressor of FLC [40]. The genes that inhibit FLC expression also include the genes involved in the autonomous pathway, including proteins that control flowering time, FCA and FPA, which decrease FLC expression and act redundantly with each other to prevent the expression of distally polyadenylated antisense RNAs at the FLC locus [41][42][43]; flowering locus K (FLK), which represses FLC expression and post-transcriptional modification [44,45]; another protein controlling flowering time, FY, which decreases FLC expression and is required for the negative autoregulation of FCA expression [46,47]; and FLOWERING LOCUS D (FLD) and LUMINIDEPENDENS (LD), which decrease FLC expression via chromatin modification [48][49][50]. Finally, there are two other genes that inhibit FLC expression: chromo domain protein LHP1 (LHP1), which is a structural component of heterochromatin involved in gene repression via methylating H3-K9 [51], and MAF1 homolog (MAF1), which is an element of the TORC10 signaling pathway and represses RNA polymerase III transcription [52].
For genes activating FLC expression, there are nine genes encoding proteins involved in the PAF1 complex: LEO1 homolog (VIP4), which is involved in histone modifications and is required for the FLC expression [53,54]; histone-lysine N-methyltransferases (ATX2, ATX4 and ATX5), which dimethylate H3K4me2 and are involved in epigenetic regulation of FLC and FT [55][56][57]; genes involved in the SWR1 complex, namely actin-related protein 6 (ARP6) and PHOTOPERIOD-INDEPENDENT EARLY FLOWERING 1 (PIE1), which are associated with transcriptional regulation of selected genes (e.g., FLC) via chromatin remodeling; ARP6, which is required for the activation of FLC and FLC/MAF genes expression through both histone H3 and H4 acetylation and methylation [58,59]; PIE1, which is required for the reactivation of FLC [60,61]. There are also genes involved in the FRIGIDA-like protein: FRIGIDA-like protein 3 (FRL3), which trimethylates H3K4, increasing FLC expression [62]; FRIGIDA-ESSENTIAL 1 (FES1), which acts cooperatively with FRI or FRL1 to promote FLC expression [63,64]; and SUPPRESSOR OF FRI 4 (SUF4), which recruits the FRI-C complex to the FLC promoter and is required for FRI-mediated FLC activation and maintains high levels of FLC expression [63,65].
There are also genes interacting with FLC expression: these include two genes involved in PcG protein: histone-lysine N-methyltransferase CLF (CLF) and PcG protein FIE1 (FIE2), which are required to maintain the transcriptionally repressive state of homeotic genes by methylating H3-K27, rendering chromatin heritably changed in its expressability [66,67]. According to the RELs of genes calculated with the FPKM value and validated using qRT-PCR, the genes involved in inhibiting FLC expression were up-regulated, while the genes involved in activating FLC expression were down-regulated at T2 vs. T1 (during vernalization); however, these genes showed the opposite expression levels at T3 vs. T1 (freezing temperature avoided vernalization) ( Figure 6; Table S3). These findings are consistent with previous studies with higher flowering rates at T2 with lower flowering rates at T3 compared with T1 ( Figure S2).

Genes Involved in Sucrose Metabolism
For vernalization to occur, sources of energy (sugars) and carbohydrate metabolism are required [33,68]. Twelve genes were found to be involved in sucrose metabolism: trehalosephosphate synthases (TPS1, TPS5, TPS7 and TPS10), which are required for vegetative growth and transition to flowering by regulating starch and sucrose degradation [69][70][71]; sucrose synthases (SUS2 and SUS3), which provide UDP glucose and fructose for various metabolic pathways [72]; sucrose-phosphate synthases (SPS and SPS1), which regulate sucrose synthesis from UDP glucose and fructose-6-phosphate [73]; phosphoglucomutases (PGMP and PGM1), which participate in both the breakdown and the synthesis of glucose [74]; and beta-amylases (BAM1 and BAM3), which play roles in circadian-regulated starch degradation and maltose metabolism [75]. According to the RELs, the genes involved in sucrose and starch degradation were up-regulated during vernalization (T2) and downregulated at freezing temperature (T3) compared with T1, while the genes involved in sucrose and starch biosynthesis showed the opposite expression levels ( Figure 7A; Table S4). This differential expression of genes is consistent with the change in soluble sugar contents at T1, T2 and T3 ( Figure 7B).

Genes Involved in Hormone Response
The GA pathway is required for early flowering by promoting the expression of the LFY gene; meanwhile, other growth hormones (e.g., ETH, IAA and CTK) can either inhibit or promote flowering [33]. Here, there were 18 genes that were found to be involved in the hormone response, including the GA receptor GID1B (GID1B), which functions as soluble GA receptor interacting with specific DELLA proteins, known as repressors of GA-induced growth and flower development [76]; ETH-responsive transcription factors (RAP2-3 and RAP2-7), which negatively regulate the transition to flowering time and cause a delay in flowering time [77]; IAA-amido synthetase (GH3.6), which is involved in auxin signal transduction [78]; auxin transporter-like protein 2 (LAX2), which is involved in proton-driven auxin influx and basipetal auxin transport [79,80]; Aux/IAA proteins (IAA8, IAA9, IAA12, IAA13 and IAA27), which function as repressors of early auxin response genes at low auxin concentrations by forming heterodimers with auxin response factors (ARFs) [81]; ARFs (ARF1, ARF4, ARF5, ARF6 and ARF19), which act as transcriptional repressors (e.g., IAA2, IAA3 and IAA7) by forming heterodimers with Aux/IAA proteins and promoting flowering [82,83]; and CTK response, i.e., histidine kinase 3 (AHK3), which the cytokinin-dependent flower-development-regulation pathway requires [84]; histidine-containing phosphotransfer protein 1 (AHP1), which functions as a two-component phosphorelay mediator between cytokinin sensor histidine kinases and response regulators (B-type ARRs) [85]; and two-component response regulator ARR2 (ARR2), which functions as a response regulator involved in the His-to-Asp phosphorelay signal transduction system [86]. Based on the RELs, the genes involved in promoting flowering were up-regulated during vernalization but down-regulated at freezing temperatures ( Figure 8A; Table S5). Meanwhile, the change in GA 3 and IAA contents was also shown to be higher at vernalization compared with freezing temperatures (Figure 8B,C).

Genes Involved in the Circadian Clock
Flowering is often triggered when plants are exposed to appropriate day lengths, which requires the circadian clock to measure the passage of time [87]. Recent studies have found that the circadian clock is also controlled by temperatures [88]. Here, two genes were found to be involved in the circadian clock: EARLY FLOWERING 3 (ELF3), which is a transcription factor part of a circadian-clock input pathway and can regulate the initiation of flowering [89,90], and cold-regulated protein (COR27), which is a negative regulator of freezing tolerance that, together with COR28, is involved in central circadian clock regulation and in flowering promotion by binding to the chromatin of clock-associated evening genes (e.g., ELF4) [91,92]. According to the RELs, the genes (ELF3 and COR27) favoring the initiation of flowering were up-regulated ( Figure 9A; Table S6), which accelerates the formation of leaves and eventually the physiological age during vernalization [33].

Genes Involved in the Downstream Floral Integrators and Meristem Identity
Flower formation occurs at the SAM and is a complex morphological event that is regulated by several genes [33]. Here, 15 genes were found to be involved in downstream floral integrators and meristem identity. Three of these genes were involved in floral integrators: FT-interacting protein 3 (FTIP3), which is required for the proliferation and differentiation of shoot stem cells in SAM [93]; FT-interacting protein 7 (FTIP7), which promotes nuclear translocation of the transcription factor homeobox 1, directly suppressing the auxin biosynthetic gene YUCCA4 [94]; and SOC1, which integrates signals from the photoperiod, vernalization and autonomous floral induction pathways [95].
Ten genes are involved in the floral meristem identity: agamous-like MADS-box protein AGL65 (AGL65), which forms a heterodimer with the MADS-box protein AGL104 [96]; squamosa promoter-binding-like proteins (SPL1 and SPL6), which bind specifically to the consensus nucleotide sequence of the AP1 promoter [97]; APETALA 2 and -like proteins (AP2, AP2-1, AP2-2, AP2-3, AIL5, ANT and At2g41710), which are involved in initiation and development of organs, including floral organs [98][99][100][101]; AP2/ERF and B3 domain-containing transcription repressor TEM1 (TEM1), which is a transcriptional repressor of flowering time in plants that prefer long days and acts directly upon FT expression [102]; and SOK2, which can influence cell division orientation to coordinate cell polarization, which is fundamental for tissue morphogenesis in multicellular organisms [103][104][105]. According to the RELs, the genes involved in floral integrators and meristem identity were up-regulated during vernalization but down-regulated at freezing temperatures ( Figure 9B; Table S7). These findings further indicate that vernalization enhances the early flowering of A. sinensis.

Genes Involved in the Cold Response
Vernalization cannot be completed in the context of short exposures to cold, which might occur during cases of fluctuating temperatures, after which plants may face and adapt to low temperatures [34]. Here, 35 genes were found to be involved in the cold response, including 2 genes involved in ABA response: ABA receptors PYL3 and PYL4 (PYL3 and PYL4), which act as positive regulators of tolerance to cold stress [106]. Genes involved in the cold response also include those involved in the ETH response: ethylene-responsive transcription factors (ERF3, ERF5, ERF010, ERF011, ERF13, ERF110 and ERF118), which regulate gene expression by stress factors and by components of stress signal transduction pathways [107,108], and ETH receptor 2 (ETR2), which acts as a redundant negative regulator of ETH signaling [109]. They also include two genes involved in cold-regulated plasma membrane proteins: cold-regulated 413 plasma membrane proteins (COR413PM1 and COR413PM2), which are up-and down-regulated, respectively, in response to low temperatures [110]. In addition, they include six genes involved in cold-responsive protein kinase: cold-responsive protein kinase 1 (CRPK1), which is a negative regulator of freezing tolerance [111]; mitogen-activated protein kinases (MPK3, MPK9, MPK16 and YDA), which are associated with the ABA-activated signaling pathway in response to oxidative stress and freezing [112,113]; serine/threonine-protein kinase SRK2E (SRK2E), which is an activator of the ABA signaling pathway [114]. They also included 15 genes that were found to be involved in protein phosphatases (PP2CA, HAB1, HAB2, KAPP, PPC6-1, PPC6-7, PP2C27, PP2C38, At2g30020, Os02g0799000, At3g15260, At3g16560, At3g62260, Os06g0651600 and At4g31860), which regulate numerous ABA responses such as cold stress [115][116][117]. In addition, transcription factor TGA2.2 (TGA2.2) is involved in the defense response [118], and cell division cycle 5-like protein (CDC5) is involved in mRNA splicing and cell cycle control and may also play a role in the response to DNA damage [119,120]. According to the RELs, most of the genes involved in cold stress were up-regulated ( Figure 10A; Table S8), and the ABA contents increased in response to vernalization and freezing temperatures ( Figure 10B), which trigger the stress response for the seedlings adapting to low temperatures.

Proposed Regulatory Networks of Flowering Genes in A. sinensis during Vernalization
Based on the functional analysis and the regulatory pathways of flowering genes in the model plant Arabidopsis [30][31][32], a schematic representation of the proposed regulatory networks of flowering genes in A. sinensis during vernalization was created and is shown in Figure 11. Briefly, during the vernalization, four pathways that regulate flowering time were observed: the vernalization pathway, the autonomic pathway, the age pathway and the GA (hormone) pathway. In the vernalization pathway, the expression of FLC is inhibited by the genes VIL1, VIL2, VRN1-L, LHP1 and MAF1 but activated by the genes encoding the PAF1 complex (e.g., VIP4, ATX2 and ATX), the SWR1 complex (ARR6 and PIE1) and FRIGIDA-like proteins (FRL3, FES and SUF4) and regulated by proteins of the polycomb group (CLF and FIE2). In the autonomic pathway, the expression of FLC is inhibited by the genes involved in RNA binding (FCA, FPA and FLK), the RNA process (FY) and chromatin modification (FLD and LD). In the age pathway, the expression of SOC1 is promoted by the genes encoding SPLs (SPL1 and SPL6), which can be positively regulated by the genes involved in sucrose metabolism (e.g., TPS1, SUS2 and SPS). In the GA (hormone) pathway, the expression of SOC1 is promoted by the genes involved in the CTK response (AHK3, AHP1 and ARR2) but inhibited by the genes involved in ETH response (RAP2-3 and RAP2-7) and IAA (e.g., GH3.6, LAX2 and IAA8); the expression of AP1 is promoted by the DELLA protein, which is inhibited by the genes involved in the GA response (GID1B). Additionally, the expression of FT (FTIP3 and FTIP7) is promoted by the genes involved in the circadian clock (ELF3 and COR27), as well as AP2 and -like proteins (e.g., AP2, AP2-1 and AIL5), but inhibited by the gene TEM1. The expression of SOC1 is also promoted by the gene AGL65. The coordinated expression of these genes during vernalization confers the transition of seedlings from vegetative growth to flowering of A. sinensis. Figure 11. Schematic representation of the proposed regulatory networks of flowering genes in A. sinensis during vernalization. The arrows indicate a promotion, with T-ends with red color indicating an inhibition, round dots at both ends mark an interaction without a known direction, and dashed lines indicate an indirect interaction.

Conclusions
The DEGs observed in A. sinensis during vernalization strongly suggest that transcriptionbased regulation occurs for the transition of seedlings from vegetative growth to flowering. The expression level of genes involved in flowering and the cold response during vernalization were almost consistent with changes in sugars and hormone contents. There are four pathways of genes that are required for regulating flowering, including the vernalization pathway, the autonomic pathway, the age pathway and the GA (hormone) pathway. While genes involved in flowering during vernalization have been mapped here, additional studies are required to determine the causative role of these genes in promoting or inhibiting the expression of the central genes (e.g., VINs, SOC1 and FT).

Informed Consent Statement: Not applicable.
Data Availability Statement: The datasets are available at NCBI, with BioSample accession SAMN24046640 to SAMN24046648; SRA accession SRR17235563 to SRR17235571 (9 objects) with every treatment three replicates; and the BioProject's metadata are available at https://www.ncbi. nlm.nih.gov/bioproject/PRJNA789039, accessed on 14 December 2021.

Conflicts of Interest:
The authors declare no conflict of interest.