Transcriptome profiling shows gene regulation patterns in ginsenoside pathway in response to methyl jasmonate in Panax Quinquefolium adventitious root

Here, we combine elicitors and transcriptomics to investigate the inducible biosynthesis of the ginsenoside from the Panax quinquefolium. Treatment of P. quinquefolium adventitious root with methyl jasmonate (MJ) results in an increase in ginsenoside content (43.66 mg/g compared to 8.32 mg/g in control group). Therefore, we sequenced the transcriptome of native and MJ treated adventitious root in order to elucidate the key differentially expressed genes (DEGs) in the ginsenoside biosynthetic pathway. Through DEG analysis, we found that 5,759 unigenes were up-regulated and 6,389 unigenes down-regulated in response to MJ treatment. Several defense-related genes (48) were identified, participating in salicylic acid (SA), jasmonic acid (JA), nitric oxide (NO) and abscisic acid (ABA) signal pathway. Additionally, we mapped 72 unigenes to the ginsenoside biosynthetic pathway. Four cytochrome P450s (CYP450) were likely to catalyze hydroxylation at C-16 (c15743_g1, c39772_g1, c55422_g1) and C-30 (c52011_g1) of the triterpene backbone. UDP-xylose synthases (c52571_g3) was selected as the candidate, which was likely to involve in ginsenoside Rb3 biosynthesis.

Panax quinquefolium L., commonly known as American ginseng, belongs to the Araliaceae family and has gained tremendous global trade and recognition as a health food supplement. The dried root powder of this plant has been used extensively for its antitumor, anti-stress, anti-ageing, anti-fatigue, cardioprotective and hepatoprotective properties 1,2 . Ginsenoside are secondary metabolites of the P. quinquefolium and its major pharmacologically active components. Ginsenosides are divided into protopanaxadiol, protopanaxatriols and oleanolic acid based on their structure.
Our laboratory previously demonstrated that ginsenoside production is enhanced by the addition of methyl jasmonate (MJ) into Panax ginseng adventitious roots 13 . The elicitors can be recognized by plant receptors which are located on the surface of the plasma membrane or endomembrane. The receptors are activated, and then in turn activate their effectors, such as ion channels, GTP binding proteins (G-proteins), and protein kinases and oxidative burst 14 . Activated effectors can promote the synthesis of signaling molecules, such as salicylic acid (SA), jasmonic acid (JA), nitric oxide (NO), abscisic acid (ABA) and so on, which transfer the elicitor signals to defense genes that have been induced by elicitor treatment, and further amplify the elicitor signal to the biosynthesis of secondary metabolites 15 .
In order to elucidate the genes involved in the ginsenoside biosynthetic pathway, we sequenced, then compared two sets of transcriptome profiles derived from MJ treated and untreated P. quinquefolium adventitious roots. In addition, signaling molecules and defense genes in response to MJ were also studied. Our results provide a foundation for genetically constructing the ginsenoside biosynthetic pathways, which will in turn aid the study of its regulation and metabolic engineering of ginsenoside compounds.

Results and Discussion
Effects of MJ on active compound content. Treatment of adventitious root cultures with MJ stimulated ginsenoside accumulation (43.66 mg/g) compared with the control group (8.32 mg/g), decreasing root biomass and polysaccharide content (Fig. S1a-d). The biomass showed a negative correlation with electrical conductivity (EC) (Fig. S1b). Finally, these results led to higher ginsenoside productivity (105.74 mg/l) compared to the control group (30.19 mg/l) (Fig. S1f).
Treatment of P. quinquefolium adventitious root cultures with MJ resuled in an increase in ginsenoside content. MJ, a derivative of jasmonic acid, is an effective elicitor that is involved in plant defense response pathways and triggers plant metabolite biosynthesis. Therefore, MJ has been used for inducing metabolite production in plant cell cultures. MJ treatment is known to activate proteinase inhibitor genes in plants, which is likely to explain the observed decrease in biomass 16 . MJ-induced NO, SA, JA and ABA accumulation. The results of this study showed that MJ can induce NO, SA, JA and ABA accumulation in P. quinquefolium adventitious roots. As shown in Fig. 1, increase of NO, SA, JA and ABA were observed, reaching the highest level (732.44 μ mol·gprotein −1 , 0.08 ng·g −1 , 1.08 ng·g −1 and 21.11 ng·mL −1 ) at 24 h, respectively. NO is a key signal molecule in plant that induces a defense response to elicitors. It has been reported that NO can participate in the secondary metabolite accumulation such as ginsenosides 17 , taxanes 18 and other bioactive compounds. Hu et al. 17 found that NO was required for oligogalacturonic acid-induced saponin synthesis in cell cultures of P. ginseng. Furthermore, elicitors also induce the accumulation of SA and JA in plant cell. Xu et al. 19 found that SA and JA have synergistic effects on regulating elicitor-induced puerarin accumulation in cell culture. ABA acts as an important signal molecule to regulate biosynthesis of secondary metabolites in some plant cell. ABA can stimulate production of indole alkaloids in C. roseus cell culture 20 and taxol production in Taxus spp cell culture 21 . In this work, MJ induced the accumulation of signal molecules (NO, SA, JA and ABA) and enhanced the ginsenosides contents.
Functional annotation and gene ontology classification. RNA samples were extracted from control and MJ treated P. quinquefolium adventitious roots. Illumina RNA sequencing technology was used to sequence the whole transcriptome of P. quinquefolium. Unigenes with sequence orientation were aligned against public databases such as the Nr, SwissProt, Nt, Pfam, COG (Clusters of Orthologous Groups of proteins), GO (Gene Ontology) and KO (KEGG Ortholog).
Unigenes with Nr annotation were further annotated and classified under GO. GO is an international standardized gene functional classification system. Among the GO classifications, assignments to the biological process class ranked highest (286,854), followed by cellular component (82,491) and molecular function (56,853). Within the biological process category, the majority of the GO terms were assigned to cellular and metabolic processes. Within the cellular components category, transcripts assigned to cell and cell parts were the most common. For molecular function, the assignments were mostly binding and catalytic activity (Fig. 2). Transcripts related to GO term binding were most abundant in the molecular function category.
COG is a classification system based on orthologous genes. Orthologous genes have the same function and a common ancestor. We annotated 12,354 unigenes to 26 groups using the COG database (Fig. 3). The largest number (1,929) of the annotated unigenes fell within the general functional prediction only (R), while the fewest number of unigenes (1) were annotated as an unnamed protein (X). Additional assignments included 593 unigenes within the secondary metabolites biosynthesis, transport and catabolism (Q) category and 79 as defense mechanism unigenes (V) (Fig. 3).
For medicinal plants, RNA-seq has been used to identify genes that are directly or indirectly involved in the biosynthetic pathways of target bioactive compounds 22 . Transcriptome analysis from only one plant organ or tissue does not provide a full transcript catalogue, even though it can serve numerous specific genetic and breeding objectives. In our current study, we generated a total of 201 million high-quality reads from control and MJ treated adventitious roots, which is significantly more than the previously reported data of 454 ESTs 3 . Functional genomics studies require highly reliable reference sequences. Therefore, the transcript library we assembled here has significant implications for functional genomic studies on P. quinquefolium due to the high sequencing depth. Transcripts related to GO term binding were most abundant in the molecular function category. This is in agreement with the previously reported GO annotation of P. ginseng adventitious roots 23 .
Differentially expressed unigenes analysis. Unigene expression was calculated using a reads per kb per million reads (RPKM) method. A total of 12,148 DEGs were identified in the two experiments (control and MJ) including 5,759 up-regulated and 6,389 down-regulated genes.
GO functional analysis was also integrated with the clustering analysis of expression patterns. Within the DEGs group, metabolic process, cell and cell part genes were found to be abundant (Fig. S2a). Within the up-regulated DEGs, GO terms related to the biosynthetic process as well as intracellular and single-organism metabolic processes were significantly enriched (Fig. S2b). By contrast, genes primarily down-regulated after induction were related to binding and heterocyclic compound binding (Fig. S2c). Approximately 29% of the DEGs were categorized as genes responsive to stimuli and stress (Fig. S2a), suggesting that the expression of a large number of genes may be altered in response to an external stimulus 24 .
Pathway assignment for all transcripts was performed based on the KEGG database. In the case of Panax, well-represented pathways included amino acid metabolism, carbohydrate metabolism, lipid metabolism 3,25 and energy metabolism 26 . In this study, MJ were chosen as exogenous precursors to increase ginsenoside production in P. quiquefolium adventitious roots. Most primary metabolic processes such as the citrate cycle, carbohydrate metabolism, and amino acid metabolism were significantly up-regulated pathways among the DGEs. Primary metabolism is essential for plant growth, plant development and plant reproduction. In cell suspension cultures, primary metabolism is essential for plant cells to propagate in liquid media 22 . Analysis of defense genes. Across all samples, several defense-related genes were identified, participating in SA, JA, NO and ABA signal pathway. pathogenesis-related protein 1 (PR1), allene oxide cyclase (AOC), 9-cis-epoxycarotenoid dioxygenase (NCED), ABA responsive element binding factor (ABF) and zeaxanthin epoxidase (ZEP) showed a pronounced up-regulation by MJ elicitation (Fig. 5).
NPR1 acts as a receptor of the SA signal and then activities PR genes expression, including PR1 27 . On the other hand, the JA signaling pathway is positively regulated by the nuclearlocalized helix-loop-helix-leucine zipper-type transcription factor MYC2 and induces plant defense related proteins, such as AOC and LOX 28 . Recent studies showed that ABA involved in a complicated network of synergistic and antagonistic interactions with other phytohormones. Signaling related genes that are modified by ABA include NCED, ABF, ZEP and sucrose nonfermenting 1-related protein kinase 2 29 . Besides, we have not identify NO-related defense genes that up-regulated significantly. In this study, the expression of PR1, AOC, NCED, ABF and ZEP up-regulated significantly, which indicated that the signal molecules, generated after MJ treatment, enhanced ginsenosides content by regulating defense genes, consisting with the results of previous studies. Cerato-platanin triggers SA-signaling pathways, as revealed by the expression of PR genes and induced the biosynthesis of camalexin 30 .
Analysis of genes involved in ginsenoside biosynthesis. Using our RNA-seq data, we inspected the expression of genes from the upstream triterpenoid precursor biosynthetic pathways, named the cytosolic mevalonic acid (MVA) pathway and the plastidial 1-deoxy-D-xylulose-5-phosphate pathway.
We observed an up-regulation of genes involved in mono-, sesqui-and tri-terpenoid metabolism. These data are in agreement with previous studies that reported transcriptional up-regulation of the precursor pathways that likely increase synthesis of terpenoid natural products 31 . Mechanistically, MJ binds to membrane receptors and activates G-proteins to trigger phospholipase A (PLA). Subsequently, PLA activates α -linolenic acid and endo-methyl jasmonate. Endogenous MJ regulates the HMGR pathway 32 and down-stream genes to produce the mono-, di-, sesqui-, and tri-terpenoid genes 23 . From our investigation, we also found that the most abundant up-regulated unigenes were assigned to HMGR and GPS in response to MJ. Regulation of the cycloartenol synthase (CAS) and lupeol synthase (LS) genes, leading to production of phytosterols and lupeol, did not significantly change. These data suggested that treatment of adventitious roots with MJ resulted in the attenuation of competitive pathways and eventually diversion of the metabolic flux to the production of the desired ginsenosides.

Analysis of putative genes involved in the late steps of ginsenoside biosynthesis.
As one of the best-characterized protein families, CYP450s are known to catalyze the oxidation function of carbon-carbon bonds as well as alkyl hydroxylation and hydroxyl oxidation reactions 33 . Our RNA-seq data revealed 7 CYP450s (c52011_g1, c48642_g1, c15743_g1, c39772_g1, c38567_g1, c35627_g1, c55422_g1) that likely involved in ginsenoside biosynthesis. GTs are another large multigene family in plants. In this study, a total of 5 GTs unique sequences (c52571_g3, c45579_g2, c47755_g1, c39632_g1, c51194_g1) were found and likely to be involved in ginsenoside biosynthesis. Thus, 7 CYP450s and 5 GTs were selected.
The phylogenetic relationship between the 7 fulllength CYP450s of P. quinquefolium adventitious root and characterized CYP450s from other plants was depicted in Fig. 6. It is noteworthy that c15743_g1, c39772_g1, c38567_g1, c35627_g1, c55422_g1 were phylogenetically close to CYP88D6, a β -amyrin 11-oxidase from G. uralensis 34 . Phylogenetic analysis also found that the obtained full-length of c52011_g1 was close to CYP72A154, a β -amyrin 30-oxidase from G. uralensis 35 .
Phylogenetic analysis showed the relationship of P. quinquefolium adventitious root GT sequences to other functionally characterized members of plant GT families (Fig. 7). Among them, alpha-1,3-glucosyltransferase (c47755_g1) and UDP-xylose synthases (c52571_g3) were regarded as a lead candidate GTs responsible for triterpene saponin biosynthesis, because of its close relation to triterpene glucosyltransferases UGTPg100, UGTPg101, UGT74AE2. Hence, there were 6 candidate CYP450 unigenes and 2 candidate GT unigenes.
Quantitative PCR analysis was performed on 8 selected CYP450 and GT genes putatively involved in the ginsenoside biosynthesis of P. quinquefolium. The qPCR results of 8 selected genes showed general agreement with their transcript abundance changes as determined by RNA-seq (Fig. 8). Four P450 genes (c15743_g1, c52011_ g1, c39772_g1, c55422_g1), and one UDP-xylose synthases (c52571_g3) showed a significant up-regulation in response to MJ.
CYP450 and GT enzymes are critical for the downstream metabolism of ginsenosides to produce protopanaxadiol and protopanaxatriol. CYP450 monooxygenases play a key role in terpenoid biosynthesis, with such activity almost invariably required for further transformation of olefinic intermediates 31 . So far, the P450 compendium that can oxidize the dammarane and amyrin backbone is expanded, now covering six positions on the triterpene backbone and including C-6 (CYP716A53v2) 8 , C-11 (CYP88D6) 34 , C-12 (CYP716A47) 8 , C-16 (CYP716Y1) 36 , C-28 (CYP716A52v2, CYP716A12) 9 and C-30 (CYP72A154) 35 . In this study, phylogenetic and PCR analysis found that the c15743_g1, c39772_g1, c55422_g1 were close to CYP88D6 and c52011_g1 was close to CYP72A154.

Effects of MJ on expression of functional genes.
Genes expression level (GPS, FPS, SS, SE, β -AS, DS, CYP716A47, CYP716A53v2, UGT74AE2, UGT94Q2, UGTPg100 and c52571-g3) in adventitious roots after MJ treatment in different time (0 h, 12 h, 24 h and 48 h) were studied. The expression levels of functional genes were up-regulated compared with untreated group. In particularly, the expression levels of UGT74AE2, UGT94Q2 and UGTPg100 that generate Rh 2 , Rg 3 and Rh 1 respectively, reaching its peak at 24 h, 12 h and 12 h respectively, which consisted with production of monomer ginsenoside Rh 2 , Rg 3 , and Rh 1 . Besides, the expression level of c52571-g3, a candidate gene responsible for Rb 3 , was also up-regulated (Fig. 9).
In Artemisia annua, MJ induced artemisinin biosynthesis by up-regulating the expression of the genes involved in artemisinin biosynthesis 40 . This study showed that exposure to MJ in adventitious roots of P. quinquefolium enhanced the production of ginsenosides through regulated the expression of functional genes involved in triterpene biosynthesis.
In summary, we report here a validated large-scale transcriptome data set of P. quinquefolium adventitious roots. This study provides an important resource for understanding the formation and accumulation of secondary metabolites, paving the way for industrialization of ginsenosides.   After pre-cultivation for four weeks, 5.0 mg L −1 MJ was added to the P. quinquefolium adventitious roots medium. The roots were taken at 0, 12, 24 or 48 h to determine signal molecule (SA, JA, NO and ABA), ginsenoside content and expression level of functional genes. Culture conditions were the same as above. Each experiment was repeated at least three times.

Methods
Library preparation and sequencing. Total RNA of control and MJ-treated adventitious roots for 12 h were isolated using the Plant RNA KitII(OMEGA, USA). RNA integrity was assessed using the RNA Nano 6000 Assay Kit of the Agilent Bioanalyzer 2100 system (Agilent Technologies, CA, USA). A total amount of 3 μ g RNA per sample was used as input material for the RNA sample preparations. Sequencing libraries were generated using NEBNext ® Ultra ™ RNA Library Prep Kit for Illumina ® (NEB, USA) following manufacturer's recommendations. Index codes were added to each sample to identify attributes for each sequence. The clustering of the index-coded samples was performed on a cBot Cluster Generation System using TruSeq PE Cluster Kit v3-cBot-HS (Illumia) according to the manufacturer's instructions. After cluster generation, the library preparations were sequenced on an Illumina Hiseq 2500 platform and paired-end reads were generated.
Transcript assembly and annotation. The left files (read1 files) from all libraries/samples were pooled and labeled as left.fq file. Similarly, the right files (read2 files) were pooled and labeled as right.fq file. Transcriptome assembly was accomplished based on the left.fq and right.fq files using Trinity 42 with min_kmer_ cov set to 2 and all other parameters left at default settings.
Gene function was annotated based on the following databases: Nr (NCBI non-redundant protein sequences); Nt (NCBI non-redundant nucleotide sequences); Pfam (Protein family); COG (Clusters of Orthologous Groups of proteins); Swiss-Prot (A manually annotated and reviewed protein sequence database); KO (Kyoto Encyclopedia of Genes and Genomes (KEGG) Ortholog database); GO (Gene Ontology). Differential expression analysis. Differential expression analysis of two conditions was performed using the DESeq R package (1.10.1). DESeq provide statistical routines for determining differential expression in digital gene expression data using a model based on the negative binomial distribution. The resulting values were adjusted using the Benjaminiand Hochberg's approach for controlling the false discovery rate. Genes with an adjusted P-value < 0.05 found by DESeq were assigned as differentially expressed.

GO and KEGG enrichment analysis. GO enrichment analysis of the differentially expressed genes (DEGs)
was implemented by the GOseq R packages, basing on Wallenius non-central hyper-geometric distribution 43 , which can adjust for gene length bias in DEGs. We used KOBAS 44 software to test the statistical enrichment of DGEs in KEGG pathways.
Quantitation of SA, JA, NO and ABA. 0.25 g of adventitious roots (control group and MJ group) were ground into powder using a mortar and pestle chilled with liquid nitrogen. Extraction and analyses of SA and JA were performed as described previously 45 . The extracts of NO and ABA were prepared by homogenizing 0.2 g of adventitious roots (control group and MJ group) in a mortar on ice, using 1.0 mL distilled water. The contents of NO and ABA were measured using commercially available kits (Nanjing Jiancheng Bioengineering Research Institute, Nanjing, China) according to the manufacturer's instructions.
Quantitative PCR. For each qRT-PCR reaction, 200 ng of total RNA was used for first strand cDNA synthesis. First-strand cDNAs were used as a template for RT-PCR reactions, which were performed as follows: 94 °C for 2 min, then 35 cycles of 94 °C for 30 s, 57 °C for 1 min, and 72 °C for 50 s; with a final 2 min extension at 72 °C. We used the ABI7500 for quantitative PCR reactions and the relative standard curve method was adopted to analyze the relative expression of genes. The PCR products were determined by agarose gel (2%) electrophoresis. The size of the fragments was estimated using a 100-bp ladder (CWBIO, China) as a size marker. All experiments were performed in triplicate. The primers of genes used in RT-PCR are shown in Table S1, Table S2. HPLC analysis. Samples were analyzed using an Agilent HPLC system containing a surveyor autosampler.
The details of the analytical procedures have been previously described 41,46 . Phylogenetic analysis. Probable entire amino acid sequences of CYP7450s and UGTs were taken from the GenBank database (http://www.ncbi.nlm.nih.gov) and evolutionary distances were computed using the Poisson correction method, and a Neighbor-Joining (NJ) tree was constructed with MEGA4. The indicated scale represents 0.1 amino acid substitutions per site. Bootstrap values obtained after 1000 replications are given on the branches.