Transcriptomic and metabolomic analyses provide insights into the biosynthesis of chlorogenic acids in Lonicera macranthoides Hand.-Mazz

Lonicera macranthoides Hand.-Mazz (L. macranthoides) is a medicinal herb that is widely distributed in South China. The developmental stage and corolla dehiscence of the flower are the important factors affecting the quality of medicinal ingredients. However, neither the regulatory mechanism controlling chlorogenic acids biosynthesis in L. macranthoides nor the molecular basis of effect of corolla dehiscence on the quality of medicinal materials is fully understood. In this study, metabolomics and transcriptomics were used to analyze the metabolic and transcriptional differences of two different cultivars closed bud type (Bt), and flowering type (Ft), as well as the effect of jasmonic acid methyl ester (MeJA) on chlorogenic acids (CGAs) biosynthesis. In total, large number of differentially expressed genes (DEGs) and differentially accumulated metabolites (DAMs) were filtered among three lines of samples. Gene metabolite correlation analyses revealed a ‘core set’ of 30 genes and 54 genes that were strongly correlated with CGAs biosynthesis and regulating the flowering, respectively. Quantitative real-time polymerase chain reaction results proved the alterations in the expression levels of genes encoding the pathways involved in CGAs biosynthesis. The ion abundances of CGAs were most significantly increased, while some of the CGAs derived and Caffeoyl-CoA-derived substances showed the most largely reduced abundances in the closed bud type (Bt) compared to the flowering type (Ft). MeJA may leads to the activation of downstream genes in CGAs biosynthesis pathway. Overall, there were significant differences in the transcriptional and metabolic levels of CGAs biosynthesis pathway in flower buds of different flowering cultivars. The redirection of metabolic flux may contribute to increased accumulation of CGAs. However, whether MeJA and flowering have direct effects on the accumulation of CGAs needs further studied. These researches effectively expanded the functional genomic library and provide new insights into CGAs biosynthesis in L. macranthoides.

Introduction of quinone chalcones in sunflower [30]. However, whether MeJA has any effect on the synthesis of CGAs in plants has not been reported.
L. macranthoides is widely distributed in South China and is highly adaptable to growth, but evolution has led to significant changes in yield and quality. For a long time, the short bud period and inability to pick buds in time were problems in the production and application of L. macranthoides. In this study, flowering type (Ft), closed bud type (Bt) and MeJAtreated closed bud type (Bt1000) (Fig 1) were used as experimental materials to conduct transcriptomic and metabolome comparative analysis of flowers to identify candidate genes involved in CGAs biosynthesis. Because of MeJA can effectively regulate corolla opening, we used MeJA-treated closed bud type (Bt1000) to test whether flowering was related to the biosynthesis of CGAs. This is the new attempt to investigate the transcriptional and metabolite differences between the two different flowering cultivars of L. macranthoides, as well as the effect of flowering on chlorogenic acids biosynthesis, and the results will enhance our understanding of the metabolic mechanisms and molecular basis of CGAs biosynthesis in L. macranthoides.

Plant materials
Ft' is a flowering type of L. macranthoides that is a local conventional cultivar in Chongqing, China. 'Bt' is a closed bud type of L. macranthoides that was bred from 'Ft' and widely cultivated in Chongqing, China. The corolla of Bt does not crack until it falls, and the Ft's bud period was shorter, followed by corolla cracking. 'Ft' and 'Bt' were planted in parallel and were all 8 years old. The experiment was conducted in the southern hilly region of Banan District, Chongqing Municipality, China (N29.34, E106.92). On sunny days, we collected yellow white buds (Ft, preflowering period) from Ft. In Bt, the yellow white buds were collected first, and then the remaining yellow white buds in the same branch were sprayed with MeJA (1000 ppm) by miniature sprayer. After 24 hours, the opened flowers (MeJA-treated closed bud type, Bt1000) were collected and used as samples. Three biological replicates were performed in each case. After washing, the samples were immediately frozen in liquid nitrogen and stored at -80˚C until further use in metabolite, RNA sequencing (RNA-Seq), and qPCR analyses.

UHPLC and ESI-Q TRAP-MS/MS conditions
LIT and triple quadrupole (QQQ) scans were acquired on a triple quadrupole-linear ion trap mass spectrometer (Q TRAP), API 4500 Q TRAP LC/MS/MS System, equipped with an ESI Turbo Ion-Spray interface, operating in positive ion mode and controlled by Analyst 1.6.3 software (AB Sciex). The ESI source operation parameters were as follows: ion source, turbo spray; source temperature 550˚C; ion spray voltage (IS) 5500 V; ion source gas I (GSI), gas II (GSII), and curtain gas (CUR) were set at 55, 60, and 25.0 psi, respectively; and the collision gas (CAD) was high. Instrument tuning and mass calibration were performed with 10 and 100 μmol/L polypropylene glycol solutions in QQQ and LIT modes, respectively. QQQ scans were acquired as MRM experiments with collision gas (nitrogen) set to 5 psi. DP and CE for individual MRM transitions were performed with further DP and CE optimization. A specific set of MRM transitions was monitored for each period according to the metabolites eluted within this period.

Metabolite identification and quantification
Metabolite data analysis was conducted with Analyst 1.6.1 software (AB SCIEX, Ontario, Canada). The supervised multivariate method, partial least squares-discriminant analysis (PLS-DA), was used to maximize the metabolome differences between a pair of samples. The relative importance of each metabolite to the PLS-DA model was checked using the parameter called variable importance in projection (VIP). Metabolites with VIP � 1 and fold change �2 or fold change �0.5 were considered differential metabolites for group discrimination [31].

RNA sequencing and annotation
Total RNA isolation and purification and cDNA library construction and sequencing were performed as previously described [32,33]. Total RNA was prepared by using the TRIzol reagent (Invitrogen) following the manufacturer's protocol. Total RNA purity and concentration were determined using the NanoPhotometer 1 spectrophotometer (IMPLEN, CA, USA) and the Qubit 1 RNA Assay Kit in Qubit 1 2.0 Fluorometer (Life Technologies, CA, USA). mRNA was isolated from total RNA using magnetic beads with oligo (dT) primers; cDNA was synthesized using a cDNA synthesis kit (TaKaRa, Dalian, China) and linking the sequencing adapter to both ends. The library preparations were sequenced on an Illumina HiSeq 4000 platform. Transcriptome assembly was accomplished based on the left.fq and right.fq using Trinity [34,35] with min_kmer_cov set to 2 by default and all other parameters set default. Gene function was annotated based on the following databases: Nr (NCBI nonredundant protein sequences); Nt (NCBI nonredundant nucleotide sequences); Pfam (Protein family); KOG/ COG (Clusters of Orthologous Groups of proteins); Swiss-Prot (A manually annotated and reviewed protein sequence database); KEGG (Kyoto Encyclopedia of Genes and Genomes Ortholog database) [36]; and GO (Gene Ontology) [37]. The fragments per kilobase million (FPKM) were calculated to quantify gene expression [38], and the DESeq2 R package [39] was adopted to determine the differentially expressed genes (DEGs). The differentially expressed genes were identified with the following parameters: corrected P-value of 0.05 and absolute fold change �2. GO and KEGG enrichment analyses of differentially expressed genes were further implemented employing the clusterProfiler R package [40].

Quantitative Real-Time PCR (qRT-PCR) analysis
For the expression analysis of the DEGs, total RNA was extracted with TRIzol reagent (Invitrogen, CA, USA) from 100 mg of all the samples. First-strand cDNA was synthesized from 2 μg of total RNA using M-MLV Reverse Transcriptase (Promega, WI, USA) according to the manufacturer's instructions. The reactions (20 μL) were terminated after 40 cycles using the CFX96TM Real-Time PCR Detection System (Bio-Rad, CA, USA), and the 18S reference gene was used as an internal control. All of the assays were repeated at least three times. The relevant PCR primer sequences are shown in S1 Table and were designed based on the CDS and ESTScan sequences of L. macranthoides by Primer Express Software (Applied Biosystems, CA, USA).

Correlation analysis between metabolome and transcriptome data
Pearson's correlation coefficients were calculated between the metabolome and transcriptome data. The coefficients were calculated from log2 (fold change) of each metabolite and log2 (fold change) of each transcript with the Excel program. Correlations with a coefficient of R2 > 0.8 were selected.

Sequence accession numbers
The whole set of transcript data can be found in the National Center for Biotechnology Information (NCBI) SRA database (SRR12727863-SRR12727871).

Plant secondary metabolomics analysis
In general, the best harvesting period for L. macranthoides was the yellow white bud stage. In this study, the flowering type (Ft), closed bud type (Bt) and MeJA-treated closed bud type (Bt1000) were used to monitor changes in the metabolite abundance associated with the biosynthesis of active substances in L. macranthoides (Fig 2). We described the metabolome of the three lines of samples using the widely targeted metabolomics approach. In total, 341 metabolites were detected and quantified from L. macranthoides, including 64 phenolic acids, 138 flavonoids, 28 lignans and coumarins, 2 tannins, 26 alkaloids, 23 terpenoids, 35 organic acids, and 25 other metabolites (Fig 2a). In the clustering heat map, the accumulation of metabolites displayed clear phenotypic variation in terms of the pattern of metabolite abundance in different types, and the Ft was clearly distinct from the Bt (Fig 2b). In this study, the concentration data of metabolites were used for the hierarchical heat map clustering analysis of samples. We observed that all the biological replicates were grouped together (top side of the figure), indicating a high reliability of the metabolome data (Fig 2b). In particular, we observed clear separations among Ft, Bt and Bt1000. The above indicates that there are obvious differences in the metabolic characteristics of them.

Comparison of metabolites produced by three lines of samples (Bt, Ft, Bt1000)of L. macranthoides
To explore the changes in secondary metabolites in Bt, Ft, Bt1000, principal component analysis (PCA) and orthogonal partial least-squares discrimination analysis (OPLSDA) were applied to the metabolites, and the three lines of samples showed significant separations (Figs 3 and 4). From the results, they were clearly distinguished in the PC1 dimension of the PCA score graph (37.74% variation) and were further differentiated at the PC2 dimension, indicating that the

PLOS ONE
Analyses provide insights into the biosynthesis of chlorogenic acids in Lonicera macranthoides Hand.-Mazz secondary metabolites varied significantly with different types and treatments. The OPLS-DA model compared the entire metabolite content of the core samples in pairs to evaluate the differences between Bt1000 and Bt (R 2 Y = 1, Q 2 Y = 0.986; Fig 4a) and between Bt and Ft (R 2 Y = 1, Q 2 Y = 0.998; Fig 4b). High predictability (Q2) of the OPLSDA models was stable and reliable and could be used to further screen for differential metabolites. The differentially accumulated metabolites (DAMs) between pairs of samples (Bt vs Bt1000 and Bt vs Ft) were determined based on the variable importance in projection (VIP) � 1 and fold change � 2 or fold change �0.5 [41]. There were 75 significantly different metabolites between Bt and Bt1000 (Bt1000 had 38 downregulated and 37 upregulated metabolites) and 103 between Bt and Ft (Ft had 43 downregulated and 60 upregulated metabolites) (Fig 4c and 4d). Furthermore, 13 common differential metabolites were detected in all three comparison groups (Fig  3b). For Ft vs Bt, 103 DAMs were detected and quantified, including 12 phenolic acids, 62 flavonoids, 12 lignans and coumarins, 1 tannins, 2 alkaloids, 2 terpenoids, 7 organic acids, and 5 others. The metabolite with the largest fold change was isorhamnetin 3-O-neohesperidoside and the top 20 fold change metabolites were mainly Flavonoids (S9 Table). For Bt vs Bt1000, 75 DAMs were detected and quantified, including 11 phenolic acids, 26 flavonoids, 5 lignans and coumarins, 2 tannins, 9 alkaloids, 2 terpenoids, 10 organic acids, and 10 others. The metabolite with the largest fold change was Nomilin and the top 20 fold change metabolites mainly fell in other categories, mainly Flavonoids and Organic acids (S10 Table).
KEGG is the main public database of metabolic pathways and is used for the research of genes, expression information, and metabolite content in a general network. The top enriched KEGG terms among the DAMs detected for Bt vs Ft were biosynthesis of secondary metabolites, phenylpropanoid biosynthesis and tryptophan metabolism (Fig 5a). Different from Bt vs Ft, the top enriched KEGG terms among the DAMs detected for Bt vs Bt1000 were biosynthesis of secondary metabolites, tryptophan metabolism and neuroactive ligand-receptor interaction (Fig 5b).

Transcriptome profiles and major transcriptomic differences
To investigate the transcriptomic changes among Bt, Ft, Bt1000, nine cDNA libraries were generated and sequenced using an Illumina deep-sequencing HiSeq™ 4000. Among all the raw reads, 97% had Phred-like quality scores at the Q20 level (an error probability of 1%). After removing low-quality sequences, 424, 085, 868 clean reads and 63.63 GB clean bases were used to assemble the transcriptome data using the Trinity method (Table 1) Table 2).
The FPKM method [35] was used to calculate the expression of all unigenes to remove the effects of length differences and sequencing depth. DEGs are defined as genes that are significantly enriched or depleted in one sample relative to another by DESeq [38] (qvalue < 0.005 and log2 (fold change) > 1). The results of transcriptome analysis showed that the gene expression levels were different in three lines of samples of L. macranthoides. The clustering heat map shows that there were obvious differences in gene expression in them. There was overlap between Bt1000 and Bt with lower expression levels, while there was almost no overlap between Bt and Ft (Fig 6a). The Venn diagram of the differential genes shows the number of common and specific genes of three lines of samples compared in pairs. From Fig 6b, the expression levels of 1, 677 genes were different in 3 different types of samples, with 13, 408 and 1, 976 having specific differences, respectively. The volcano map (Fig 6c) revealed that the expression levels of 10, 836 and 4, 249 DEGs were higher or lower in Ft than in Bt, while the expression levels of 2, 954 and 699 DEGs were higher or lower in Bt1000 than in Bt. The KEGG enrichment analysis of the DAMs indicated that most of the identified DEGs act on multiple metabolic processes related to terpenoid backbone biosynthesis, phenylpropanoid biosynthesis, sesquiterpenoid and triterpenoid biosynthesis and stilbenoid, diarylheptanoid and gingerol biosynthesis in Bt1000 vs Bt (Fig 6e). The DEGs of Ft vs Bt were enriched on terpenoid backbone biosynthesis, pentose and glucuronate interconversions and starch and sucrose metabolism (Fig 6f). The phenylpropanoid biosynthesis which we focused on enriched 73 DEGs. In total, 6 DEGs were selected for qRT-PCR quantification experiments to verify the

PLOS ONE
Analyses provide insights into the biosynthesis of chlorogenic acids in Lonicera macranthoides Hand.-Mazz accuracy of transcriptome sequencing. The results indicated that the expression patterns from qRT-PCR testing were well correlated with sequencing results (S1 Fig), indicating the reliability of the transcriptomic data.

Differentially expressed genes related to flowering in three lines of samples (Bt, Ft, Bt1000) of L. macranthoides
Ft is prone to corolla dehiscence in the florescence period, and this period is relatively short, generally 3-7 days. If the florescence period is overcast and rainy, it is easy to fall flowers,  Table), indicating that these FRGs mainly respond to MeJA treatment through positive feedback regulation. In our experiment, the transcriptional abundance of flowering-related genes was significantly altered under MeJA treatment. Some genes in the CDF family, especially i1_LQ_LMFt_c43539/f1p0/1681 and i1_LQ_LMFt_c60085/f1p0/1664, which were upregulated by 125-fold and 106-fold, respectively, were heavily programmed and highly induced under MeJA treatment (S2 Table) Table), respectively. The above results indicate that MADS-Box and other genes play an important role in regulating the flowering of L. macranthoides, and further research will be carried out.

DEG and DAM in the chlorogenic acids biosynthesis pathway between Bt and Ft
It has been reported that CGAs is the most interesting of the biologically active ingredients of L. macranthoides [1]. As the content of CGAs in Ft and Bt was significantly different over the course of this analysis, KEGG pathways involved in CGAs biosynthesis were further analyzed in this study. CGAs biosynthesis mainly involves the phenylalanine biosynthesis pathway, which is very complex and involves the regulation of various enzymes, such as PAL, 4CL, C4H, HCT and C3H (Fig 7a). According to KEGG enrichment results and FPKM values, 30 key enzyme genes differentially expressed in CGAs biosynthesis were successfully screened from DAMs, including PAL (10), 4CL (8), C4H (5), HCT (4), and C3H (3) (S4 Table). The data indicated that all CGAs biosynthesis pathway genes showed significantly higher levels of gene transcripts in flowering type (Ft). Six of the ten PAL genes involved in the synthesis of cinnamic acid were not expressed in closed bud type (Bt), two of which exhibited a higher level of expression in Ft. In 4CL genes, i1_LQ_LMFt_c23140/f1p7/1891 was not expressed in Bt. The expression level of most of the 4CL genes in Ft was upregulated, with the highest difference being a 24-fold increase. The C4H-2 was not expressed in Bt, while the other C4H were upregulated in Ft. Especially in HCT genes, there were two genes with high differential expression multiples, up to 288 and 775 times (Fig 7b). However, interestingly, although these genes were upregulated in Ft, the ion abundance of most CGAs in Ft were significantly lower than that of Bt (Table 3). We know that in the phenylalanine biosynthesis pathway, the synthesis reaction does not stay on CGAs, and the next step product is synthesized under the catalysis of other enzymes. In our experiment, a total of 8 DAMs were confirmed based on UPLC-MS, as shown in Fig 7c. Apart from cinnamic acid, the 7 DAMs are downstream of the phenylalanine biosynthesis pathway. Upstream of CGA synthesis, the ion abundance of eriodictyol and naringenin in Ft was significantly higher than that in Bt. Compared with Bt, the ion abundance of Rutin, Coniferin and Scopolin in Ft was up-regulation (Fig 7c). To some extent, although the upstream genes of CGAs biosynthesis are upregulated in Ft, due to the generation of caffeoyl-COA and CGAs are catalyzed by corresponding enzymes to synthesize downstream substances, resulting in the ion abundance of CGAs in Ft being significantly lower than that of Bt.

DEG and DAM in the chlorogenic acids biosynthesis pathway between Bt and Bt1000
In the preliminary study of our group, it was found that MeJA treatment at a certain concentration could effectively promote the corolla opening of Bt. To determine whether MeJA affects the synthesis of CGAs and the relationship between flowering and synthesis of CGAs, we treated the buds of Bt with MeJA and carried out studies by means of metabolomics and transcriptomics analysis methods. In our study, the genes related to CGAs biosynthesis in Bt showed an upregulation under the treatment with MeJA. Among these DEGs, i2_LQ_LMFt_c5476/f1p0/2700 (annotated as PAL) and i1_LQ_LMFt_c48921/f1p17/1795 (annotated as C3H) were upregulated by 17-fold and 37-fold (S5 Table), respectively. We detected that the upstream substance in Bt1000 was higher than that in Bt. However, the ion abundance of CGAs was slightly lower than that of Bt. We found that the alterations of Bt1000 compared to Bt were increased in cinnamic acid and coniferin, and decreased in caffeic acid and coniferyl alcohol (S6 Table). From the above results, MeJA leads to the activation of downstream genes in the phenylalanine biosynthesis pathway, and more substances flow downstream. However, whether MeJA and tube cracking have direct effects on the biosynthesis of CGAs needs to be further studied.

Correlation analysis between DEGs and DAMs in other metabolic pathways
According to the results of differential metabolite analysis in this study and combined with the conjoint analysis of DEGs and DAMs results, the genes and metabolites in the same group were simultaneously mapped to the KEGG pathway map to better understand the relationship between genes and metabolites. We found that the alterations of Ft compared to Bt were also associated with 12 metabolic pathways, including pentose and glucuronate interconversions, ascorbate and aldarate metabolism, ubiquinone and other terpenoid-quinone biosynthesis, phenylpropanoid biosynthesis and so on (S7 Table). The alterations of Bt compared to Bt1000 were also associated with 20 metabolic pathways (S8 Table). In the correlation analysis, Dgalacturonic acid, cinnamic acid and succinic acid were highly correlated with the corresponding metabolic pathways, which also provides a new perspective for our subsequent research.

Discussion
In general, L. macranthoides's buds or open flowers are used as medicine [2]. The flowering type (Ft)'s bud period was shorter, and after corolla dehiscence, the flowers were easy to drop in cloudy and rainy days which brings considerable trouble to their harvest [42]. The closed bud type (Bt) is a new cultivar with non-cracking corolla, and the CGAs abundance is significantly higher than that in Ft; therefore, Bt is recommended in production. With the flowering process, CGAs content decreases significantly. The abundance of active components and volatile compounds are closely correlated with floral developmental stages [43]. Moreover, MeJA can effectively regulate corolla opening. The quality of herbal medicine has been very difficult to control and evaluate primarily because of the complexity and incomplete knowledge of the active medicinal compounds [42]. To explore the synthesis mechanism of CGAs in different flowering types and MeJA-treated Bt, a comprehensive analysis of the transcriptome and metabolome was performed. CGAs is the most interesting of the biologically active ingredients of L. macranthoides [10]. CGAs are the primary phenylpropanoid generated from the shikimic acid pathways with high antioxidant activities. In our results, 341 metabolites were detected and quantified by UHPLC-ESI-MS/MS with a self-built database. We profiled the chemical compositions of L. macranthoides in a systematic and comprehensive way, which provided a reference for their utilization in the future. Studies on L. japonica, coffee, artichoke and other plants have revealed that PAL, C4H and 4CL are the common enzymes in the upstream CGAs metabolic pathway [2,41,44,45]. DEGs related to CGAs biosynthesis were identified between Bt and Ft, including ten PAL, eight 4CL, five C4H, four C3H and four HCT genes. Among them, we revealed that most of the enzyme genes in the CGAs biosynthetic pathway were highly active in Ft. Similarly, in terms of metabolites, the upstream substances of CGAs, coumaric acid and cinnamic acid, were upregulated in Ft. In contrast, the content of CGAs in the phenylalanine synthesis pathway was significantly different. The CGAs in Bt was significantly higher than that in Ft. On the basis of previous research, we drew a putative road map of the phenylalanine biosynthesis pathway. Transcriptomics integrated with metabolomics revealed that most components of downstream material showed drastically higher abundances in Ft than in Bt, accompanied by significant reductions in the CGAs. The above results show that in Ft, the particle flow does not stay at the CGAs stage but is synthesized downstream at a high rate. Similar results have also been reported in ginkgo and peanut species [46,47]. An organism is complex. A metabolite may be regulated by multiple genes, and a gene may regulate multiple metabolites, so there may not be one-to-one relationship between metabolites and genes.
The correlation between flowering and chlorogenic acids biosynthesis is also discussed. MeJA is an endogenous growth regulator in plants and has a wide range of physiological functions. From previous research results, MeJA can regulate the synthesis of flavonoids in plants [44,47]. However, the effect of MeJA on the biosynthesis of phenolic acids is still unclear. At the same time, whether there is any correlation between flowering and chlorogenic acids biosynthesis has not been reported. Most of the 16 FRGs were upregulated under MeJA treatment, indicating that these FRGs mainly respond to MeJA treatment through positive feedback regulation. A total of 38 DEGs, including the MADS-Box, FRI, FCA, ELF and SVP families, were found in the typical flower type (Ft) and the close bud type (Bt) of L. macranthoides. Among these DEGs, 28 DEGs were upregulated in Ft, and 13 DEGs were not expressed in Bt. MADSbox family genes in plants are important transcriptional regulators and are involved in the regulation of flower organ development and flowering time [48][49][50]. Different MADS-Box family genes showed different expression trends in Ft and Bt, indicating that this gene family had both positive and negative regulation on the flowering of L. macranthoides. From the results of this study, we screened 14 genes related to 'early flowering', one of which was mostly upregulated in Ft. In particular, the expression of this kind of gene also showed an upregulation in MeJA-treated samples, indicating that the study of this kind of gene will provide us with a new idea for the analysis of L. macranthoides flowering time. The above results indicated that the expression of genes regulating flowering was significantly different in different samples. However, by analyzing the metabolites of the three lines of samples, no direct effect of flowering genes on chlorogenic acids biosynthesis was found. We usually think that flowering affects the content of CGAs in L. macranthoides, but no direct evidence was found to prove this. Studies have shown that the endogenous hormone ethylene plays a key role in the flowering and aging process of plants [51]. In subsequent studies, the correlation between hormone metabolism, flowering and biosynthesis of active substances will be studied by combining transcriptomics and metabolites.
Currently, studies have shown that genes mentioned in this article are affected by MeJAtreated, and the abundance of metabolites also varies, but the mechanism of action is unclear. Therefore, we plan to focus on the molecular mechanism of MeJA regulating flowering and CGAs biosynthesis. At the same time, the molecular regulation mechanism of CGAs biosynthesis will be further explored to provide a better theoretical basis for the breeding of L. macranthoides varieties and the preservation of excellent genes.

Conclusion
Through conjoint metablomic and transcriptomic analysis of Ft, Bt and Bt1000, we concluded that: (1) 17, 061 DEGs and 341 DAMs among Bt, Ft and Bt1000 were identified. (2) Based on correlation analyses between DEGs and DAMs, we found 30 candidate genes for the accumulation of CGAs. We further identified approximately 54 candidate genes for regulating the flowering. (3) The redirection of metabolic flux may contribute to increased accumulation of CGAs in Bt compared to Ft. MeJA treatment may leads to the activation of downstream genes in the phenylalanine biosynthesis pathway. However, whether the MeJA and flower tube cracking have direct effects on the accumulation of CGAs needs to further studied. These results provide novel insights into CGAs biosynthesis. The candidate genes for CGAs biosynthesis and flowering presented here represent a valuable data set for future functional studies.
Supporting information S1