Transcriptomic analyses reveal molecular mechanisms underlying regulation of oral attributes in Lonicera japonica Thunb.

Background: Lonicera Japonica Thunb. is a perennial, semi-evergreen and twining vine in the family of Caprifoliaceae, which is widely cultivated in Asia. Thus far, L. japonica is often used to treat some human diseases including COVID-19, H1N1 influenza and hand-foot-and-mouth diseases, however, the regulatory mechanism of intrinsic physiological processes during different floral developmental stages of L. japonica remain largely unknown. Results The complete transcriptome of L. japonica was de novo-assembled and annotated, generating a total of 195850 unigenes, of which 84657 could be functionally annotated. 70 candidate genes involved in flowering transition were identified and the flowering regulatory network of five pathways was constructed in L. japonica . The mRNA transcripts of AGL24 and SOC1 exhibited a downward trend during flowering transition and followed by a gradual increase during the flower development. The transcripts of AP1 was only detected during the floral development, whereas the transcript level of FLC was high during the vegetative stages. The expression profiles of AGL24 , SOC1 , AP1 and FLC genes indicate that these key integrators might play the essential and evolutionarily conserved roles in control of flowering switch across the plant kingdom. We also identified 54 L. japonica genes encoding enzymes involved in terpenoid biosynthesis pathway. Most highly expressed genes centered on the MEP pathway, suggesting that this plastid pathway might represent the major pathway for terpenoid biosynthesis in L. japonica . In addition, 33 and 31 key genes encoding enzymes involved in the carotenogenesis and anthocyanin biosynthesis pathway were identified, respectively. PSY transcripts gradually increased during the flower development, supporting its role as the first rate-limiting enzyme in carotenoid skeleton production. The expression level of most anthocyanin biosynthetic genes was dramatically decreased during the flower developmental stages, consistent with the decline in the contents of anthocyanin. Conclusion These results identified a large number of potential key regulators controlling flowering time, flower color and floral scent formation in L. japonica , which improves our understanding of the molecular mechanisms underlying the flower traits and flower metabolism, as well as sets the groundwork for quality improvement and molecular breeding of L. japonica .


Introduction
Lonicera japonica Thunb. is a perennial, semi-evergreen and twining vine in the family of Caprifoliaceae. Since it has double-tongued owers that open white and fade to yellow, this ornamental plant is also called "Jin Yin Hua" in Chinese (literally "gold and silver ower"). L. japonica is widely cultivated in Asian countries such as China, Japan, and Korea. However, due to its strong viability, L. japonica is believed to be invasive to the ecology of some western counties [1]. L. japonica has been used as traditional Chinese Medicine for over thousands of years for its anti-bacterial, anti-viral, anti-endotoxin, anti-in ammatory and anti-pyretic properties and has been recorded in the "Shen Nong Ben Cao Jing" and "Ben Cao Gang Mu" (Compendium of Materia Medica) as early as the 17th century [2]. The commercial value of L. japonica in the herbal medicine trading market has increased rapidly in recent years. More than 30% of current traditional Chinese medicine prescriptions have been reported to contain extracts from different parts of L. japonica [3]. Thus far, L. japonica is often used to treat some human diseases such as acute respiratory syndromes, H1N1 in uenza, hand-foot-and-mouth diseases, pancreatic cancer [4,5]. In the current, a pandemic of coronavirus disease 2019 (COVID-19) has spread globally. COVID-19 is caused by severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2). Recently, emerging data indicate that L. japonica-based substances may be effective in the prevention and treatment of SARS-CoV-2 infection [6][7][8][9]. Chlorogenic acid (CGA) and luteoloside are the main active components used for evaluating the quality of L. japonica. The contents of CGA and luteoloside are highest accumulation in alabastrum, while lower in other developmental stages [10,11]. Although some transcriptomic studies have been conducted to analyze the biosynthesis of pharmaceutically active ingredients during the ower development [10][11][12], the regulatory mechanism of intrinsic physiological processes during different oral developmental stages of L. japonica remain largely unknown.
Flowering is an intricate biological process during the angiosperm life cycle. In recent years, the molecular framework for integrating different signaling pathways relevance to oral traits including owering time, oral color and oral aroma have been studied in plants such as rice, gourds, potato and sorghum, poplar, citrus, and especially in model plant Arabidopsis thaliana [13,14]. Flowering time are controlled by ve major pathways, including the photoperiod pathway, the autonomous pathway, the vernalization pathway, the gibberellin pathway and the circadian clock pathway [15,16]. More than 300 owering genes have been established a comprehensive network to respond to endogenous or external environmental cues in Arabidopsis [15]. These genes activate or repress oral transformation through a small number of ower integrator genes, such as FLOWERING LOCUS T (FT), FLOWERING LOCUS C (FLC), SUPPRESSOR OF OVEREXPRESSION OF CONSTANS 1 (SOC1), AGAMOUS-LIKE 24 (AGL24), APETALA1 (AP1) genes [17,18]. Although the molecular aspect of owering time control in Arabidopsis has been widely studied, the genes associated with L. japonica owering time are almost unknown. L. japonica also has commercial value for its fragrant volatiles that are widely used in manufacture of perfumes, scented oils, bathing products, and other household products. Terpenes with relative low molecular weight (such as monoterpenes and sesquiterpenes) account for the largest proportion of volatile organic compounds in L. japonica owers [19]. The biosynthesis of plant volatile terpenes involves three steps: (1) production of the universal C5 isoprenoid precursors, isopentenyl diphosphate (IPP) and dimethylallyl diphosphate (DMAPP), (2) formation of geranyl diphosphate (GPP, C10) and farnesyl diphosphate (FPP, C15) by condensation of IPP with DMAPP, and (3) conversion of GPP and FPP substrates into monoterpenes and sesquiterpenes respectively. In the cytosol, IPP and DMAPP are derived from acetyl-CoA through the mevalonic acid (MVA) pathway and catalyzed to FPP by farnesyl diphosphate synthase (FPPS). Then FPP is converted to the various sesquiterpenes by the cytosolic/mitochondrial terpene synthases (TPSs) [20]. In the plastids, IPP and DMAPP are synthesized from pyruvate and glyceraldehyde-3-phosphate (G3P) by the 2-C-methyl-D-erythritol 4-phosphate (MEP) pathway and converted into GPP by GPP synthase (GPPS). The plastidial TPSs catalyze GPP into the various monoterpenes [20]. The monoterpenes and sesquiterpenes may then be utilized to synthesize various terpene derivatives by cytochrome P450 (CYP450) enzymes, dehydrogenases and reductases, dehydrogenases and reductases [21]. For example, CYP76C1 and CYP76C3 catalyze linalool to 8hydroxylinalool in Arabidopsis [22].
Flower color is directly determined by distinct colorful pigments. During L. japonica development, the ower color transits from the initial pale green to white and nally yellow. The developmental color change is mainly due to the increased contents of carotenoid as well as the decreased contents of chlorophyll and anthocyanin [23,24]. The main pathways of anthocyanin and carotenoid biosynthesis have been well understood in model plants, and it is known to be fairly well conserved among the plants [25,26]. Phytoene synthase (PSY) is the rate-limiting enzyme in carotenoid skeleton production. Phytoene then undergoes transformation catalyzed by phytoene desaturase (PDS), 15-cis-ζ-carotene isomerase (Z-ISO), ζ-carotene desaturase (ZDS), and carotenoid isomerase (CRTISO) to produce lycopene. Lycopene is transformed to δ-carotene and γ-carotene by lycopene ε-cyclase (LYCE) and lycopene β-cyclase (LYCB), respectively. Anthocyanin is one kind of water soluble natural pigment which belongs to avonoid family.
The previous study of L. japonica mainly focused on the biosynthesis of pharmaceutically active ingredients [10-12, 27, 28], yet relevant researches for exploring the physiological and developmental regulation during L. japonica ower growth remain scarce. Here, we used RNA-seq technology to analyze the transcriptome of L. japonica leaf, ower bud, white ower and yellow ower tissues. Our detailed analysis identi ed a large number of potential key regulators controlling owering time, ower color and oral scent formation. The expression patterns of the genes related to oral traits were analyzed and their dynamic regulatory networks were constructed in L. japonica. The results of this study improve our understanding of the molecular mechanisms underlying the ower traits and ower metabolism in L. japonica, as well as set the groundwork for quality improvement and molecular breeding of this species.

Plant materials
The fresh leaves, buds, white and yellow owers were collected 5-year-old L. japonica plants from Fengqiu cultivation base (35°05'N and 114°25'E), Henan province, China. Samples with similar morphology (i.e., young leaves ~2 cm width, buds not yet bloomed into a full-size ower, white owers and yellow owers) were manually collected separately from three plants excluding subtending bract and bracteole. Samples of 3 plants at the same stage were combined and regarded as one biological replicate that representing each stage, and three independent replicates were performed. These samples were immediately frozen in liquid nitrogen following collection and stored at -80°C.
RNA extraction, cDNA library construction and transcriptome sequencing Each sample was grounded to powder in liquid nitrogen using a mortar and pestle. Total RNA was extracted by using RNeasy Plant Mini kit (Qiagen, Hilden, Germany). The quantity and quality were determined using Agilent 2100 Bioanalyzer (Agilent Technologies, Santa Clara, CA, USA). The poly (A) mRNA was isolated from total RNA samples with Magnetic Oligo (dT) Beads (Illumina, San Diego, CA, USA) and mRNAs were added with fragmentation buffer to shear mRNA into short fragments (200 ~nt).
First strand cDNA was synthesized using random hexamer primer and M-MuLV Reverse Transcriptase (RNase, H-), then second strand cDNA was synthesized using DNA polymerase I and RNase H. The AMPure XP system (Beckman Coulter, Beverly, USA) was used to select cDNA fragments preferentially 150-200 bp in length. Products were puri ed using AMPure XP system. The twelve libraries were sequenced using Illumina Hiseq X Ten platform.

De novo assembly and functional annotation
Clean, high quality reads were obtained by removing the adaptor sequences, reads with ambiguous "N" bases and reads with more than 10% Q<20 bases. At the same time, Q20, Q30 and GC-content was of the clean data were calculated. De novo assembly of the transcriptome was performed Trinity software [29] with default parameters and no reference sequences. Brie y, reads were assembled into contigs using an optimized k-mer length (25-mer) by the Inchworm program at the rst step. The minimally overlapping contigs were clustered into connected components by the Chrysalis program, and then the transcripts were constructed by the Butter y program. Finally, the transcripts were further clustered based on nucleotide sequence similarity, and the longest transcripts in each cluster was regarded as a unigene. For detection of transcription pro les in different tissues, clean reads from each library were initially assembled separately. To obtain a uniform transcriptome reference across the samples, all clean reads from twelve libraries were pooled together and de novo assembled to generate unigenes for assembly evaluation, gene annotation, and expression analysis. To identify putative functions of L. japonica unigenes, the unigenes were annotated with the NCBI non-redundant databases (NR), Pfam, Swiss-Prot, Gene Ontology (GO), Clusters of Orthologous Group (COG) and Kyoto Encyclopedia of Genes and Genomes (KEGG) databases using BLASTX searches (E-value ≤1.0e -5 ). The Blast2GO software package was used to compare and determine the unigene GO annotations [30]. Finally, WEGO2.0 was used to obtain the GO functional classi cations of all annotated unigenes [31]. Based on KEGG mapping, unigenes were assigned to multiple pathways using BLASTx, thereby retrieving KEGG Orthology (KO) information.

Analysis of differentially expressed genes (DEGs)
To calculate the amount of gene expression of each unigene, we used the abundance of reads derived from each locus to estimate gene expression and calculate transcripts per kilobase million (TPM) values with the program RSEM (RNA-Seq by Expectation Maximization) [32]. By pairwise comparisons of the different libraries, the DEGs were performed with EdgeR package based on the read counts for each gene at different libraries [33]. In addition, the false discovery rate (FDR) control method was used to identify the threshold of the P-value in the signi cance tests. Signi cance of DEGs was determined based on FDR<0.01 and fold change ≥2 (log2 Ratio≥1).

Validation by RT-qPCR
Twelve selected DEGs involved in ower timing, oral scent and ower color were determined by RT-qPCR. Total RNA was extracted using the same method as described above, and then the rst strand cDNA was synthesized using the Superscript III First Strand cDNA Synthesis System (Invitrogen, USA) following the manufacture's protocol. The primers were designed with the Primer premier 5.0 software and synthesized by Sangon Biotech Co., Ltd. (Shanghai, China). LjActin was used as an internal standard. The relative expression levels of genes were calculated using the 2−ΔΔCt method. The qPCR was conducted using the LightCycler® 480 (Roche, Germany) with the SYBR Green real-time PCR Master Mix. Real-time quanti cation was performed in a reaction mix containing 2 ul cDNA, 0.2uM gene-speci c forward and reverse primers, 1*SYBER Green master mix in a nal 20 μl reaction. qPCR cycling conditions were for 40 cycles, of 15s at 95℃, 30s at the annealing temperature at 60 ℃, 30s at the extension at 72 ℃ with 5min at 72 ℃ for nal extension.

RNA-seq and de novo assembly yielded 195,850 unigenes
To deeply understand the molecular mechanism of ower development in L. japonica, twelve libraries were constructed and sequenced by using an Illumina HiSeq X Ten platform with 150 paired-end reads. After cleaning and quality checking, about 107.40 Gb total clean reads were obtained and each library produced no less than 6 (Table S1). Subsequently, de novo assembly generated 408409 transcripts and 195850 unigenes, with the N50 (covering 50% of all the nucleotide sequences of the largest unigene length) of 2099 nt and 1310 nt, respectively. In total, 91643 unigenes (46.79%) were between 201 and 400 nt; 35435 unigenes (27.20%) were between 401 and 800 nt; 18007 unigenes (9.19%) were between 801 and 1200 nt; 9403 unigenes (8.02%) were between 1201 and 2000 nt; 17231 unigenes (8.80%) were longer than 2000 nt (Table 1).

Transcriptome pro les of L. japonica
In order to understand the relationship between the transcriptome database, we used read count data to perform unsupervised principal component analysis (PCA) (Fig. S1a). Along the PC1 axis, we observed three clusters, including opened ower tissues (white owers (WF) and yellow owers (YF)) within the rst group, oral buds (BD) and leaves (LF) as the other two groups. Along the PC2 axis, the opened owers were further separated into two groups. Correlation analysis based on Pearson's distance matrix using the entire transcriptome dataset were shown in Fig. S1b. The PCA analysis and correlation matrix showed a good correlation between the replicate sets for each of the four tissue samples. These results suggested that WF and YF showed the most similarity of transcriptomes, whereas LF showed the greatest difference to other tissues. These results also indicated that the presence of signature unigenes associated with and speci c to each tissue.

Differential expressional genes (DEGs) analysis in different tissues
Following transcriptome assembly and annotation, DEGs analysis was carried out with a criteria of fold change ≥2 and FDR < 0.05 between different tissues. Substantial transcriptional differences were seen in the pairwise comparisons between different tissues (Fig. 3a). The number of DEGs was high in the vegetative tissue when compared to the reproductive tissues. In the LF, 2525, 4453, 4703 transcripts were up-regulated and 2790, 2816, 3059 transcripts were down-regulated compared to their expression in the BD, WF and YF respectively, which indicated a large shift in the transcriptional program. WF and YF had the smallest number of DEGs, with 437 up-regulated and 559 down-regulated. This result was consistent with the PCA and hierarchical correlation data that the opened owers (WF and YF) showed the most similarity of transcriptomes (Fig. S1). To further investigate gene expression pro les, hierarchical clustering of all DEGs were performed as shown in Fig. 3b.
The enriched biological process GO terms of the DEGs were analyzed between different two samples. Although many of the enriched GO terms were common, a few were unique in different sets of genes. For example, the GO terms related to shoot axis formation, second shoot formation, oral organ development, fruit morphogenesis were enriched in LF vs. BD (Fig. S2a), regulation of seed dormancy process, positive regulation of mitotic cell cycle, positive regulation of seed germination, mitochondria RNA metabolic process, cellular response to cold in LF vs. WF (Fig. S2b), response to cytokine, regulation of response to red or far-red, regulation of mitotic metaphase/anaphase, photosynthetic electron transport in photosynthesis in LF vs. YF (Fig. S2c), terpenoid catabolic process, regulation of phosphate process, regulation of nucleocytoplasmic transport, regulation of G2/M transition of mitotic process in BD vs. WF (Fig. S2d), terpenoid transport, regulation of timing of meristematic phase, regulation of programmed cell death, isoprenoid transport in BD vs. YF (Fig. S2e), response to temperature stimulus, avin adenine dinucleotide binding, gibberellic acid mediated signaling pathway, RNA export from nucleus in WF vs. YF (Fig. S2f) .
The owering time pathway in L. japonica In the model plant Arabidopsis thaliana, genes involved in owering transition have been well characterized [34,35]. In this study, 70 unigenes to owering control pathway were assigned based upon sequence homology search and a detailed owering regulatory network of L. japonica was constructed (Fig. 4a). These genes were grouped into various owering control pathways, which include photoperiod pathway genes such as PHYA oral meristem identify genes LEAFY (LFY) and AP1 etc. All these were identi ed in our L. japonica RNAseq database and their expression level in different tissues was analyzed in Fig. 4b. The transcript levels of key genes of owering pathway and oral molecular networks were further analyzed in ower development by RT-qPCR (Fig. 4c). AGL24 is a dosage-dependent activator for owering. The expression of LjAGL24 mRNA was relatively high in LF, a drastic decrease in BD, and followed by gradually increase during ower development (Fig. 4b, 4c). SOC1 also functions as a owering promoter, together with AGL24, activating downstream targets [36]. LjSOC1 had the lowest expression in BD and increased 4-fold during ower development (Fig. 4b, 4c). AP1, a MADS-box ower speci c gene, is involved in owering transition and ower development. The expression of LjAP1 was only detected in the in orescence (BD, WF and YF) but not in vegetative tissues (LF) (Fig. 4b, 4c), suggesting that the function of AP1 in the speci cation of oral meristem is conserved [37]. FLC acts as a repressor of owering and is a convergence point for environmental and endogenous pathways that regulate owering time. The expression of LjFLC1 is much higher than that of LjFLC2, indicating that LjFLC1 may play the dominant role in regulating owering time. LjFLC1 exhibited the similar expression pro les of AtFLC [38], which is the highest expression level in LF and a gradual decrease during the ower maturation (Fig. 4b, 4c). Except of AGL24, SOC1, AP1 and FLC1, some other genes associated with owering time show dynamic change during ower development, such as PHYA, LFY, FT, SVP, GI, CCA1, TOC and ZTL (Fig. 4b), reinforcing the studies that these genes form a complex network to regulate the ower transition.
The oral aromas pathway in L. japonica The volatile organic components of L. japonica are often used in foods, cigarettes and cosmetics [39]. Most of the oral volatile compounds in L. japonica belonged to monoterpenes and sesquiterpenes, leading us to analyze the KEGG annotation results regarding the biosynthesis of terpenoid backbone. In this study, we identi ed 54 candidate genes involved in the MEP and MVA pathways (Fig. 5a). Most of them were highly expressed, which might indicate the terpenoid biosynthesis pathway is active and dynamic in L. japonica. DXS and AACT are the rst enzymes of the MEP and MVA pathways, respectively. The expression of some LjDXS (Trinity_DN49357_c3_g1) and LjAACT (Trinity_DN49155_c0_g3) genes were very high (FRKM value>3000) in owers, indicating their vital roles in terpene volatiles in L. japonica. Some of the most highly expressed genes centered on the MEP pathway, such as DXS, HDS, IDS, GGPPS (Fig. 5b), suggesting that this plastid pathway, might represent the major pathway for terpenoid biosynthesis in L. japonica. Within different tissues, MEP pathway genes, including DXS, HDS, IDI, GGPPS are relatively higher expressed in reproductive tissues than vegetative tissues, whereas genes involved in MVA pathway, such as HMGS, HMGR, MVK, and MVD were highest expressed in BD (Fig. 5b), at which the emission of volatile sesquiterpenoids has been shown to be high in L. japonica. To con rm the expression level of terpenoid-related genes in the MVA and MEP pathways, the expression of MVK, MVD, HDS and IDS genes at four tissues in L. japonica were tested by RT-qPCR. The results showed that the transcripts of MVK and MVD had the strongest expression in BD, whereas the transcript levels of DXS and DXR are relatively lower in leaves than in opened owers (Fig. 5c). Eleven TPS unigenes were identi ed in the current transcriptome, including unigenes for monoterpenes linalool_synthase (LIS), αterpineol synthase (TPS-Cin), geraniol synthase (GES), and sesquiterpenes (E)-nerolidol synthase (NES).
The transcriptional patterns of TPS genes were very different in L. japonica (Fig. 5b), which was reasonable to explain the divergent volatile compound emissions in various tissues.
Genes related to ower color in L. japonica Previous studies have shown that the content of chlorophyll and anthocyanin decreased dramatically during the early L. japonica ower development from the bud to the opened owers, whereas the carotenoid went up substantially [24,40]. Especially during the transformation from white owers to yellow owers, the content of carotenoid increased more than three times [40,41]. In this study, we focused on the carotenoid and anthocyanin biosynthetic pathway in L. japonica. First, we identi ed 33 genes encoding enzymes involved in carotenoid biosynthetic pathway (Fig. 6a). PSY is the rst and the rate-limiting enzyme in carotenoid skeleton production. We found that the level of LjPSY transcripts gradually increased during the developmental stage from BD to YF (Fig. 6b, 6c), which is consistent with the content of carotenoid in different developmental stages, indicating that LjPSY plays the vital role in carotenoid biosynthesis. Interestingly, the expression of some genes such as LYCE and NSY were decreased from BD to YF stage (Fig. 6b). Engineering of anthocyanin biosynthetic pathway in L. japonica is also detected. A total of 31 unigenes involved in anthocyanin biosynthetic pathway has been found to be expressed in L. japonica tissues (Fig. 6a). Some important enzymes were annotated more than 5 unigenes, such as CHS, PAL, C4H, HCT, CYP73A (Fig. 6a, 6b). The transcripts of PAL, C4H, CHS, F3'H showed the highest expression in LF, a dramatic decrease in BD and extremely low expression in the opened owers. In contrast, the transcripts of 4CL, ANS, DFR showed the highest expression in BD (Fig.  6b). To validate these data, the expression level of PAL, PSY, CHS and ANS were selected for RT-qPCR analysis. The relative expression level of these four genes were consistent with their relative abundances (Fig. 6c).

Sequencing and functional annotation
In our study, we sequenced twelve cDNA libraries from L. japonica leaves, oral buds, white ower and yellow ower with three biological replicates using the Illumina HiSeq X Ten platform, yielding 408409 transcripts and 195850 unigenes. Previous transcriptomic studies on L. japonica were mainly focused on the biosynthesis of chlorogenic acid and luteolosides, the major two bioactive metabolites [3,12,42,43], while the mechanism that control ower phenotypes of L. japonica is almost unknown. Our data focused on the expression of candidate genes involved in L. japonica oral transition and developmental processes. This result provided a genetic and molecular frame work for integrating different signaling pathways relevant to oral traits, contribution to the molecular breeding of the genus Lonciera and functional characterization of oral-related genes.
Genes related to the regulation of owering time through ve major pathways L. japonica owering is a highly coordinated, genetically programmed process by the interaction of oral transition genes, endogenous signals, and environmental factors. Previous studies have shown that the owering time was controlled by ve pathways in Arabidopsis, including photoperiod pathway, vernalization pathway, autonomous pathway, circadian clock pathway and GA pathway [15,16]. In this study, the transcripts of LjAGL24 and LjSOC1 exhibit similar pattern with a relative high expression in LF, a drastic decrease in BD, and followed by a gradually increase during ower development (Fig. 4b, 4c). In Arabidopsis, AGL24 and SOC1 functions as owering activators [36]. The reduction of AGL24 activity by RNAi results in late owering, whereas overexpression of AGL24 causes precocious owering [44]. AGL24 mRNA is present in all of Arabidopsis tissues with the strongest expression in stem, a valley in ower bud and the gradually increased expression during ower development [44]. Similarly, AtSOC1 was found to show a relatively high expression during vegetative stage, a downward trend to very low level during the transition to owering, and a gradual increase during the ower development [38]. Our results show that the expression pattern of LjAGL24 and LjSOC1 are consistent with the previous studies of AtAGL24 and AtSOC1 [38,45], indicating that AGL24 and SOC1 might play the essential and evolutionarily conserved roles in control of owering switch across the plant kingdom. AP1, a MADS-box ower speci c gene, is involved in owering transition and ower development. The expression of LjAP1 was only detected during the oral development but not in the vegetative tissues (Fig. 4b, 4c), similar to the expression of AP1 in Arabidopsis and grapevine [37,38], suggesting that the function of AP1 in the speci cation of oral meristem is evolutionarily conserved. FLC is a key regulator of owering time in Arabidopsis. The expression of AtFLC was high during the vegetative meristem development, maintained at low levels during the owering transition, and continued the downward trend to a very low level during the ower development [38]. In the present study, the expression pro les of LjFLC1 is similar to the expression pro les of AtFLC, suggesting its evolutionarily conserved roles in modulation of owering time. Genes in the circadian clock/photoperiod pathway, play critical roles by regulating the expression of CO and oral integrator FT [46]. The mRNA level of LjCO and LjFT were at a low level (FPKM < 25), since we harvest the samples in the morning and these genes were accumulated in the evening [46]. In the present study, 70 candidate genes involved in owering time pathway were identi ed (Fig. 4a) and their expression levels were analyzed (Fig. 4b), which enabled us to infer the regulatory network of owering time genes in L. japonica. More detailed analyses of additional developmental stages or tissues would be bene cial for future investigations.
Volatile terpenoid metabolism genes in L. japonica Terpenes, especially monoterpenes such as linalool, geraniol and nerol, but also some sesquiterpenes such as farnesene, nerolidol, and caryophyllene, are common constituents of oral scent in L. japonica [19]. The oral scent biosynthesis of terpenes have been studied in many plants, including Syringa oblate [47], Salvia o cinalis [48], Polianthes tuberosa [49] and Arabidopsis thaliana [50]. The previous GC/MS analysis showed that the main volatile constituents of L. japonica are farnesol (16.2%) and linalool (11.0%) for the ower fraction, hexadecanoic acid (16.0%) and linalool (8.7%) for the leaf fraction, and hexadecanoic acid (31.4%) for the stems [19]. The regulation of terpenoid biosynthesis is complex, and the pathway uxes are mostly controlled at the transcript level [51]. In this study, we identi ed 54 L. japonica genes encoding enzymes involved in terpenoid biosynthesis pathway. Our results showed that the genes involved in MVA pathway including AACT, HMGCS, MVK, and MVD had the strongest expression in oral buds, which was consistent with a strong emission of farnesol from these ower organs [19]. Moreover, our results support that farnesol is synthesized through MVA pathway as shown in Fig. 5a, and its biosynthesis and emission are closely correlated to the expression levels of those genes. Linalool was found to be an abundant volatile compound from L. japonica reproductive and vegetative tissues, accounting for 11.0% in owers and 8.7% in leaves [19]. In this study, genes involved in the MEP pathway, such as DXS, gcpE, ispS, GGPS showed very high expression levels and broad expression patterns in four tissues (Fig. 5b), which may be attributable to the relatively high linalool contents in most L. japonica tissues.
The transcription regulation of carotenoid and anthocyanin biosynthetic pathways in L. japonica During the process of L. japonica oral development, the ower color transits from the initial green to white and nally yellow in a short time frame. Metabolic analyses of pigment in L. japonica revealed that the content of chlorophyll and anthocyanin decreased drastically from bud to white owers, whereas the carotenoid pigments went up substantially from white to yellow owers [23,24], suggesting that the early ower colors were mainly caused by the decreased level of anthocyanins and chlorophylls, and late ower color primarily by the increased contents of carotenoids [24]. There is increasing evidence that carotenogenesis in plant tissues is predominantly regulated at the transcriptional level [52]. For example, it was reported that a higher level of PSY and DXS might be responsible for the marigold color development from pale-yellow to orange [53]. In Camellia nitidissima, carotenoid synthesis pathway genes PSY and crtZ exert synergistic effects in carotenoid biosynthesis during ower development [54]. Previous studies showed that the expression level of L. japonica PSY1, PDS1, PDS3 and LCYB were signi cantly higher during yellow ower stages than during early ower developmental stages [40]. In this study, 33 L. japonica unigenes were identi ed related to the carotenoid biosynthesis pathway. The mRNA level of LjPSY transcripts gradually increased during the developmental stage from oral bud to yellow owers (Fig. 6c), which is consistent with previous reports [40]. In contrast, the expression of LjLYCE and LjNSY were decreased from BD to YF stage (Fig. 6b). These results can be plausibly explained by the fact that the amount of carotenoids in the tissues is not only attributed to the synthesis of carotenoids, but also to carotenoid degradation and sink capacity [55]. Indeed, the L. japonica genes encoding carotenoid cleavage dioxygenases (LjCCD4 and LjCCD1b) exhibited stable expression during the ower development and played important roles in carotenoid degradation [40]. In L. japonica, the previous studies have identi ed 11 DEGs involved in the anthocyanin biosynthetic pathway, and all of them were most highly expressed in the early ower development [24]. In this study, we identi ed 31 unigenes encoding 10 putative enzymes involved in anthocyanin biosynthesis. The transcript level of PAL, C4H, CHS, F3'H showed the highest expression in LF, while the transcripts of 4CL, ANS, DFR showed the highest expression in BD. The expression level of most anthocyanin biosynthetic genes was dramatically decrease during the ower developmental stages (Fig. 6b), which is consistent with the sharp decline in the contents of anthocyanin from bud to the opened owers.

Conclusion
In the present study, the complete transcriptome of L. japonica was de novo-assembled and annotated, generating a total of 195850 unigenes, of which 84657 could be functionally annotated. The candidate genes involved in owering transition were identi ed and the owering regulatory network of ve pathways was constructed in L. japonica. The expression pro les of some oral integrator genes indicate that these key integrators might play the essential and evolutionarily conserved roles in control of owering switch across the plant kingdom. We also identi ed 54 L. japonica genes encoding enzymes involved in terpenoid biosynthesis pathway. Most highly expressed genes centered on the MEP pathway, suggesting that this plastid pathway might represent the major pathway for terpenoid biosynthesis in L. japonica. In addition, the key genes involved in the carotenogenesis and anthocyanin biosynthesis pathway were identi ed, which may elucidate the regulatory mechanisms of ower color transitions in L. japonica. These results provide insights into the molecular mechanisms underlying the owering time, oral scent and ower color, and will be useful for molecular breeding of L. japonica. Availability of data and materials All the data generated or analysed during this study are included in this published article and its supplementary information les. The datasets used and/or analysed during the current study are available from the corresponding author (s.y0001@163.com) on reasonable request.

Consent for publication
Not applicable.
Competing interests