Transcriptomic Analysis Reveals LncRNAs Associated with Flowering of Angelica sinensis during Vernalization

Angelica sinensis is a “low-temperature and long-day” perennial plant that produces bioactive compounds such as phthalides, organic acids, and polysaccharides for various types of clinical agents, including those with cardio-cerebrovascular, hepatoprotective, and immunomodulatory effects. To date, the regulatory mechanism of flowering under the photoperiod has been revealed, while the regulatory network of flowering genes during vernalization, especially in the role of lncRNAs, has yet to be identified. Here, lncRNAs associated with flowering were identified based on the full-length transcriptomic analysis of A. sinensis at vernalization and freezing temperatures, and the coexpressed mRNAs of lncRNAs were validated by qRT-PCR. We obtained a total of 2327 lncRNAs after assessing the protein-coding potential of coexpressed mRNAs, with 607 lncRNAs aligned against the TAIR database of model plant Arabidopsis, 345 lncRNAs identified, and 272 lncRNAs characterized on the SwissProt database. Based on the biological functions of coexpressed mRNAs, the 272 lncRNAs were divided into six categories: (1) chromatin, DNA/RNA and protein modification; (2) flowering; (3) stress response; (4) metabolism; (5) bio-signaling; and (6) energy and transport. The differential expression levels of representatively coexpressed mRNAs were almost consistent with the flowering of A. sinensis. It can be concluded that the flowering of A. sinensis is positively or negatively regulated by lncRNAs, which provides new insights into the regulation mechanism of the flowering of A. sinensis.


Plant Material
The seedlings (root-shoulder diameter 0.4-0.5 cm; Figure S1) of Angelica sinensis (cultivar Mingui 1) were stored at 0 (vernalization temperature) and −3 • C (freezing temperature), respectively. After storage at 0 • C for 14 (T1) and 60 days (T2), as well as at −3 • C for 125 days (T3), the shoot apical meristem (SAM) was cut from the root shoulder of the seedlings for transcriptomic analysis and qRT-PCR validation. Three biological replicates were performed for each treatment of T1, T2, and T3. Herein, the treatment of T1, T2, and T3 represents uncompleted, completed, and avoided vernalization, respectively, based on the EBF rate ( Figure S2) when the stored seedlings were cultivated and grown in a long-day condition.

Full-Length Isoform Sequencing and Analysis
Total RNA of the SAM samples was extracted using Trizol reagent (Omega Bio-Tek, Norcross, GA, USA). The integrity of the RNA was determined using an Agilent 2100 Bioanalyzer (Agilent Technol., California, CA, USA) and agarose gel electrophoresis, and the purity and concentration of the RNA were determined using a microspectrophotometer (NanoDrop Technol., Wilmington, DE, USA). The high-quality RNAs were sequenced on a Pacific Biosciences Sequel platform (Gene Denovo Biotechnology Co., Ltd., Guangzhou, China). Raw reads of cDNA library were analyzed using a SMRT Link (V8.0.0) [26]. Briefly, high-quality CCS were extracted from the subreads BAM file; the integrity of transcripts (full-length sequences) was judged based on whether CCS reads contained primers (5 and 3') and polyAs; then, FLNC reads were generated by removing primers, barcodes, and polyAs; finally, FLNC reads were assembled to obtain the entire isoform [27].

Analysis of Long Noncoding RNAs (lncRNAs)
Isoforms that were not annotated against the four databases-NR, Swiss-Prot, Kyoto KEGG, and KOG-were used for the analysis of lncRNAs. The isoform that was assessed as a noncoding transcript by both CNCI and CPC software was finally confirmed as a lncRNA [28,29].

Characterization of LncRNAs
To date, the genome of A. sinensis has not been sequenced. Thus, the lncRNA analysis of A. sinensis was performed via a BLAST search with an E-value cut-off of ≤1 × 10 −5 against the known lncRNAs from the TAIR database (https://www.arabidopsis.org accessed on 30 March 2022) [30]. The function of lncRNAs was annotated based on their coexpression mRNAs [31][32][33]. Herein, the biological functions of the coexpressed mRNAs were searched on the UniProt database (https://www.uniprot.org accessed on 30 March 2022).

qRT-PCR Validation
Based on the coding sequences (CDS) of coexpressed mRNA of lncRNA, 49 primer sequences of representatively coexpressed mRNAs (Table 1) were designed using the NCBI primer-blast tool. First-strand cDNA synthesis and qRT-PCR reaction were carried out using SuperRealPreMix Plus (FP205; Tiangen Biotech., Beijing, China) according to the manufacturer's instructions; specifically, the cDNA was synthesized successively with one cycle (95 • C, 15 min) and 40 cycles (95 • C, 10 s; 60 • C, 20 s; and 72 • C, 30 s), and the qRT-PCR reaction was incubated successively at 95 • C for 15 s, 60 • C for 1 min and 95 • C for 1 s. The Actin (ACT) gene was used as a reference control gene with forward: TGGTATTGTGCTGGATTCTGGT and reverse: TGAGATCACCACCAGCAAGG (amplicon size 109 bp) [34]. 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) was built to correct the gene expression level ( Figure S3 and Figure S4), and the expression levels of the 49 candidate genes and their standard deviations for every variant were added to the Supplementary Materials (Table S1). The REL of coexpressed mRNAs was calculated using the 2 − Ct method [35] according to the following formula.

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 [36]. A t-test in SPSS 22.0 was performed for independent experiments, with p < 0.05 as the basis for statistical differences.

LncRNAs Analysis
In total, 2327 lncRNAs were obtained after assessing the protein-coding potential of coexpressed mRNAs based on the two software programs, CNCI and CPC ( Figure 1A), with 607 genes aligned against the known lncRNAs from the TAIR database of model plant Arabidopsis (Figure 1B), and 345 lncRNAs with coexpressed mRNAs of A. sinensis identified ( Figure 1C) based on the SwissProt database. Based on the biological functions, the 272 characterized lncRNAs ( Figure 1D) were divided into six categories: chromatin, DNA/RNA and protein modification (29); flowering (36); stress response (24); metabolism (117); biosignaling (23), and energy and transport (43) ( Figure 1E). The base sequences of the 272 lncRNAs are shown in Table S2.

LncRNAs Linked with Flowering and Expression Levels of Their Coexpressed mRNAs
In total, 36 lncRNAs were linked with flowering based on the biological functions of their coexpressed mRNAs, with 12 lncRNAs directly associated with flowering, namely SRR1, PHL, PHYA, AGL62, AGL79, ATJ3, BBX29, CLE13, CLE44, MXC17.10, At1g06515, and BHLH30 (Table 3), and 24 lncRNAs indirectly associated with flowering, such as cell division, embryo development, and cell wall organization (Table S3). The expression Figure 2. The expression levels of coexpressed mRNAs of lncRNAs linked with chromatin, DNA/RNA and protein modification in A. sinensis at T2 vs. T1 and T3 vs. T1, as determined by qRT-PCR. T1, T2, and T3 represent uncompleted, completed, and avoided vernalization, respectively. The "*" represents a significant difference at p < 0.05 level between T2 vs. T1 and T3 vs. T2 for the same gene.

Discussion
Vernalization is a process considered to be an epigenetic switch whereby flowering is promoted by prolonged exposure to cold (0 to 10 • C); meanwhile, it can be lost at high temperatures or avoided below freezing temperatures [37,38]. Epigenetic regulation involves diverse molecular mechanisms including chromatin remodeling, DNA methylation, histone modification, histone variants, and ncRNAs [39]. Studies on Brassica rapa found that 127 differentially expressed lncRNAs were coexpressed with 128 differentially expressed genes, indicating that lncRNAs played an important role during vernalization [40]. In this study, 272 characterized lncRNAs were identified from A. sinensis and divided into six categories, namely (1) chromatin, DNA/RNA, and protein modification; (2) flowering; (3) stress response; (4) metabolism; (5) biosignaling; and (6) energy and transport, based on their coexpressed mRNAs.
FLC is a MADS-box transcriptional regulator that acts as a potent repressor of flowering [38]. In Arabidopsis, the epigenetic silencing of the floral repressor gene FLC is a well-characterized key event of vernalization [41]. In this study, 29 lncRNAs linked with chromatin, DNA/RNA, and protein modification were identified in A. sinensis during vernalization. For the chromatin modification, two coexpressed mRNAs, HMGB2 and HMGB3, are involved in binding preferentially double-stranded DNA and up-regulated in response to cold stress [42]. For the DNA/RNA modification, At1g05910 is involved in DNA demethylation and the negative regulation of chromatin silencing [43]; RID2 is involved in rRNA methyltransferase activity [44]; and At4g26600 is involved in RNA methylation [45]. For protein modification, 24 lncRNAs were involved and the roles of nine coexpressed mRNAs were represented as follows. H2AV plays a central role in regulating transcription, repairing DNA, replicating DNA, and stabilizing the nucleus chromosome [46]; At2g28720 and HTR2 are involved in compacting DNA into chromatin [47,48]; HOS15 promotes the deacetylation of histone H4 [49]; REF6 is involved in demethylating 'Lys-27' of histone H3, regulating flowering time by repressing FLC expression and interacting with the NF-Y complex to regulate SOC1 [50,51]; SKP1A and ASK21 are involved in ubiquitination and form an SCF E3 ubiquitin ligase complex together with CUL1, RBX1, and an F-box protein [52,53]; and SRK2G and SRK2H are involved in protein phosphorylation [54,55]. In previous studies, HMGB2 and HMGB3 were up-regulated in response to cold stress but down-regulated in response to drought and salt stresses [42]; H2AV, At2g28720 and HTR2 exhibited a high level in response to osmotic and drought stresses [56]; HOS15 was found to act as a repressor of cold stress-regulated gene expression and played a role in gene regulation for plant acclimation and tolerance to cold stress [49]; SKP1A and ASK21 were overexpressed in the host stress response [57]; and SRK2G and SRK2H were found to be positive regulators in stress responses such as drought, salt, and cold [54,55,58]. Currently, although most of the lncRNAs have been reported to be involved in the stress response in many plants, their roles in the regulation of flowering time have been studied in model plants [23]. Previous studies have found that FLC in Arabidopsis is epigenetically regulated by lncRNAs COOLAIR and COLDAIR [59][60][61]. In this study, a lower expression level was noted for almost all coexpressed mRNAs of lncRNAs involved in chromatin, DNA/RNA, and protein modification at vernalization (T2, 0 • C) compared with freezing temperature (T3, −3 • C), as well as down-regulation at T2, indicating that these lncRNAs participate in epigenetic silencing by transferring euchromatin to heterochromatin and confer early bolting and flowering of A. sinensis. Representatively, the down-regulation of mRNAs HMGB2 and HMGB3 during vernalization transfers the heterochromatin to the euchromatin [42], which generates the ability of flowering; in contrast, their up-regulation at freezing temperatures inhibits flowering by keeping the heterochromatin, which indicates that their coexpressed lncRNAs play positive roles in regulating the flowering of A. sinensis. The down-regulation of mRNA REF6 below 0 • C weakens the demethylation of histone H3 and delays flowering time by inducing FLC expression [50,51]; meanwhile, there is an increased expression level with decreased temperatures, which indicates that this coexpressed lncRNA plays a negative role in regulating the flowering of A. sinensis.
Flower formation occurs at the SAM and is a complex morphological event that is required not only for the circadian clock to measure the passage of time but also the regulation of meristem identity genes [42]. For the 6 representative coexpressed mRNAs of the 12 lncRNAs directly linked with flowering, SRR1 is involved in a circadian clock input pathway and regulation of the expression of clock-regulated genes such as CCA1 and TOC1 [62]; PHL is involved in the circadian rhythm and the regulation of the timing of transition [63]; PHYA is involved in the regulation of flowering time and expression of its own gene as negative feedback [64]; AGL62 is required for promoting the nuclear proliferation of early endosperm [65]; AGL79 is involved in positively regulating the transition of the meristem from the vegetative to reproductive phase [66]; and ATJ3 plays a continuous role in plant development, such as in photoperiodism, flowering, and positive regulation of flower development [67]. Previous studies found that mRNAs (e.g., miR156, miR169 and miR172) play a crucial role in developmental processes in rice, wheat, and maize, especially in the formation of the floral meristem, with miR172 controlling AP2-like genes [23,68,69]. Studies on the flowering of Chenopodium quinoa found that pivotal flowering homologs, including photoreceptor genes PHYA and CRY1, as well as genes associated with florigen-encoding genes (FT and TWIN SISTER of FT) and circadian clockrelated genes (ELF3, LHY, and HY5), were specifically affected by night-break and competed with the positive-and negative-flowering lncRNAs [70]. In this study, down-regulation of all the coexpressed mRNAs involved in circadian clock and meristem identity genes was observed, which can be considered acceptable and reasonable, because these mRNAs are often highly expressed at the plant development stage (at photoperiod), while their expression levels were examined during vernalization. In addition, increased expression with decreased temperatures indicates that their coexpressed lncRNAs play negative roles in regulating the flowering of A. sinensis.
The expression of numerous lncRNAs has been demonstrated to be significantly affected by various stresses [23,30]. During vernalization, plants have to face and adapt to low temperature [38]. To date, extensive studies have reported that lncRNAs participate in defense responses associated with plant immunity and adaptation to the environment [22]. Heat-responsive lncRNAs have been found to be differentially expressed in Brassica juncea, and cold-responsive lncRNAs have been identified in grape and Arabidopsis [71][72][73]. lncR-NAs could regulate HSP family genes (HSP82 and HSP83) in response to heat stress in Populus x canadensis Moench, and HSP18.1 in response to Cd stress Betula platyphylla [74,75]. For the eight representative coexpressed mRNAs of the 14 lncRNAs linked with the temperature response, ACBP6 confers resistance to cold and freezing [76]; ENO2 acts as a positive regulator in response to cold stress [77]; ADH1 is required for survival and acclimation in response to abiotic stresses (e.g., cold, salt, and drought) [78,79]; CSP2 contributes to the enhancement of cold and freezing tolerance [80]; and HSP17.8, HSP70-3, HSP70-10, and HSP90-3 play vital roles in adapting to biotic and abiotic stresses [81,82]. In this study, increased expression for cold-tolerated mRNAs (ACBP6, ENO2, ADH1, and CSP2), and decreased expression for heat-tolerated mRNAs (HSP17.8, HSP70-3, HSP70-10, and HSP90-3), were observed with decreased temperatures, which indicates that their coexpressed lncRNAs play positive roles in adapting to low temperatures.
For vernalization to occur, sources of energy (sugars) and carbohydrate metabolism are required [37]. In recent years, the roles of lncRNAs in regulating metabolism in cancer, insulin, and chicken have been reported [83][84][85], while, in plants, studies are still limited. For the five representative coexpressed mRNAs of the 19 lncRNAs linked with carbohydrate metabolism, CSY4 is involved in the synthesis of socitrate from oxaloacetate [86]; GAPCP2 plays a critical role in glycolysis and exhibits up-regulation under drought stress [87,88]; At3g52990 is involved in the synthesis of pyruvate from D-glyceraldehyde 3-phosphate [89]; PGM2 participates in the synthesis of glucose [90]; and BAM1 is required for starch breakdown [91]. In this study, the differential expression of these coexpressed mRNAs regulating sucrose and starch metabolism provided energy for the morphogenesis of seedlings and adaptation to low temperatures during vernalization. Representatively, the decreased expression of metabolite-synthesized mRNAs (CSY4, At3g52990, and PGM2), and increased expression of energy-produced mRNA GAPCP2 and metabolite-degraded mRNA BAM1, were observed with decreased temperatures, which indicates that their coexpressed lncR-NAs play positive roles in carbohydrate metabolism.
Endogenous hormones such as gibberellin, auxin, cytokinin, brassinosteroid and abscisic acid can either inhibit or promote flowering [37]. In previous studies, a pre-miRNA of miR393 was identified in Brassica rapa during vernalization, and the overexpression of an miR393-resistant form of TIR1 (mTIR1) could enhance auxin sensitivity, thus leading to pleiotropic effects on plant development [92]. For the 8 representative coexpressed mRNAs of the 13 lncRNAs linked with hormone signaling, ARF1 is involved in the recruitment of COPI and GDAP1 to membranes and various auxin-dependent developmental processes [93]; CUL1 participates in forming a SCF complex together with SKP1, RBX1, and a F-box protein and is involved in floral organ development, the auxin signaling pathway and ethylene signaling [94]; SOFL4 and SOFL5 are involved in cytokinin-mediated development [95]; GRF2 and GRF11 participate in the brassinosteroid (BR)-mediated signaling pathway [96,97]; ERF3 is found to be differentially expressed in response to stresses and also regulates other ERFs [98]; and SF1 is required for development and is involved in the alternative splicing of FLM pre-mRNA [99]. In this study, the differential expression of these mRNAs involved in hormone signaling also played certain roles in regulating the flower-bud differentiation of seedlings and cold tolerance during vernalization. In previous studies, GRF2 and ERF3 were found to be down-regulated, while GRF11 was upregulated in response to cold stress [100,101]; here, contrary findings for mRNAs GRF2 and GRF11 were observed with temperatures decreased, which indicate that their coexpressed lncRNAs may play negative roles in hormone signaling. The interaction between SF1 and FLM pre-mRNA controls flowering time in response to temperature [102]; here, decreased expression of mRNA SF1 was observed with decreased temperatures, which indicates that this coexpressed lncRNA may play a positive role in hormone signaling.
Energy generation and the transport of energy, nutrients and metabolites play essential roles in growth and development and stress tolerance [103]. For the 8 representative coexpressed mRNAs of the 43 lncRNAs linked with energy and transport, PURU1 is involved in photorespiration and participates in preventing the excessive accumulation of 5-formyl tetrahydrofolate [104]; MES16 is involved in chlorophyll breakdown by its demethylation [105]; GPT1 is involved in NADPH generation through a series of processes including Glc6P transport, starch biosynthesis, fatty acid biosynthesis, and oxidative pentose phosphate [106]; ABCF5 belongs to the ABC transporter superfamily and is involved in protein transport [107]; SECA2 is involved in protein export coupling ATP hydrolysis [108]; VPS26A plays a role in vesicular protein sorting and is essential in endosome-to-Golgi retrograde transport [109]; CML19 is a potential calcium sensor that binds calcium and is involved in the early response to stress [110]; and NHX6 is involved in trafficking to the vacuole and exchanging the low-affinity electroneutral Na(+) or K(+)/H(+) [111]. In this study, the differential expression of these coexpressed mRNAs associated with energy and transport provided energy, nutrients, and metabolites for A. sinensis seedlings to obtain the capacity for vernalization and, meanwhile, to adapt to low temperatures.
Based on the above functional analysis of lncRNAs identified in this study, flowering pathways proposed in the previous literature [14][15][16] and a general model of stressresponsive regulation by regulatory lncRNAs [23], a model of vernalization-induced bolting and flowering by regulatory lncRNAs in A. sinensis is proposed ( Figure 8). Briefly, the vernalization of seedlings firstly triggers the differential expression of lncRNAs; then, the lncRNAs either act as a precursor of miRNAs or as a miRNA target mimic, which further binds their related targets; then, the binding of targets regulates the expression of their downstream mRNAs that are involved in various biological processes, including the temperature response, flowering pathways (i.e., epigenetic modification, flowering induction, carbohydrate metabolism, and hormone signaling), as well as energy and transport; finally, these biological processes promote the transition of the meristem from the vegetative phase to the bolting and flowering of A. sinensis. sion of lncRNAs; then, the lncRNAs either act as a precursor of miRNAs or as a miRNA target mimic, which further binds their related targets; then, the binding of targets regulates the expression of their downstream mRNAs that are involved in various biological processes, including the temperature response, flowering pathways (i.e., epigenetic modification, flowering induction, carbohydrate metabolism, and hormone signaling), as well as energy and transport; finally, these biological processes promote the transition of the meristem from the vegetative phase to the bolting and flowering of A. sinensis.

Conclusions
From the above observations, we found that the lncRNAs positively or negatively regulated the expression of their downstream mRNAs through epigenetic changes at the level of transcription and post-transcription for the flowering of A. sinensis during vernalization. This coexpressed mRNA analysis of lncRNAs focused on five pathways, namely (1) chromatin, DNA/RNA, and protein modification; (2) floral development; (3) temperature response; (4) carbohydrate metabolism; and (5) hormone signaling. While several candidate lncRNAs were identified, their causative roles require further investigations.

Conclusions
From the above observations, we found that the lncRNAs positively or negatively regulated the expression of their downstream mRNAs through epigenetic changes at the level of transcription and post-transcription for the flowering of A. sinensis during vernalization. This coexpressed mRNA analysis of lncRNAs focused on five pathways, namely (1) chromatin, DNA/RNA, and protein modification; (2) floral development; (3) temperature response; (4) carbohydrate metabolism; and (5) hormone signaling. While several candidate lncRNAs were identified, their causative roles require further investigations.