Enhanced Biosynthesis of Chlorogenic Acid and Its Derivatives in Methyl-Jasmonate-Treated Gardenia jasminoides Cells: A Study on Metabolic and Transcriptional Responses of Cells

Chlorogenic acid and its derivatives (CQAs) are considered as important bioactive secondary metabolites in Gardenia jasminoides Ellis (G. jasminoides). However, few studies have investigated the biosynthesis and regulation of CQAs in G. jasminoides. In this study, methyl jasmonate (MeJA) was used to enhance CQAs accumulation in cultured G. jasminoides cells. Moreover, the possible molecular mechanism of MeJA-mediated accumulation of CQAs is also explored. To this end, time-course transcriptional profiles of G. jasminoides cells responding to MeJA were used to investigate the mechanism from different aspects, including jasmonate (JAs) biosynthesis, signal transduction, biosynthesis of precursor, CQAs biosynthesis, transporters, and transcription factors (TFs). A total of 57,069 unigenes were assembled from the clean reads, in which 80.7% unigenes were successfully annotated. Furthermore, comparative transcriptomic results indicated that differentially expressed genes (DEGs) were mainly involved in JAs biosynthesis and signal transduction (25 DEGs), biosynthesis of precursor for CQAs (18 DEGs), CQAs biosynthesis (19 DEGs), and transporters (29 DEGs). Most of these DEGs showed continuously upregulated expressions over time, which might activate the jasmonic acid (JA) signal transduction network, boost precursor supply, and ultimately stimulate CQAs biosynthesis. Additionally, various TFs from different TF families also responded to MeJA elicitation. Interestingly, 38 DEGs from different subgroups of the MYB family might display positive or negative regulations on phenylpropanoids, especially on CQAs biosynthesis. Conclusively, our results provide insight into the possible molecular mechanism of regulation on CQAs biosynthesis, which led to a high CQAs yield in the G. jasminoides cells under MeJA treatment.


INTRODUCTION
Gardenia jasminoides Ellis (G. jasminoides) is an evergreen shrub of family Rubiaceae, which is widely distributed in Southern China. Its fruits are commonly used in foods and herbal medicines in China (Debnath et al., 2011). Up to now, various bioactive components have been isolated and identified from the fruits of G. jasminoides, including polyphenols, iridoid glycosides, triterpenes, flavonoids, and essential oil (Zhou et al., 2010;Han et al., 2015). Among them, chlorogenic acid and its derivatives (CQAs), belonging to depsides of certain transcinnamic acids and quinic acid, are considered as important bioactive components in G. jasminoides. In addition, CQAs exhibit various biological activities, such as anti-inflammatory, hypolipidemic, antibiotic, and antioxidant properties, which mark them as being medically intrinsic as well (Nishi et al., 2013;Ali et al., 2017;Huang et al., 2017).
Due to the lack of reference genome information, it is time-consuming and laborious to investigate gene function by the traditional method in non-model organisms (Dai et al., 2011). However, RNA sequencing (RNA-seq), as a technology to obtain almost all transcriptional information of samples under certain physiological conditions, can prospectively clarify the gene function and metabolic regulation mechanism in nonmodel organisms due to its characteristics of "high throughput, low cost, covering a multitude of low abundance gene sequencing depth, and high sensitivity" (Todd et al., 2016). Rai et al. (2017) conducted RNA-seq on nine tissues of L. japonica and identified potential candidate genes that may participate in the biosynthesis of CQAs. Comparative transcriptome analysis was performed between Eucommia ulmoides with high and low CQAs contents, which looked into potential structural genes and transcription factors (TFs) involved in CQAs biosynthesis (Ye et al., 2019).
Previously, a few studies on RNA-seq of G. jasminoides have been reported. Some of them focused on the physiological changes in petal senescence or melatonin treatment (Tsanakas et al., 2014;Zhao et al., 2017), and others focused on the biosynthetic pathway of carotenoids or crocins Ji et al., 2017). Up to now, there have been few studies on CQAs biosynthesis in G. jasminoides, especially on the mechanism of MeJA-mediated CQAs biosynthesis.
In the current study, MeJA (200 µM) was utilized to increase the yield of CQAs in the cultured G. jasminoides cells. Additionally, the cells treated for 0 h (G0h), 8 h (G8h), 20 h (G20h), and 40 h (G40h) were selected for RNA-seq to investigate the time-course transcriptional profiles in response to MeJA elicitation. Moreover, the abundant transcriptional information was analyzed in detail to gain insight into the transcriptional changes of cells treated by MeJA. In short, this work would provide valuable information to understand the molecular mechanisms of MeJA elicitation for high CQAs production in the cultured G. jasminoides cells.

Suspension Culture of G. jasminoides Cells
The embryogenic G. jasminoides callus used in this study was obtained from Jiangxi Key Laboratory of Natural Products and Functional Food (Jiangxi Agricultural University, Nanchang, China). G. jasminoides cell suspension culture was established according to the method we reported previously (Liu et al., 2018). The detailed procedures were as follows: Murashige and Skoog (MS) medium was supplemented with 0.3 mg L −1 kinetin, 0.5 mg L −1 α-naphthylacetic acid, and 30 g L −1 sucrose (pH = 5.8) for the suspension cell culture. Six grams (fresh weight) of callus was inoculated into a 100-ml culture flask containing 40 ml MS medium, subsequently cultivated in a rotating shaker with 115 r min −1 at 28 ± 1 • C under continuous light (illuminance = 1,000 l×), and routinely subcultured every 10 days. After 10 times of subculture, a homogeneous and rapidly growing suspension culture was obtained and used for further experiments.

Cell Elicitation Method for High CQAs Yield
Methyl jasmonate was dissolved in 75% ethanol and sterilized by filtration prior to elicitation. G. jasminoides cells were cultivated in the same medium and conditions as mentioned above. On the fifth day after inoculation, the aforesaid MeJA solution was added into the culture medium at a final concentration of 200 µM, and then the cells were cultured until the scheduled harvest time (0-80 h). The harvested cells were frozen immediately in liquid nitrogen and stored at −80 • C prior to use. Each treatment was carried out in triplicate.

RNA Extraction, Library Construction, and Illumina Sequencing (RNA-seq)
Four samples, in which cells were treated by MeJA for 0, 8, 20, and 40 h, were prepared by mixing three replicate samples of the harvested cells before RNA extraction. Total RNA was extracted using a Total RNA Extractor (SanGon, Shanghai, China) according to the manufacturer's instructions. mRNA was enriched by Oligo(dT) and used to construct a nonstranded RNA-seq transcriptome library. After fragmentation, reverse transcription, 3 add A, adaptor ligation, and PCR amplification, the cDNA library was obtained and then subjected to the Illumina HiSeq X Ten (PE150) platform with an average insert length of 350 bp for paired-end sequencing with one technical replicate. In addition, total cDNA libraries of cell samples exposed to MeJA for 0, 8, 20, and 40 h were designated as "G0h, " "G8h, " "G20h, " and "G40h", respectively, in the current study.

De novo Assembly and Gene Functional Annotation
To obtain the clean reads, raw reads were filtered by removing adapter sequences, low-quality reads (quality value < 20), and the reads with a length less than 35 bp. Prior to assembly, each library mentioned in section "RNA Extraction, Library Construction, and Illumina Sequencing (RNA-seq)" was pooled. Then, de novo assembly was conducted by software Trinity (version number: 2.4.0) with a min_kmer_cov parameter of 2. High-quality reads were subjected to module Chrysalis to obtain unique contigs, which were clustered to generate Bruijn. Module Butterfly was used to process Bruijn to obtain fulllength transcripts with alternative splicing and to distinguish paralogs. After transcript de-redundancy, the longest ones of each transcript cluster were regarded as isoforms representing unigenes (Grabherr et al., 2011).
All assembled unigenes were used as query sequences to search against publicly available protein databases, including National Center for Biotechnology Information (NCBI) non-redundant protein sequences (NR) and Swiss-Prot, for homology comparison using NCBI Blast+ program (version number: 2.28) with identity set at >30% and a cutoff E-value of 10 −5 (Conesa et al., 2005). By using WEGO software, the assembled unigenes were annotated on Gene Ontology (GO) and Clusters of Orthologous Groups of proteins (COG/KOG) for gene function prediction and classification (Ye et al., 2006). The KEGG Automatic Annotation Server (KAAS, version number: 2.1) was used to assign unigenes based on the Kyoto Encyclopedia of Genes and Genomes (KEGG) database for the secondary metabolic pathway annotation (Kanehisa et al., 2007). Annotation of TFs and transcriptional regulators was conducted by blasting the assembling unigenes against PlnTFDB V3.0 1 with E-values < 1e−5. When a conflict occurred among the database results, the priority order of alignments was Nr, Swiss-Port, KEGG, GO, and then COG.

Analysis of Differentially Expressed Genes
To normalize the gene expression levels, the fragment per kilobase of exon model per million mapped reads (FPKM) method was adopted on the RSEM software (version number: 1.0) (Mortazavi et al., 2008). The differential expression analysis of any two groups was analyzed using the previously reported method (Audic and Claverie, 1997) with Benjamini and Hochberg false discovery rate (FDR). Log 2 [fold change (FC)] > 1 (FDR < 0.001) was set as the threshold to identify upregulated differentially expressed genes (DEGs), and Log 2 (FC) < −1 (FDR < 0.001) was set as the threshold to identify downregulated DEGs (Reiner et al., 2003). To further analyze DEGs involved in secondary metabolism, KEGG significant enrichment was performed by comparing DEGs to the entire background of all KEGG-associated unigenes with a hypergeometric test (p < 0.05, the p-value of the hypergeometric test represents the pathway significance; in other words, the lower the p-value, the more significant the pathway changes), and then significantly overrepresented metabolic and signal transduction pathways were identified and ranked . Comparative analyses between G0h and G8h, G20h, or G40h were designated as "G0h vs G8h, " "G0h vs G20h, " and "G0h vs G40h, " respectively, in the current study.

Validation of DEG Expression Level by qRT-PCR
Total RNA was isolated from G. jasminoides cells elicited by MeJA for different periods (0,8,20, and 40 h) using TRIzol R reagent (Invitrogen) following the manufacturer's recommendations. After treatment with DNase I (SanGon, China), 1 µg of RNA was used in reverse transcription with Revert Aid Premium Reverse Transcriptase kits (Thermo Scientific TM , China). qRT-PCR was performed using SG Fast qPCR Master Mix (ABI, United States) and ran on Step One Plus Real-Time PCR Detection System (ABI, United States). The cycling conditions were initial denaturation at 95 • C for 3 min, followed by 45 cycles of 95 • C for 3 s, 60 • C for 30 s, and 72 • C for 15 s. Fifteen selected unigenes were validated by qRT-PCR experiments, and qRT-PCR reactions were carried out in triplicate. Gene-specific primer pairs (Supplementary Table 1) were designed by Primer 6.0 software, and a polyubiquitin (UBQ5) gene (c15827_g1) served as the reference gene because it was relatively stable across all time points. Gene relative expression was calculated by the 2 − CT method (Livak and Schmittgen, 2001), and all data were expressed as means ± SD after normalization. To evaluate the quality of RNA-seq data, correlation analyses of gene relative expression between qRT-PCR and RNA-seq were performed with Pearson correlation coefficient (two-tailed test, confidence interval = 95%).

CQAs Accumulation in MeJA-Treated G. jasminoides Cells
The effect of MeJA on intracellular accumulation of CQAs was examined at 80 h after elicitation ( Figure 1A). After 8 h, MeJA began to stimulate a rapid increase in CQAs. These CQAs remained in relatively high concentrations until the 72nd hour of elicitation, the maximum of which was 21.7 times (26.26 mg g −1 ) more than that of control, whereas total CQAs content varied in a little range (0.71-1.21 mg g −1 ) without MeJA elicitation (CK). The selected cell samples were exposed to MeJA for 0, 8, 20, and 40 h, and both total CQAs content and its fold changes were highest in G40h (19.59 mg g −1 , 17.65 times), followed by G20h (6.81 mg g −1 , 6.49 times), G8h (1.05 mg g −1 , 1.57 times), and G0h (0.711 mg g −1 , 1 time). Based on the results, it was observed that MeJA made a differential contribution to individual CQAs production ( Figures 1B-I).

RNA-seq and de novo Assembly
In order to comprehensively elucidate the time-course transcriptome of G. jasminoides cells under MeJA elicitation, four cDNA libraries from the cells stimulated by MeJA for 0, 8, 20, and 40 h were sequenced via the Illumina HiSeq × Ten platform, which generated 36. 94, 40.39, 40.92, and 38.96 million raw reads, respectively ( Table 1). After the quality control of raw reads, the clean reads were obtained with GC percentage ranging from 51.44 to 55.22% (Table 1). Subsequently, the Trinity program was used to assemble all clean reads into 63,724 transcripts with an N50 of 1,153 bp and an average length of 703 bp, which were then joined into 57,069 unigenes with an N50 of 1,051 bp and an average length of 655 bp ( Table 2). As shown in Figure 2A, all unigenes were longer than 200, and 200-300 bp was the most prominent length distribution interval. Among all unigenes, there were 35,782 (62.70%) unigenes with lengths shorter than 500 bp, 10,517 (18.50%) unigenes with lengths ranging from 500 to 1,000 bp, and 10,770 (18.88%) unigenes with lengths more than 1,000 bp. Overall, the quality of assembly results was high enough to conduct further analyses.

Functional Annotation and Classification of Unigenes
In order to investigate the functions of assembled unigenes, the entire unigenes were searched against the NR, Pfam, Swiss-Prot, GO, KOG, and KEGG databases. The annotation ratios of each database are shown in Figure 2B. A total of 46,054 unigenes (80.7%) were aligned to homologous sequences in those public databases.

Identification of DEGs During the MeJA-Treated Process
The pathway enrichment analysis was implemented by searching for DEGs against the KEGG database to get insight into mechanisms of biological functions and detect the genes responsible for cross-talk between jasmonate (JAs) signaling and key biochemical pathways. To identify DEGs and analyze their time-course expression profile during the MeJA-treated process, we analyzed the following comparisons: G0h vs G8h, G0h vs G20h, and G0h vs G40h, in which 3,611 (477 upregulated/3,134 downregulated), 3,405 (466 upregulated/2,939 downregulated), and 7,926 (860 upregulated/7,066 downregulated) DEGs were identified, respectively (Figures 3A-C).
The expression profiles based on the log 2 (fold change) values revealed a total of 8,264 DEGs in G0h vs G8h, G0h vs G20h, and G0h vs G40h ( Figure 3D). Furthermore, KEGG enrichment and cluster analyses on all DEGs were conducted to investigate the major biochemical metabolism and signal transduction pathways in which DEGs participated . All DEGs were individually enriched into 296 KEGG pathways, and the top 33 most significantly enriched pathways were displayed in Figure 3E. KEGG enrichment results showed that the top five most significantly enriched pathways were "Ribosome, " "Oxidative phosphorylation, " "Phenylpropanoid biosynthesis, " "Plant hormone signal transduction, " and "Steroid biosynthesis." Besides, other pathways involved in primary metabolism or secondary metabolism, e.g., "Fatty acid biosynthesis, " "Steroid hormone biosynthesis, " "Flavonoid biosynthesis, " and "Indole alkaloid biosynthesis, " were significantly enriched as well.
Generally, accumulation of the secondary metabolites was enhanced as a result of the increasing metabolic flux. Thus, the upregulated DEGs were subjected to the KEGG database for significant enrichment to further clarify the molecular mechanism of MeJA elicitation on the secondary metabolism (Figures 3F-H). Pathway enrichment on upregulated DEGs was significantly assigned to "Phenylpropanoid biosynthesis" and "Phenylalanine metabolism, " which were closely bound up with CQAs biosynthesis, in all comparisons. Other key biosynthesis pathways, including "Flavonoid biosynthesis, " "Terpenoid backbone biosynthesis, " "Sesquiterpenoid and triterpenoid biosynthesis, " "Indole alkaloid biosynthesis, " and "Isoquinoline alkaloid biosynthesis, " which provide precursors for CQAs biosynthesis or other therapeutic secondary metabolites, were significantly enriched as well. Moreover, the pathways in response to MeJA signaling (e.g., "alpha-Linolenic acid metabolism" and "Plant hormone signal transduction") were also shown to be highly enriched. Consequently, we make further efforts to focus on JAs biosynthesis and its signal transduction pathway, key precursor biosynthesis pathway, phenylpropanoid biosynthesis pathway, transporters, and TFs.
In terms of DEGs in the plant hormone signal transduction pathway (Supplementary Figure 1), a total of 87 DEGs were identified. Overall, the DEGs in the signal transduction [auxins (AUX), cytokinins (CK), brassinosteroids (BR), gibberellins (GA), abscisic acid (ABA), ethylene (ET), and salicylic acid (SA)] except for JA signaling showed a downward trend in G8h, G20h, and G40h compared to that in G0h (Supplementary Figure 2). As for JA signaling, six DEGs comprising MYC2 (three DEGs) and JAZ (three DEGs) showed a significant increase of transcript abundance from G0h and remained in extremely high expressions in G8h, G20h, and G40h. For instance, a MYC2 DEG (c10788_g1) and a JAZ DEG (c1385_g1) were most highly expressed (Figure 4C), in which the gene expression (FPKM value) in G40h could be 105.8 and 22.41 times those in G0h, respectively. Moreover, one DEG encoding Col 1 protein displayed a downward trend in G8h, G20h, and G40h compared to that in G0h. The results revealed that the exogenous MeJA treatment could modulate the expression of genes for endogenous JA biosynthesis and trigger the JAs signal transduction network.

Analysis of DEGs Engaged in Biosynthesis of Precursors for Phenylpropanoids
Aromatic amino acids, including phenylalanine, tyrosine, and tryptophan, were generated via the shikimate pathway in plants and microbes (Figure 5A), and phenylalanine was a vital precursor for phenylpropanoid biosynthesis. Our transcriptome data showed that 18 DEGs were annotated to the biosynthetic pathway for aromatic amino acids under MeJA treatment ( Figure 5B). Among them, there were 10 DEGs found to be upregulated and involved in phenylalanine and tyrosine biosynthesis, namely, 3-deoxy-7-phosphoheptulonate synthase (aroF/G/H, five DEGs), 3-dehydroquinate dehydratase I (aroD, one DEG), 3-phosphoshikimate 1-carboxyvinyltransferase (aroA, one DEG), arogenate/prephenate dehydratase (ADT/PDT, one DEG), and tyrosine aminotransferase (TAT, two DEGs). Another eight DEGs showed downregulated expressions, in which seven DEGs participated in tryptophan biosynthesis, namely, anthranilate synthase component II (TrpG, one DEG), anthranilate phosphoribosyltransferase (TrpD, two DEGs), tryptophan synthase beta chain (TrpB, three DEGs), and tryptophan synthase alpha chain (TrpA, one DEG). The above-mentioned results indicated that MeJA could stimulate the expression of genes for the phenylalanine and tyrosine biosynthesis and, meanwhile, downregulated the

Analysis of DEGs Related to Phenylpropanoid Biosynthesis and Transporter
Chlorogenic acids derived from phenylpropanoid biosynthesis pathway are a class of depsides of certain trans-cinnamic acids and quinic acids ( Figure 6A). A total of 19 DEGs (Figures 6A-C) were identified to be involved in CQAs biosynthesis, and 17 of these DEGs were found to be significantly upregulated in G8h, G20h, and G40h compared to those in G0h, namely, phenylalanine ammonia-lyase (PAL, nine DEGs), trans-cinnamate 4-monooxygenase (C4H, three DEGs), 4-coumarate-CoA ligase (4CL, two DEGs), shikimate O-hydroxycinnamoyltransferase (HCT, two DEGs), and coumaroyl quinate/coumaroyl shikimate 3 -monooxygenase (C3H, one DEG). In particular, the expression of the PAL family showed the most significant upregulation, and PAL (c19588_g1) was upregulated over 28.89-fold in G8h compared with that in G0h. However, none of the hydroxycinnamoyl CoA quinate hydroxycinnamoyl transferase (HQT) gene was even annotated by the BLAST program in our transcriptome data. These results implicated a correlation between the increasing metabolic flux and CQAs accumulation owing to MeJA elicitation.
In plants, the transporter system and subcellular localization are major elements in the process of biosynthesis regulation, transport, and accumulation of secondary metabolites (de Brito Francisco and Martinoia, 2018). The ATP-binding cassette (ABC) transporter and multidrug and toxin extrusion (MATE) protein were reported to be the main transporters for phenylpropanoid. In our transcriptome data (Figure 6D), 20 DEGs were subjected to ABC transporters, while nine DEGs were classified into MATE. Among these DEGs, four ABC transporters and four MATEs were upregulated in "G0h vs G8h, " "G0h vs G20h, " and FIGURE 5 | Putative biosynthesis pathways (A) and expression heatmap of DEGs (B) for aromatic amino acids (L-phenylalanine, L-tyrosine, and L-tryptophan); numbers in parentheses following each gene name indicate the number of corresponding DEGs. Some structural genes that were not DEGs were not shown in A.

Identification and Classification of DEGs Associated With TFs
Due to powerful regulation, TFs play a pivotal role in plant reproduction, bioactive metabolite synthesis, and response to environmental stress and thus have become a research hotspot recently (Sun et al., 2013). Further identification of TFs will help us gain a better understanding of gene regulatory networks, especially for CQAs biosynthesis in MeJA-treated G. jasminoides cells.
Some TF families (e.g., MYB, bHLH, NAC, ERF, and WRKY families), especially MYB, usually participate in the regulation on the secondary metabolism. Thus, DEGs in the above-mentioned TF families were classified for further analysis (Figure 7B).
Among these DEGs, there were 63 (9 up/54 down), 62 (12 up/50 down), and 147 (25 up/122 down) DEGs in the MYB family of "G0h vs G8h, " "G0h vs G20h, " and "G0h vs G40h, " respectively. Furthermore, a total of 38 DEGs, which could play positive/negative regulation on CQAs biosynthesis, were identified and listed in Figure 7C. In the positive regulators, the MYB85 family contained the maximum number of DEGs (17 unigenes), and nine of those DEGs displayed an upregulated expression in G8h, G20h, and G40h compared with those in G0h. Besides, MYB111 (one DEG), MYB20 (two DEGs), and MYB15 (one DEG) were found to be upregulated as well. In terms of the negative regulators (Figure 7C), 14 DEGs were identified, in which 11 downregulated DEGs were classified into MYB60 (one DEG), MYB52 (five DEGs), and MYB4 (five DEGs) families. To this end, these results revealed that MeJA could significantly modulate the expression of MYB TFs, which might participate in phenylpropanoid biosynthetic regulation.   that qRT-PCR expression patterns of the selected DEGs were highly correlated with those in RNA-seq. Thus, our RNA-seq data were reliable. Conclusively, this study proves that the MeJA treatment could trigger extensive transcriptional reprogramming and ultimately enhances CQAs accumulation.

Regulation of JAs Biosynthesis and
Signal Transduction in MeJA-Treated G. jasminoides Cells Jasmonic acid and its derivatives, collectively termed as JAs, are a class of potent lipid-derived phytohormones in plants that mediate responses to many biotic and abiotic stresses and perform various functions of reproduction and metabolic regulation via the signal transduction pathway (Kazan and Manners, 2008). Earlier studies indicated that JAs enhanced the formation of various phenylpropanoids (Prieto and Corchete, 2014;Liu et al., 2018). However, little is known about the specific role of JAs on CQAs accumulation in G. jasminoides cells. Alpha-linolenic acid originating from chloroplast membranes is considered as an important precursor for JAs biosynthesis. Undergoing condensation, cyclization, and β-oxidation reactions, JA were ultimately formed by LOX, AOS, AOC, OPR, ACX, ACCA1, and MFP2 (Figure 4A; Schaller and Stintzi, 2009). Most of these DEGs showed the upregulated expression in "G0h vs G8h, " "G0h vs G20h, " and "G0h vs G40h." Specially, AOC, which participates in the construction of the characteristic cyclopentanone ring structure of JAs (Jacobo-Velázquez et al., 2015), is pivotal for the JAs biosynthetic pathway and stimulates the accumulation of various secondary metabolites. Overexpression of AOC could stimulate endogenous JA accumulation and transcript levels of genes (PAL, C4H, 4CL, COMT, and CAD) related to phenylpropanoid biosynthesis in rice. Meanwhile, higher accumulation of phenolic acids, e.g., cinnamic acid, ferulic acid, and benzoic acid, was also found in the AOC-overexpressed transgenic line (Fang et al., 2016). Hence, DEGs in the JAs biosynthetic pathway could modulate endogenous JAs biosynthesis and further trigger the JA signal transduction network.
The bHLH TF MYC2 is a central regulator involved in the transcriptional regulation of JA-responsive gene expression in the JA signal transduction pathway (Figure 4B), and its expression is repressed by a complex repressor of JASMONATE ZIM DOMAIN (JAZ) proteins in the absence of a JA signal. Upon sensing of the external stimuli (e.g., the application of exogenous MeJA), the SCF COI1 -ubiquitin-proteasome pathway is activated and then leads to the degradation of the JAZ repressor, which paves the way for MYC2 to regulate the JA-responsive gene expression (Kazan and Manners, 2008). In our study, we found that six DEGs encoding putative MYC2 (three DEGs) and JAZ (three DEGs) were significantly upregulated by MeJA treatment (Figure 4C). However, DEG of Col 1 (c20119_g3) showed a downward trend in the three comparisons. It was reported that there were positive and negative feedback regulatory loops to regulate JA signaling. On the one hand, the external stimuli promoted the degradation of the JAZ repressor to activate transcription of genes encoding MYC2. Moreover, the external stimuli subsequently facilitate the transcription of genes encoding JAZ, together with the downregulation of Col 1, resulting in resetting the signal transduction pathway and avoiding a runaway response (Kazan and Manners, 2008;Jacobo-Velázquez et al., 2015). In addition, MYC2 was reported to be a regulator of phenylpropanoid biosynthesis. Salvianolic acid B, a phenolic acid derived from phenylpropanoid, was strongly induced in transgenic Salvia miltiorrhiza that overexpressed SmMYC2, and expressions of genes for the phenylpropanoid biosynthesis (e.g., PAL, C4H, and 4CL) were upregulated in the overexpression lines (Yang et al., 2017).
The complex network of communication among plant hormone signaling pathways is often referred to as hormone cross-talk, which is employed in many plant processes (Shigenaga et al., 2017). Beside the JA signaling pathway, other hormone signaling, such as AUX, CK, BR, GA, ABA, ET, and SA, were also responsive to MeJA treatment in G. jasminoides cells (Supplementary Figures 1, 2). Whether hormone cross-talk participates in regulation of CQAs remains to be validated in future studies.

Candidate Genes Involved in CQAs Biosynthesis and Effects of DEGs on CQAs Accumulation in MeJA-Treated G. jasminoides Cells
Chlorogenic acid and its derivatives are derived from the phenylpropanoid biosynthesis pathway in plants, and CQAs in G. jasminoides cells consisted of three classes, namely, mono-CQAs (Figures 1C-E), di-CQAs (Figures 1F,G), and its derivatives (Figures 1H,I). It is believed that there might be three distinct pathways for CQAs biosynthesis ( Figure 6A). In the initial step, phenylalanine is catalyzed by PAL and converted into cinnamic acid. Subsequently, CQAs are synthesized through three putative pathways, which are briefly described as follows: (1) generating CQAs via an activated intermediate (caffeoylglycoside); (2) synthesis of CQAs from caffeoyl-CoA and quinic acid by HQT; and (3) hydroxylation of p-coumaroyl-quinic acid by C3H to form CQAs (Villegas and Kojima, 1986;Moglia et al., 2016). The last two pathways for CQA biosynthesis were extensively studied in recent years. Moreover, Lallemand et al. (2012) revealed that the synthesis of di-CQAs could be mediated by HCT with mono-CQA and CoA, which was validated in vitro. However, the structural genes for di-CQA derivative biosynthesis are still unknown and require further research work ( Figure 6A).
In our RNA-seq data, a number of genes encoding PAL, C4H, 4CL, HCT, and C3H were annotated, while none of the sequence in the RNA-seq data was aligned to HQT. Thus, the primary biosynthetic pathway for CQAs in G. jasminoides cells was speculated as follows: the mono-CQAs biosynthesis might be mainly performed on the basis of pathway (3) mentioned above (involving PAL, C4H, 4CL, HCT, and C3H). Subsequently, di-CQAs were generated by HCT from mono-CQA and ultimately converted to di-CQA derivatives in the presence of several unknown structural genes.
A total of 19 DEGs were identified to encode enzymes involved in different steps of mono-CQA biosynthesis, such as PAL (nine DEGs), C4H (three DEGs), 4CL (three DEGs), HCT (three DEGs), and C3H (one DEG), 17 of which showed significantly upregulated expressions over time, and the increased expressions of structural genes were consistent with the mono-CQA accumulation in MeJA-treated G. jasminoides cells (Figures 1C-E). Furthermore, the continuously upregulated expression of C3H (C18639_g1 and C4961_g1) further converted mono-CQA into di-CQAs (Figures 1F,G). Although the specific structural genes for di-CQA derivatives were still unclear, the enhancement of the precursor supply (mono-and di-CQAs) could be one of the reasons for their increasing accumulation (Figures 1H,I) owing to MeJA elicitation.

DEGs Involved in Transporter
In plants, the transporter system played a critical role in the process of biosynthesis regulation, transport, and accumulation of the secondary metabolites (de Brito Francisco and Martinoia, 2018). At the cellular level, CQAs were synthesized primarily in chloroplasts and cytoplasm and finally transferred to the vacuole for long-term storage, in which the transport process was involved . The transport of the secondary metabolites in plants mainly comprised membrane transportermediated and vesicle-mediated mechanisms (de Brito Francisco and Martinoia, 2018;Li et al., 2019), etc. As for the membrane transporter-mediated mechanism, ABC and MATE families' transporter played roles in the sequestration of various phenolic compounds in the vacuole, such as lignin, flavonoid, anthocyanin, and their precursors. An ABC-type gene (MtABCG10) was responsible for the membrane translocation of p-coumaric acid, which was a key precursor for CQAs' and other phenylpropanoid component's biosynthesis, and silencing the expression of MtABCG10 resulted in a lower accumulation of phenylpropanoid-derived medicarpin and its precursors (Biała et al., 2017). Furthermore, the expression of ABC transporters could be induced by MeJA elicitation. In Silybum marianum, ABC-type transporters were proposed to be implicated in the secretion of flavonolignans, and MeJA treatment could stimulate the extracellular flavonolignan accumulation as well as the expression of ABC-type transporters (Prieto and Corchete, 2014). MATEs, some of which are JAs-responsive secondary metabolite transporters, usually act as transporters for phenolic and other compounds. In V. vinifera, a MATE-type gene (VvMATE2) was assigned as the putative proanthocyanidin transporter expressed during berry development (Pérez-Díaz et al., 2014), and MATE1 from Medicago truncatula was engaged in flavonoid transport (Zhao and Dixon, 2009). Recently, Li et al. (2019) reported a vesicle-mediated mechanism for the CQAs transport in L. japonica. Interestingly, some MATE-type transporters were required for the vesicle trafficking, and, for instance, the vesicle trafficking of anthocyanidin into vacuole in V. vinifera was mediated by two MATE-type transporters (anthoMATEs) (Gomez et al., 2011). A total of 29 DEGs were identified in the ABC and MATE superfamilies ( Figure 6D) in our RNA-seq data. The differential expression and putative roles of ABCs and MATEs revealed that they might be associated with response to MeJA and secondary metabolite accumulation, especially for CQAs, in G. jasminoides cells. Because information about the specific CQAs transport is still limited, the further functional characterization of the transporters needs to be made out. Thus, our results might provide some useful information concerning transport and accumulation of CQAs in G. jasminoides cells.

Regulation of DEGs Associated With TFs on CQAs Accumulation
Transcription factors, also known as sequence-specific DNA binding proteins, emerge as one of the key factors that modulate the expression of specific genes and the accumulation of messenger RNA at the transcriptional level through sequencespecific DNA binding and protein-protein interactions (Liu et al., 2015). A number of candidate DEGs of TFs were identified in MeJA-treated G. jasminoides cells (Figure 7B), some of which (e.g., MYB, bHLH, ERF, and WRKY families) could regulate JAs-induced accumulation of the secondary metabolites (Zhou and Memelink, 2016). TFs of MYB families were proposed to be involved in the positive/negative regulation on biosynthesis of various secondary metabolites, especially on phenylpropanoid biosynthesis (Dubos et al., 2010;Liu et al., 2015). The phenylpropanoid pathway is a complex metabolic network with many shared substrates and branches, and some MYB TFs can simultaneously modulate multiple compounds rather than a single compound, such as chlorogenic acids, flavonoids, and lignins, by regulating the expression of genes in the shared upstream pathway (Tuan et al., 2014).
In terms of the positive regulator from the MYB family ( Figure 7C), MYB85 and MYB20 TFs from Arabidopsis thaliana were developmentally associated with the secondary wall thickening via inducing the expression of 4CL and HCT, which catalyzed the common biosynthetic step of chlorogenic acids and lignin ( Figure 6A; Geng et al., 2020). It was reported that MYB15 could drive lignification through activating the expression of genes related to lignin biosynthesis, e.g., PAL, C4H, HCT, COMT, CCoAOMT, and CAD genes (Chezem et al., 2017). Among them, PAL, C4H, and HCT were the common upstream genes in chlorogenic acids and lignin biosynthesis pathways ( Figure 6A). In addition, MYB103 displayed a positive regulation on the F5H expression and S-lignin biosynthesis in A. thaliana (Ohman et al., 2012). Moreover, MYB111, a R2R3-MYB protein, specially controlled the early steps of the flavonoid biosynthetic pathway catalyzed by chalcone synthase (CHS), chalcone isomerase (CHI), flavanone 3-hydroxylase (F3H), and flavonol synthase (FLS), and its expression could stimulate flavonoid accumulation (Stracke et al., 2007).
In this work, 36 DEGs were identified and classified to several subgroups of the MYB TF family from A. thaliana, including four kinds of positive regulators (MYB85, MYB20, MYB15, and MYB111) and three kinds of negative regulator (MYB52, MYB4, and MYB6) on phenylpropanoid biosynthesis. Most of the positive regulators displayed a continuously upregulated expression, but most of the negative regulators showed a downtrend in expression after MeJA treatment ( Figure 7C), which might be one of the reasons for the upregulation of structural genes for CQAs and lignin biosynthesis (Figure 6). Though the functional characterization of those MYBs and other TFs needs to be verified in the future, the above findings might bring insight into the molecular mechanisms regulating CQAs contents in G. jasminoides cells under MeJA elicitation.

CONCLUSION
Our work indicated that the application of exogenous MeJA could effectively stimulate the CQAs accumulation in G. jasminoides cells. The stimulation mechanism of MeJA was further investigated from the perspective of differential expression of genes, including genes related to JAs biosynthesis, signal transduction, biosynthesis of precursor for CQAs, CQAs biosynthesis, transporters, and TFs. Our RNA-seq analysis of MeJA-mediated transcriptional changes indicated that numerous unigenes in these above-mentioned pathways were differentially expressed, which might be implicated in CQAs biosynthesis and regulation (Figure 9). In MeJA-treated G. jasminoides cells, MeJA triggers the expression of genes involved in endogenous JAs biosynthesis (LOX, AOS, AOC, OPR, ACX, and ACCA1) and JAs signal transduction (MYC2, JAZ, and Col 1). Then through the signal transduction network, genes related to biosynthesis of aromatic amino acids, namely, precursors for CQAs, showed an increased expression (aroF/G/H, aroD, aroA, ADT/PTD, TAT, and hisC), which might boost supply of the precursors. Meanwhile, TFs, especially for the MYB family, showed significant response to MeJA treatment and might display positive (MYB85, MYB20, MYB15, and MYB111) or negative (MYB52, MYB4, and MYB6) regulations on CQAs biosynthesis, resulting in increased expression of key genes for CQAs biosynthesis (PAL, C4H, 4CL, HCT, and C3H). Ultimately, the accumulation of CQAs in G. jasminoides cells was significantly increased. Moreover, our data will surely provide a massive genetic resource for further investigation of CQAs biosynthesis and lay the foundations for genetic engineering to boost the yield of these important compounds in G. jasminoides cells.

DATA AVAILABILITY STATEMENT
The data presented in the study are deposited in the (NCBI Sequence Read Archive) repository, accession number (PRJNA672865).

AUTHOR CONTRIBUTIONS
ZL: methodology, validation, investigation, data curation, visualization, writing-original draft, and writing-review and editing. AM: formal analysis and writing-review and editing. ZW: methodology, validation, and writing-review and editing. XZ: methodology and validation. YZ and LC: supervision and validation. MG: supervision, project administration, and writing-review and editing. ZY: supervision, funding acquisition, project administration, and writing-review and editing. All authors contributed to the article and approved the submitted version.