Transcriptome Analysis of Lycoris chinensis Bulbs Reveals Flowering in the Age-Mediated Pathway

Lycoris is a summer bulbous flower that commonly needs to go through a long period of vegetative growth for 3 to 5 years before flowering. Plant flowering is regulated by a complex genetic network. Compared with most perennial flowers, knowledge on the molecular mechanism responsible for floral transition in bulbous flowers is lacking, and only a few genes that regulate flowering have been identified with few reports on the floral transition in Lycoris. In this study, we identified many differentially expressed genes (DEGs) and transcription factors (TFs) by RNA-Seq in L. chinensis bulbs of different ages, including one- to four-year-old nonflowering bulbs and four-year-old flowering bulbs. Some DEGs were enriched in Gene Ontology (GO) terms between the three- and four-year-old bulbs, and there most genes were enriched in terms of metabolic process and catalytic activity. In the four-year old bulbs, most of the DEGs that may be involved in flowering were classified under the GO term biological process, which was a totally different result from the vegetative bulbs. Some DEGs between flowering and nonflowering bulbs were enriched in plant hormone signal transduction, including the hormones auxin, cytokinin, abscisic acid, and ethylene, but no DEGs were enriched in the gibberellin pathway. Auxin is the main endogenous phytohormone involved in bulb growth and development, but cytokinin, abscisic acid, and ethylene were shown to increase in flowering bulbs. In addition, energy-metabolism-related genes maintain a high expression level in large bulbs, and some positive regulators (SPL, COL, and AP1) and early flowering genes were also shown to be highly expressed in the meristems of flowering bulbs. It suggested that sugar molecules may be the energy source that regulates the signal transduction of flowering by connecting with phytohormone signaling in Lycoris. A total of 1911 TFs were identified and classified into 89 categories, where the top six families with the largest gene numbers were C2H2, NAC, AP2/ERF-ERF, C3H, MYB-related, and WRKY. Most DEGs were in the AP2/ERF-ERF family, and most of them were downregulated in 4-year-old flowering bulbs. A number of families were reported to be involved in plant flowering, including NAC, AP2/ERF, MYB, WRKY, bZIP, MADS, and NF-Y. These results can act as a genetic resource to aid in the explanation of the genetic mechanism responsible for the flowering of Lycoris and other bulbous flowers.


Introduction
The lifecycle of angiosperms includes four stages: embryo/juvenile, adulthood, flowering and reproduction, and senescence [1]. Flowering, an important stage in the lifecycle of plants, is not only beneficial for plant reproduction and species evolution, but is also a prerequisite for crop seed or fruit harvesting and is even an important trait of ornamental plants. Plant flowering is a process from flower formation induction to flower initiation, and flower organogenesis. Flower formation induction is the stage at which a plant transforms from vegetative to reproductive growth, and it determines the flowering time of the plant [2]. It is regulated by a complex gene network that integrates a variety of exogenous environmental factors (such as light and temperature) and plant endogenous factors (such as plant hormones and age) to ensure an accurate flowering time [3].
Previous genetic and molecular studies on the model plant Arabidopsis thaliana showed that there are five main pathways that regulate plant flowering: the photoperiod, vernalization, age, GA, and the autonomous pathway [4]. For example, the photoperiod pathway regulates FLOWERING LOCUS T (FT) and SUPPRESSOR OF OVEREXPRESSION OF CO1 (SOC1) through CONSTANS (CO) to upregulate expression and promote flowering. The interaction of CO and different proteins is very important for regulating the accurate transcription of downstream genes. Vernalization and autonomous pathways inhibit the activity of FLOWERING LOCUS C (FLC) and increase the expression of FT and SOC1, thereby promoting flowering [5]. In the aging pathway, miR156 and its target genes, SQUAMOSA PROMOTER BINDING PROTEIN-LIKE (SPL) and miR172, are key components of this pathway, regulating the transition from plant vegetative growth to reproductive development [6]. From the juvenile to the adult stage, the expression level of miR156 gradually decreases, and the expression of the negative regulatory target SPL increases, which promotes an increase in miR172 expression; miR172 downregulates the flower-forming inhibitor APETALA2 (AP2) and promotes the transcription of flowerforming integrins (FT, SOC1, and AGL24) and floral meristem recognition genes (FUL, LFY, AP1), finally inducing plant flowering [7]. Under the single or combined action of these pathways, the upstream genes precisely regulate the expression of the two key downstream flowering genes, FT and SOC1, at the transcription and translation levels, thereby affecting flowering.
In these regulatory pathways, some genes show relative conservation in the regulation of flowering in perennials, such as SOC1 and FT, which are among the key transcription factors that integrate signals from almost all genetic pathways. The key genetic regulators of floral transition have been well characterized in model species, and increasing numbers of related genes are being identified in crops [8,9]. However, the main regulatory pathways are varied in different species. For example, many annual herbaceous species germinate in the autumn and over winter and then flower in spring, having typical light and temperature responses during the flowering initiation stage [10]. In these taxa, low-temperature vernalization and long-day photoperiods are often required, which helps to relieve repressors of flowering and activate floral inducers [11]. However, in many perennials, such as most woody fruit trees [12,13] and bulb flowers [14,15], the acquisition of competence to flower is usually strongly delayed due to the high dependency of the miR156/SPL system [16], where the age pathway is the main factor that regulates flowering.
Bulbous flowers are rich in variety with strong stress resistance and are often used as garden flowers, cut flowers, and pot flowers. Similar to other perennial flowers, bulbous flowers also need to undergo a transition from vegetative growth to reproductive development before they can bloom. However, due to the bulb size and environmental factors, most bulbous flowers often take several years to flower [17]. Having a long vegetative growth period leads to higher production costs. Lycoris is a genus containing about 20 species of bulb-producing perennial plants from the family Amaryllidaceae, which are native to eastern and southern Asia and are mainly distributed in China and Japan [18]. They bloom from summer to autumn, and then the leaves grow in autumn, winter, and spring. The flowers have a special spidery shape and rich colors. In addition, their bulbs are abundant in Amaryllidaceae alkaloids. Among these, the component galantamine is used in medications to treat Alzheimer's-type dementia [19,20]. Thus, they have both high ornamental and medicinal value. However, it takes a long time for them to transition from seed to bloom, often as long as 3-5 years, so they are often propagated by the division of bulbs. It also takes 2-3 years of vegetative growth for the small asexually reproduced bulbs to complete the flower transformation process [21], which not only slows the breeding process of Lycoris, but also increases the production cost of commodity bulbs. Therefore, it is important to explore the factors that regulate the flowering transition of Lycoris.
Several environmental factors that may affect the flowering of Lycoris have been focused on in previous research, such as the low temperature [22], moderate shade, and water requirement [23]. These studies involved physiological and biochemical measurements, and provide some evidence that these endogenous and exogenous factors can promote or inhibit the flowering of Lycoris to a certain extent [22][23][24]. For example, earlier heating may lead to earlier flower bud development in spring [22], long days and drier environments increase reproductive growth, and moderate shade improves the yield and quality of cultivated L. radiata flowers [23]. However, it needs to be mentioned that all of these studies were based on bulbs with the ability to flower. The large Lycoris bulbs with flowering ability commonly experience a period of summer dormancy and then flower without leaves. The exogenous environmental factors may have some effect on the flowering time, but they are not the key genetic factors that determine whether the bulbs can flower or not. The main genetic factors that affect bulb development, apical meristems on flowering, and floral transformation are largely unknown. In the present study, we aimed to reveal the differentiation of gene expression and explore the key regulatory genes involved in the process of bulb development and flowering by high-throughput RNA sequencing (RNA-Seq) technology and to study the molecular regulation mechanism involved in flowering at the molecular level in Lycoris. The results provide abundant genetic resources that can be used in research on the flowering regulation of Lycoris and other bulbous flowers.

Plant Material and Sample Collection
One-to four-year-old L. chinensis bulbs were used in this study. Some four-year bulbs have the ability to flower, and some do not, and these were collected separately in this study. The bulbs were planted in the Nanjing Botanical Garden Mem. Sun Yat-Sen (Nanjing, China). The 1-4-year-old bulbs were selected (P1-P5) and cut along the middle of the bulb using a scalpel, and then the shoot apical meristem (SAM) was carefully cut off. Samples P1-P5 were 1-3-year-old bulbs, 4-year-old nonflowering bulbs and 4-year-old flowering bulbs. Fifteen meristem samples (3 replications for each sample) were collected and immediately frozen in liquid nitrogen and stored at −80 • C for RNA-Seq and qRT-PCR validation analyses.

RNA Extraction, Library Construction and Transcriptome Sequencing
Total RNA from 15 meristem samples was extracted using the Quick RNA Isolation Kit (Huayueyang Biotechnology Co., Ltd., Beijing, China) following the manufacturer's protocol. The purity and quality of RNA were detected by 1.0% (w/v) agarose-gel electrophoresis, and the quality and quantity were further measured using the Agilent 2100 bioanalyzer. Then, the cDNA libraries for Illumina sequencing were constructed by mRNA enrichment, double-stranded DNA synthesis, sequencing adapter modification and fragment size selection, PCR amplification, and library detection. The qualified libraries were sequenced on the Illumina HiSeq2500 platform at Novogene Bioinformatics Technology Co., Ltd. (https://www.novogene.com/, accessed on 20 January 2018, Tianjin, China). Raw reads were deposited in the NCBI database (https://www.ncbi.nlm.nih.gov/, accessed on 8 June 2022) under BioProject number PRJNA847051 (SRA accession: SRR19578267-SRR19578281).

Raw Data Processing, Transcriptome Assembly, and Gene Functional Annotation
In order to fit the information analysis, the raw reads were processed by removing the adaptor sequences. Reads with a ratio of more than 0.1% N (N means that the base information cannot be determined) and the low-quality reads were removed. Then, the clean reads were assembled de novo by Trinity software (r20140413p1, http: //trinityrnaseq.github.io, accessed on 30 January 2018, Broad Institute and Hebrew University of Jerusalem) [25]. The longest transcript sequence of each gene was taken as the unigene for subsequent analysis. For the functional annotation of all the unigenes, BLASTx (E-value ≤ 10 -5 ) was carried out on seven major public databases, including the Nr (NCBI nonredundant protein sequences), Nt (NCBI nucleotide sequences), Pfam (Protein family, http://pfam.sanger.ac.uk/, accessed on 30 January 2018), KOG (euKaryotic Ortholog Groups), Swiss-Prot (A manually annotated and reviewed protein sequence database, http://www.ebi.ac.uk/uniprot/, accessed on 30 January 2018), KEGG (Kyoto Encyclopedia of Genes and Genomes, http://www.genome.jp/kegg/, accessed on 30 January 2018), and GO (Gene Ontology, http://www.geneontology.org/, accessed on 30 January 2018). The GO and KEGG annotated genes were classified. Then, the GO annotated genes were classified into three major categories (BP, biological process; CC, cellular component; MF, molecular function). KOG is divided into 26 groups, and the KO annotated genes were classified according to the KEGG metabolic pathways they participate in. The iTAK program and database were used to identify the transcription factors [26].

Differentially Expressed Gene (DEG) Analysis
The clean reads from each sample were mapped to the reference sequence assembled by Trinity using RSEM software [27] with the bowtie2 default parameter mismatch 0. According to the statistics of the bowtie comparison, the readcount number of each sample mapped to each gene was obtained, and further FPKM values were conversed and the gene expression level was analyzed. The readcount was input for the differentially expressed gene (DEG) analysis. In this study, three replications of each sample were used. We used DESeq [28] for the DEGs analysis with a screening threshold of padj < 0.05. When the log 2 FoldChange of samples was more than 1.5, it was screened for DEG analysis. In order to better understand the function of DEGs, we performed GO and KEGG enrichment analyses for each comparison. GOseq [29] was used for the GO enrichment analysis, and the number of significantly enriched genes in each GO term was displayed in the form of a histogram. The KEGG enrichment analysis was performed by KOBAS (2.0) with the parameters -fdr to BH and FDR ≤ 0.05. The 20 most enriched pathway entries are displayed in a scatter plot. If there are less than 20 enriched pathway entries, all pathway entries were displayed. In addition, the pathway maps are also shown for the enriched differential genes. The clustered heat maps were drawn using Tbtools software [30] with a log2 scale.

qRT-PCR Validation
To validate the accuracy and reliability of the transcriptome results, 20 DEGs involved in flowering were randomly selected for the quantitative real-time polymerase chain reaction (qRT-PCR), and EXP1 was selected as the internal standard control [31]. The genespecific primers for qRT-PCR (Table S1) were designed using Primer3 Release 2.3.4 [32] (pp. 365-386) and synthesized by General Biosystems Co., Ltd. (Chuzhou, China). qRT-PCR was conducted on a StepOnePlus™ real-time PCR system (Applied Biosystems, Foster City, CA, USA). The first cDNA strand was synthesized with the PrimeScript TM reagent kit (Takara, Dalian, China). The mRNA expressions were quantified using the SYBR Premix Ex Taq II kit (TaKaRa, Dalian, China). Three biological and technical replicates for each gene were employed. The mean ± standard deviation of three biological replicates for each sample was calculated. The relative RNA expression levels were analyzed using the 2 −∆∆CT method [33]. The STDEVA function in Excel was used to assess the standard errors of deviation.

Overview of Transcriptome Sequencing Data Analysis
In order to detect the key factors that regulate flowering in bulbs of different ages in L. chinensis, five samples (three biological replicates for each sample) were selected for the construction of 15 cDNA libraries and high-throughput RNA-Seq. They were 1-3-year-old bulbs that cannot flower, and 4-year-old bulbs, both nonflowering and flowering (Figure 1a-c, P1-P5). A summary of the RNA-Seq results is shown in Table 1. In each library, around 50 million clean reads were generated and about 60% of the reads were mapped in reference sequences. For each library, 7.0 Gb of data were obtained. The Q30 percentages for all libraries exceeded 91.5%, and the GC (Guanine and Cytosine) contents ranged from 45.65% to 47.05%. A total of 200,902 unigenes were identified (Table S2). To identify the differentially expressed genes affecting bulb development and flowering regulation, we compared the global transcriptomic profiles of five samples (P1-P5). The Venn diagram in Figure 1d shows the number of DEGs obtained from four group comparisons, where 101 genes were commonly and differentially expressed during the five stages. During the total vegetative growth stage (P1-P3), the number of DEGs was about 4000, of which upregulated and downregulated genes accounted for about half of these ( Figure 1e,f, Table S3). During the transition from vegetative growth to reproductive growth (P3 to P4 and P4 to P5), the number of DEGs increased significantly, reaching more than 5500 (Figure 1g,h, Table S3). The volcano plots for all combinations are shown in Figure 1e-n.

Gene Function Annotation
In order to understand the functions of the unigenes, gene function annotations were conducted in seven major databases, including Nr, Nt, Pfam, KOG/COG, Swiss-Prot, KEGG, and GO. A total of 29.62% of them were annotated in at least one database, and the percentage of annotated gene functions in NR, NT, GO, PFAM, and Swiss-Prot was 22.71, 10.34, 17.21, 17.05, and 15.04, respectively ( Table 2). After annotation of GO, KOG, and KEGG, the unigenes were categorized into different groups. For GO annotation, three categories were classified: BP (biological process), CC (cellular component), and MF (molecular function). Most unigenes were annotated with the terms 'cellular process', 'metabolic process', and 'binding' (Figure 2a). A total of 25 classifications were annotated with KOG functions. The genes related to 'posttranslational modification', 'protein turnover', and 'chaperones' were the most abundant, accounting for more than 15% of the total number. This was followed by the classifications of 'general function prediction only' and 'translation, ribosomal structure, and biogenesis', which accounted for more than 10% of the total genes. The least common categories were 'cell motility', 'extracellular structures', and 'defense mechanisms' (Figure 2b). After the KOG annotation, the genes were classified by KEGG according to the metabolic pathways they participate in. They genes were found to be involved in a total of 19 pathways. Of them, the 'translation' pathway enriched the largest number of genes, followed by the 'folding, sorting, and degradation' and 'carbohydrate metabolites' pathways. The number of genes involved in these three classifications was more than 1300, which is significantly higher than that in other pathways (Figure 2c).

GO and KEGG Enrichment of Differentially Expressed Genes (DEGs)
The GO enrichment of DEGs identified the biological functions significantly related to these genes. In the pairwise comparison enrichment analysis, enriched GO terms were shown in bulbs aged more than three years, including comparisons for P3 vs. P4 and P4 vs. P5 (Figure 3c,d). These pairs had significantly more genes than P1 vs. P2 and P2 vs. P3, where the number of DEGs in most GO terms was less than 25, and most genes were enriched in the oxidation-reduction process (Figure 3a,b). During the transition from vegetative growth to reproductive growth, most genes were classified by the GO term molecular function and two processes, metabolic processes (biological process), and catalytic activity (molecular function), which were enriched in more than 700 genes (Figure 3c). In four-year-old bulbs, most genes that determine whether to flower were classified under the GO term biological process, and these terms were totally different from the DEGs identified in the other three comparisons (Figure 3d).

Key Genes Involved in Bulb Growth and Flowering in L. chinensis
Because of the genetic conservation of plant flowering, many genes have been shown to be involved in different flower-forming pathways and have conserved regulatory functions in model plants. In this study, we classified the key up-and downregulated genes that may be involved in flowering between P4 and P5 samples. In flowering bulbs (P5), there were more upregulated genes than in four-year-old nonflowering bulbs (P4), for example, squamosa promoter-binding-like protein 5 (SPL5), agamous-like MADS-box protein AGL14/AGL61, zinc finger protein CONSTANS-LIKE 15-like, zinc finger CCCH domaincontaining proteins, AP2/ERF, energy-metabolism-related genes, histones, hormone-related genes (cytokinin response regulator protein, abscisic acid receptor, indole-3-acetic acidamido synthetase GH3.1, gibberellin-regulated protein), and some transcription factors. The downregulated genes in P5 were the terminal flower 1-like protein, UPSTREAM OF FLClike, AP2-11, auxin-responsive and induced protein, ethylene insensitive 3, E3 ubiquitinprotein, and F-box protein and transcription factors (IBH1-like and NAC) ( Table 3). Based on 469 DEGs potentially involved in flowering, principal component analysis (PCA) was performed to explore the relationships between gene expression and samples ( Figure S1). Results from the PCA analysis were able to distinguish samples of different ages together with flowering and nonflowering bulbs, indicating that the biological replicates of sam-ples were good and the genes involved in bulb development and flowering had specific expression patterns in each sample.
Based on the genes identified as being involved in flowering regulation, hormone signaling is considered the main mediating factor in the processes of vegetative to reproductive transformation in plants. In the KEGG enrichment of DEGs, some genes were shown to be enriched in the main pathways of plant hormone signal transduction, including auxin, cytokinine, abscisic acid, and ethylene, but no DEGs were enriched in the gibberellin pathway ( Figure 5). In the auxin signal transduction pathway, two unigenes (AUU/IAA), annotated as the auxin-responsive protein IAA and the auxin-induced protein, were downregulated in the flowering bulb. Four unigenes encoding the auxin-response protein SAUR and two unigenes encoding indole-3-acetic acid-amido synthetase GH3 were upregulated in the flowering bulbs. They work together to regulate cell enlargement and plant growth and further affect Lycoris flowering (Figure 5a). In the cytokinin signal transduction pathway, one unigene encoding a type-b response regulator (B-ARR) was downregulated in the flowering bulbs, showing the same pattern as unigenes related to AUU/IAA biosynthesis (Figure 5b). In the abscisic acid signal transduction pathway, four abscisic acid receptors PYL4 were upregulated and shown to inhibit protein phosphatase 2C synthesis, but two unigenes of protein phosphatase 2C were also upregulated, showing the negative regulatory function of stomatal closure and seed dormancy. This indicates that abscisic acid may play a negative role in the Lycoris floral transition (Figure 5c). In the ethylene biosynthesis pathway, one ethylene-insensitive protein and two ethylene-responsive transcription factors showed higher expression levels in the flowering bulbs, suggesting that they play positive roles in Lycoris flowering (Figure 5d).  (a-d) The auxin, cytokinin, abscisic acid, and ethylene signal transduction pathways. The borders of KO nodes of upregulated genes are marked in red, and the downregulated genes are marked in green. The expression level of the related unigenes from P3 to P5 are shown in the heatmap beside the nodes, the color from blue to red means that the gene expression level is from low to high. In addition, transcription factor (TF) genes play various roles in the regulation of plant growth and development and the defense response to adversity. Multiple TFs families were differentially expressed in different samples. Here, we identified 1911 TFs using the iTAK program through the RNA-Seq of L. chinensis bulbs. They were classified into 89 categories, and 48 families with more than 10 unigenes are shown in Figure 6a. A number of families were reported to be involved in plant flowering, including NAC, AP2/ERF, MYB, WRKY, bZIP, MADS, and NF-Y. In these families, the top six with the largest numbers of members (more than 80 in each family) were C2H2, NAC, AP2/ERF-ERF, C3H and MYB-related, and WRKY. In the AP2/ERF-ERF family, there were the greatest number of differentially expressed genes in the samples (Figure 6d). The differentially expressed genes in four families other than C3H (C2H2, NAC, MYB-related, and WRKY) are shown in heatmaps (Figure 6b,c,e,f).

qRT-PCR Validation of the RNA-Seq Results
In order to explore the main genetic factors related to bulb age and flowering regulation, we summarized the key genes involved in regulating bulb development and the transition from vegetative growth to reproductive development and verified the accuracy of the transcriptome by qRT-PCR. A total of 20 key genes mediating the flowering time were selected for qRT-PCR validation, including FT, CONSTANS-LIKE 1, CONSTANS-LIKE 15, APETALA1, APETALA2, AGL72, NF-YB3, and some genes of the SPL family. The expression patterns of the candidate unigenes indicated a high degree of

qRT-PCR Validation of the RNA-Seq Results
In order to explore the main genetic factors related to bulb age and flowering regulation, we summarized the key genes involved in regulating bulb development and the transition from vegetative growth to reproductive development and verified the accuracy of the transcriptome by qRT-PCR. A total of 20 key genes mediating the flowering time were selected for qRT-PCR validation, including FT, CONSTANS-LIKE 1, CONSTANS-LIKE 15, APETALA1, APETALA2, AGL72, NF-YB3, and some genes of the SPL family. The expression patterns of the candidate unigenes indicated a high degree of consistency between qRT-PCR and RNA-Seq, indicating the reliability of the RNA-Seq data (Figure 7).

Discussion
It is well known that the flowering transition of plants is a complex morphological and physiological change that occurs in response to endogenous factors (hormones, age, autonomous, and sugar signals) and exogenous environmental factors (light and temperature). Plants integrate these responses to act on downstream floral integrators to regulate flowering [34,35]. For nonmodel plants and crops without a reference genome, a large amount of RNA-Seq data facilitates gene mining to regulate the plant stress responses, growth, and development [36] (pp. 175-185) [37,38], which also has been widely used to explore the genetic regulation of flowering in many plants, such as bamboo orchid flowering [39] and the flowering transition of saffron [40] and roses [41]. Compared with many perennial flowers, data on gene identification and molecular mechanisms that regulate flowering in bulbous flowers are insufficient, and only a few genes have been identified. For example, FT and LFY regulate the flowering of Narcissus tazetta [42], and COL and FT regulate the flowering of Lilium [43,44]. The reproductive initiation gene TGSQA (a homolog of Arabidopsis APETALA1) acts as a flowering activator to regulate flowering in Tulipa gesneriana [45]. In bulbous flowers, in the process of transitioning from vegetative growth to reproductive development and acquiring the flowering ability, the genes and substances in the meristem at the base of the bulb change as the bulb size increases and flower when the environment is suitable [15]. Therefore, bulb size and age are considered key factors affecting flowering in bulbous flowers. In this study, a total of 200,902 unigenes were assembled from 1-4-year-old vegetative and reproductive L. chinensis bulbs, and many transcription factors and key genes that may regulate plant flowering were identified. The result provides genetic information for the flowering regulation in L. chinensis.
In bulbs of different ages, a large number of DEGs were identified. These included many key genes and TFs that may be involved in bulb development and flowering, such as the key up-and downregulated genes between four-year-old nonflowering and flowering bulbs (Table 3). According to the pattern of these DEGs, the overall conclusion is that metabolic processes and energy accumulation are enhanced with increasing age and enlarged bulb size, and auxin is the main endogenous phytohormone involved in bulb growth and development during this process. Once the age and size of bulbs reach the bloom standard, many key genes and factors perform functions to regulate flowering, including positive regulators (SPL, COL, and AP1, high expression in flowering bulbs), negative regulators (FLC and AP2, terminal flower gene, low expression in flowering bulbs), early flowering genes, and histones. Meanwhile, energy metabolism and hormones also contribute to the floral transition and flowering. The energy-metabolism-related genes maintain high expression levels, and the endogenous phytohormones, such as cytokinin, abscisic acid, and ethylene also increase in flowering bulbs. In contrast, the genes related to auxin synthesis are downregulated (Figure 8). It is suggested that sugar molecules may be the energy sources that regulate flowering signal transduction by connecting with phytohormone signaling in Lycoris.
The enrichment analysis of DEGs found that the DEGs were enriched in different KEGG pathways only in the three-year-old and four-year-old bulbs, suggesting that these pathways may be the key pathways involved in the floral transition in L. chinensis. From three-year-old to four-year-old nonflowering bulbs, two pathways, 'metabolic process' and 'catalytic activity', enriched the most DEGs. In the P3 vs. P4 and P4 vs. P5 comparisons, two pathways 'starch and sucrose metabolism' and 'anthocyanin biosynthesis' enriched DEGs, but the DEGs were enriched in the plant hormone signal transduction pathway in the P4 vs. P5 group. Sugar molecules function as energy sources and osmotic regulators and have been recognized as important regulators of floral signal transduction that connect with phytohormone signaling in many plants [46][47][48]. Proteomic studies of the flowering and dormancy of L. radiata revealed a large number of differentially expressed proteins involved in the energy and sugar metabolic processes by GO and pathway enrichment [49]. Therefore, our results indicate that there may be active starch and sucrose metabolism, anthocyanin biosynthesis, and catalytic activity in the meristem with the enlargement of the bulbs, and the plant hormone signal is the main factor that induces the floral transition in four-year-old L. chinensis bulbs. In the metabolism or response to phytohormones, the gibberellic acid (GA) pathway did not enrich DEGs between two samples of four-year-old bulbs. However, in the auxin, cytokinin, abscisic acid, and ethylene pathways, there were some differentially expressed genes that may affect the flowering regulation of L. chinensis; for example, the genes regulating auxin synthesis had lower expression levels in flowering bulbs, but the genes involved in cytokinin, abscisic acid, and ethylene were commonly upregulated for flowering (Table 3, Figure 5). In complex floral transition regulatory networks, the five regulatory pathways are not independent but often interact with each other to achieve their functions. Many key genes, transcription factors, and miRNAs have been reported to regulate flowering in different pathways, including integrin genes (FT, SOC1, COL, LEAFY, AP1, AP2), several transcription factor families, and miRNA156/172. For example, COL is regulated by many factors related to the meristem expression pattern, including temperature, light, transcription factors, and phytohormones (auxin, gibberellin, and abscisic acid) [41,50]. In the meristem of L. chinensis bulbs, 89 TFs families were identified, covering almost all families involved in plant growth and development. Most of the families were found to control the floral transition, including C2H2 [51], NAC [52], MYB [53], WRKY [54,55], AP2/ERF [56,57], NF-Y [58], and MADS-domain transcription factor [59]. Transcriptome and qRT-PCR analyses confirmed that different transcription factors have different expression patterns during bulb development and floral transition. For example, the MYB-related MYB1R1, NAC48, WRKY38, ERF12, and CRF1 have the highest expression levels in the four-year-old flowering bulbs, suggesting positive regulation of floral transition. In contrast, MYBAS2, NAC7, WRKY70, and ERF110 showed the lowest expression levels and may negatively regulate floral transition in L. chinensis. We also verified the expression of other core genes by qRT-PCR. SPL14 and SPL16 were shown to have the highest expression levels in the flowering bulbs. In the age pathway of flowering, SPL was shown to positively regulate the expression of miR172, and miR172 negatively regulates AP2 (inhibitor of FT), finally leading to the regulation of plant flowering [4]. In this study, AP2 was downregulated in the flowering bulbs, which is consistent with the mechanistic models of plant flowering by age pathway. Therefore, this study showed that the age is the main factor that regulates the flowering in L. chinensis, but the floral transition is a complex process involving physiological and biochemical changes, and it may be regulated by the sugar metabolism and hormone signal transduction pathways. Here, we provided abundant candidate differentially expressed TFs and genes associated with the bulb development and floral transition. Our data provide a genetic foundation for further research on the flowering regulatory mechanism of Lycoris.

Conclusions
In this study, we identified a large number of differentially expressed genes and transcription factors by RNA-Seq in L. chinensis bulbs of different ages. These DEGs were enriched in different GO and KEGG pathways in nonflowering and flowering bulbs, suggesting that they may regulate flowering by pathways including age, hormones, and vernalization in Lycoris. Starch and sucrose metabolism, anthocyanin biosynthesis, and catalytic activity are enhanced with increasing age and enlarged bulb size, and auxin is the main endogenous phytohormone involved in bulb growth and development. Once the bulb size is sufficient for flowering, the energy-metabolism-related genes maintain the high expression level, and the endogenous phytohormones (cytokinin, abscisic acid, and ethylene) also increase in flowering bulbs, causing some positive regulators (SPL, COL, and AP1) and early flowering genes in the Lycoris meristem to have high expression levels. These results suggest that sugar molecules may be the energy sources that regulate the floral signal transduction by connecting with phytohormone signaling in Lycoris. In these regulatory pathways, some key genes and TFs shown to contribute to plant flowering have been identified, but they have rarely been revealed in relation to Lycoris flowering regulation. Thus, our results provide an abundant genetic resource that can be used for further studies on flowering transition and time.
Supplementary Materials: The following supporting information can be downloaded at: https: //www.mdpi.com/article/10.3390/biom12070899/s1, Figure S1: Principal component analysis of 469 DEGs in five samples of Lycoris chinensis. Each dot represents a sample, and colors code for stages P1-P5; Table S1: qRT-PCR premiers used in the study; Table S2: All unigenes identified in the study; Table S3: Differentially expressed genes in four comparison groups.
Author Contributions: F.Z. and Z.W. conceived and designed the experiments; G.C. and X.S. prepared the samples and performed the experiments. F.Z. and N.W. analyzed the data. F.Z. prepared the draft; Z.W. and F.Z. reviewed and edited the manuscript. All authors have read and agreed to the published version of the manuscript.