Integrated Transcriptomic and Metabolomic Analysis of Five Panax ginseng Cultivars Reveals the Dynamics of Ginsenoside Biosynthesis

Panax ginseng C.A. Meyer is a traditional medicinal herb that produces bioactive compounds such as ginsenosides. Here, we investigated the diversity of ginsenosides and related genes among five genetically fixed inbred ginseng cultivars (Chunpoong [CP], Cheongsun [CS], Gopoong [GO], Sunhyang [SH], and Sunun [SU]). To focus on the genetic diversity related to ginsenoside biosynthesis, we utilized in vitro cultured adventitious roots from the five cultivars grown under controlled environmental conditions. PCA loading plots based on secondary metabolite composition classified the five cultivars into three groups. We selected three cultivars (CS, SH, and SU) to represent the three groups and conducted further transcriptome and gas chromatography-mass spectrometry analyses to identify genes and intermediates corresponding to the variation in ginsenosides among cultivars. We quantified ginsenoside contents from the three cultivars. SH had more than 12 times the total ginsenoside content of CS, with especially large differences in the levels of panaxadiol-type ginsenosides. The expression levels of genes encoding squalene epoxidase (SQE) and dammarenediol synthase (DDS) were also significantly lower in CS than SH and SU, which is consistent with the low levels of ginsenoside produced in this cultivar. Methyl jasmonate (MeJA) treatment increased the levels of panaxadiol-type ginsenosides up to 4-, 13-, and 31-fold in SH, SU, and CS, respectively. MeJA treatment also greatly increased the quantity of major intermediates and the expression of the underlying genes in the ginsenoside biosynthesis pathway; these intermediates included squalene, 2,3-oxidosqualene, and dammarenediol II, especially in CS, which had the lowest ginsenoside content under normal culture conditions. We conclude that SQE and DDS are the most important genetic factors for ginsenoside biosynthesis with diversity among ginseng cultivars.


INTRODUCTION
Panax ginseng is a medicinal plant that has long been used as a tonic agent throughout Asia, including Korea and China (Jia et al., 2009). Ginseng is rich in useful metabolites, such as polyacetylenes, polysaccharides, phenolic compounds, and terpenes. Triterpene saponins, also known as ginsenosides, are the major compounds in ginseng showing anti-oxidant, antiinflammatory, and anti-cancer activity (Kim, 2012;Zhang et al., 2013). Dozens of ginseng accessions have been bred through pure line selection from three local landrace populations (Jakyung, Chungkyung, and Hwangsook) and registered as cultivated varieties (cultivars) in Korea . These cultivars exhibit different morphological and physiological characters Lee et al., 2015), as well as different metabolite accumulation patterns Kim et al., 2009;Lee et al., 2011;Cho et al., 2012).
Candidate genes responsible for ginsenoside biosynthesis in the Panax genus have been identified in P. ginseng Jung et al., 2014;Kim et al., 2014;Jayakodi et al., 2015;Wang P. et al., 2015;Wei et al., 2015), P. notoginseng Rai et al., 2016), and P. vietnamensis (Zhang et al., 2015). Transcriptome analysis of P. notoginseng identified several candidate unigenes encoding cytochrome P450 and glycosyltransferase and revealed that these genes were conserved across all Panax species (Rai et al., 2016). However, the regulatory mechanism for the ginsenoside biosynthesis pathway still remains unknown. Integrated omics analysis of the transcriptome and metabolome is an efficient tool for obtaining a precise understanding of biosynthetic pathways. This technique allows identification of novel functional genes and characterization of the regulatory mechanism and key factors of biosynthetic pathways (Saito et al., 2008;Mounet et al., 2009). The regulatory genes in the tryptophan (Dubouzet et al., 2007), flavonoid (Cho et al., 2016), flavanol , and glucosinolate biosynthetic pathways (Hirai et al., 2004(Hirai et al., , 2007 have been explored using integrated analysis, and potential target genes have been identified. In this study, we applied an integrated metabolomics and transcriptomics approach to explore the ginsenoside biosynthetic pathway using adventitious roots from five ginseng cultivars. Specifically, we compared the expression levels of genes involved in ginsenoside biosynthesis, the accumulation patterns of intermediates, and the composition and quantity of ginsenosides in adventitious roots with or without methyl jasmonate (MeJA) treatment. This analysis allowed us to uncover the dynamic regulatory mechanism for the ginsenoside biosynthetic pathway controlled by genetic factors and MeJA treatment. , and the roots were cultured as previously described (Jayakodi et al., 2014). Adventitious roots (42 days old) were treated with 200 µM MeJA. For metabolites profiling, the adventitious roots were harvested after 1 week of MeJA treatment, with five biological replicates for both nontreated and MeJA-treated adventitious roots. For quantification of gene expression using RT-PCR, adventitious roots were collected at 0, 12, 24, and 48 h after MeJA treatment, with three biological replicates. For transcriptome analysis, adventitious roots were harvested from plants without MeJA treatment, with three biological replicates. The adventitious roots were immediately stored at −80 • C until use.

Sample Preparation
For UPLC-ESI/MS analysis, 10 mg of lyophilized adventitious root tissue was sonicated with 1 mL of 70% methanol for 30 min at room temperature. The extracts were centrifuged at 13,000 × g for 5 min, and the supernatants were filtered through a 0.2 µm PTFE syringe filter (Toyo Roshi Kaisha, Japan). For LC-UV/DAD analysis, 10 mg of lyophilized adventitious root tissue was sonicated three times using 2 mL of 70% methanol for 30 min at room temperature. The extracts were dried using N 2 gas and resuspended in 200 µL of 80% methanol.
The micrOTOF-QII was run in negative mode, and the source conditions were as follows: capillary −4.5 kV, nebulizer pressure at 1.2 bar, dry gas flow rate of 8 L/min, and dry gas temperature at 200 • C. The ion transfer and collision stages were set as follows: funnel 1 RF 400 Vpp, funnel 2 RF 400 Vpp, hexapole RF 400 Vpp, quadrupole ion energy 15 eV, collision energy 10 eV, collision RF 400 Vpp, transfer time 100 µs, and pre-pulse storage 5 µs. Highpurity nitrogen was used as the nebulizer and dry gas, and argon was used as the collision gas. Lithium formate was used as an internal standard, with a scan range from 50 to 1500 m/z.

LC-UV/DAD Analysis
The nine ginsenosides [Rg 1 , Re, Rf, Rg 2 (20S), Rg 2 (20R), Rb 1 , Rc, Rb 2 , and Rd] were quantified using an Agilent 1260 infinity HPLC system (Agilent, United States) with a discovery C18 column (2.1 mm × 100 mm, 5 µm, Sigma-Aldrich, United States). The mobile phase was composed of A: water + 0.1% formic acid and B: acetonitrile + 0.1% formic acid, as used for UPLC-ESI/MS. The flow rate was 0.3 mL/min. The gradient conditions were as follows: 80% A held for 20 min, followed by 70% A for 35 min, 55% A for 45 min, 10% A for 55 min, holding for 5 min, and a post-run with 80% A for 15 min. The UV detector was set at 203 nm, and the injection volume was 5 µL. For quantification, calibration curves of each ginsenoside were obtained, and intra-and interday precision and accuracy parameters were determined for validation (Supplementary Tables S1, S2).

GC-MS Analysis
The squalene, 2,3-oxidosqualene, farnesol, and dammarenediol II contents in CS, SH, and SU were analyzed by GC-MS as described previously Lee et al., 2012). Briefly, KOH saponification-based extracts were analyzed using GCMS-QP2010 (Shimadzu, Japan) with a DB-5MS column (30 m × 0.25 mm, 0.25 µm, Agilent technologies, Germany). The carrier gas was helium, the flow rate was 1.0 mL/min, and the injection volume was 1 µL, in split mode (1:2). The mass detector was in electron impact mode of 70 eV, and the ion source temperature was 200 • C. The mass detection range was 40-600 m/z. A single peak of each compound was identified using selected single ion monitoring mode and the retention times of authentic compounds. The areas of the compounds were normalized to that of the internal standard.

Data Processing and Multivariate Analysis
Raw data generated from UPLC-ESI/MS were processed using MZmine software version 2.10 1 . The data were subjected to peak deconvolution using a Savitzky-Golay filter and aligned through the RANdom Sample Consensus algorithm. The relative abundance of each peak was centered and scaled to unit variance, with base weight set at 1/standard deviation. The filtered peak lists were imported into SIMCA-P+ (version 12.0, Umetrics, Sweden), and principal component analysis (PCA) and partial least square-discriminant analysis (PLS-DA) were performed. In the PLS-DA model, metabolites with variable influence on the projection (VIP) value above 1 and P-value under 0.05 were selected for the VIP lists.

RNA Isolation and Transcriptome Sequencing Using the Illumina Platform
Total RNA was extracted from adventitious roots of cultivar SU using an RNeasy Plant Mini Kit (Qiagen, Germany) according to the manufacturer's protocol. The quality and quantity of extracted RNA was checked using an ND-1000 (NanoDrop Technologies Inc., United States) and formaldehyde-agarose gel electrophoresis. An RNA-Seq library with a 300-bp insertion size was then constructed and sequenced on the Illumina sequencing platform (NextSeq 500, LabGenomics Co., South Korea). The three replicated RNA-Seq reads were deposited in the NCBI Sequence Read Archive 2 (SRA) under accession numbers SRR2134236, SRR2134277, and SRR2134367. RNA-Seq data for adventitious roots from SH and CS were obtained from a previous study (Lee Y.S. et al., 2016). 1 http://mzmine.sourceforge.net/ 2 http://www.ncbi.nlm.nih.gov/Traces/sra

Digital Expression Analysis of Genes Related to Ginsenoside Biosynthesis
For expression profiling among cultivars, candidate genes involved in ginsenoside biosynthesis were selected based on a previous report (Jayakodi et al., 2015). The expression levels of genes related to the ginsenoside biosynthetic pathway were determined based on fragments per kilobase of exon per million fragments (FPKM) values from RNA-seq data for CS, SH, and SU using RSEM (Li and Dewey, 2011). Genes showing significantly different expression among the three cultivars were identified by analysis of variance (ANOVA, P-value < 0.05) using R (version 3.1.0 3 ) and the R packages mvtnorm (Genz et al., 2008) and multcomp (Hothorn et al., 2008). Mean FPKM values of the genes were converted to log2 values, which were used for visualization using the heatmap package in R. An integrative analysis between expression of the unigenes and the relative quantity of 21 metabolites was conducted with a bidirectional multivariate regression method based on orthogonal projection to latent structures (OPLS) (Bylesjö et al., 2007). The regression method of gene expression profiles (mean-centered, x variables) and metabolite profiles (mean-centered and Pareto-scaled, y variables) was analyzed using O2PLS model implemented in SIMCA tool (version 14.1, trial, Umetrics, Sweden).

RT-PCR Analysis
For RT-PCR analysis, cDNA was produced from total RNA of CS, SH, and SU using a SMART cDNA Synthesis kit (Clontech, Japan) according to the manufacturer's protocol. The synthesized cDNA was diluted 1/10 and used as template for RT-PCR. The RT-PCR was performed on a LightCycler 480 (Roche, Germany) under the following thermal cycling conditions: 95 • C for 5 min, followed by 40 cycles of 95 • C for 15 s, 58 • C for 10 s, and 72 • C for 10 s. Primer sets for each gene were designed using Primer3 (Rozen and Skaletsky, 1999), and the specificity of the primer sets was confirmed by sequencing PCR amplicons on an ABI 3730 XL DNA Analyzer (Applied Biosystems, United States). All primer sets used in this study are listed in Supplementary Table S3.

Secondary Metabolite Profiles in Adventitious Roots of Five Ginseng Cultivars
We performed secondary metabolite profiling in adventitious root tissue from five ginseng cultivars (CP, CS, GO, SH, and SU) using UPLC-ESI/MS to investigate metabolic diversity among ginseng cultivars (Figure 2). A previously reported method using a reverse-phase column (Lee G.J. et al., 2016;Wang et al., 2016) was applied to selectively separate secondary metabolites including ginsenosides from adventitious root extracts. A total of 21 metabolites were identified via comparisons with authentic standards or data found in previous reports (Xie et al., 2008;Qi et al., 2012;MacCrehan and White, 2013). Among these, 16 metabolites with VIP FIGURE 2 | Representative chromatograph of adventitious roots from five ginseng cultivars (CP, CS, GO, SH, and SU) analyzed by UPLC-ESI/MS in negative mode; CP, Chunpoong; CS, Cheongsun; GO, Gopoong; SH, Sunhyang; and SU, Sunun. The identity of each peak is shown in Table 1. values above 1 and P-value below 0.05 were characterized as biomarkers able to discriminate among adventitious roots from different ginseng cultivars. These biomarkers were PPD/PPTtype ginsenosides, except for four metabolites identified as chikusetsusaponin IVa, gypenoside XVII, calenduloside E, and an isomer of ginsenoside Ro (Table 1). These results indicate that PPD/PPT ginsenosides are the main contributors to the metabolic differences in adventitious roots among these cultivars.
We subjected the UPLC-ESI/MS data to multivariate analysis to visualize the metabolic differences among cultivars. The five ginseng cultivars were divided into three groups according to the PCA loading plots (R 2 X: 0.809 and d Q 2 : 0.825): Group 1: CS; Group 2: SH; and Group 3: SU, CP, and GO. This grouping indicated that the accumulation pattern of metabolites varied among cultivars, but that some cultivars shared similar metabolic profiles (Figure 3).
UPLC-ESI/MS analysis also revealed that the relative abundance of representative ginsenosides in adventitious roots differed among the five ginseng cultivars (Figure 4). SH roots were rich in PPD-type ginsenosides such as Rb 1 , Rb 2 , Rc, and Rd compared to the other cultivars, although the amount of PPT-type ginsenoside in SH was comparable to that in other cultivars except CS. CS showed the lowest amounts of both PPD and PPT ginsenosides, although it was relatively rich in Rg 1 . The abundances of PPD-type and PPT-type ginsenosides were not significantly different among cultivars CP, GO, and SU.
Transcriptome Analysis of Adventitious Roots from CS, SH, and SU UPLC-ESI/MS analysis revealed differences in the accumulation patterns of metabolites among ginseng cultivars, suggesting that the ginsenoside biosynthetic pathway might be differentially regulated in these cultivars (Table 1 and Figure 3). We selected one cultivar to represent each group in the PCA loading plots (CS, SH, and SU; Figure 3) and compared these three cultivars via transcriptome analysis ( Table 2). We obtained transcriptome data for adventitious roots from CS and SH from a previous study (Lee Y.S. et al., 2016) and generated transcriptome data for SU in the current study. An average of 15 million (M) reads from three replicates were produced from SU adventitious roots. After trimming, an average of 14 M reads was subsequently used for expression profiling of genes involved in the ginsenoside biosynthetic pathway.

Comparison of Ginsenoside Biosynthesis Pathways in Adventitious Roots from CS, SH, and SU
We previously identified 23 candidate unigenes related to the MVA and triterpene biosynthetic pathways from the P. ginseng root transcriptome database (Jayakodi et al., 2015). To identify genes associated with the differences in ginsenoside contents among cultivars, we analyzed the expression levels of the 23     117.19 ± 10.06 * Significantly differentially expressed unigenes among adventitious roots from CS, SH, and SU are marked with an asterisk (P < 0.05), and significant differences are represented by different letters. CS, Cheongsun; SH, Sunhyang; SU, Sunun. The full names of genes are shown in Figure 1. unigenes based on FPKM values and compared the expression patterns among cultivars ( Table 3).
Among the seven candidate genes in the MVA pathway, unigene Pg_Root110920_c0_seq2 (encoding IDI) showed the highest expression while unigene Pg_Root123699_c0_seq7 (encoding HMGR) showed the lowest expression in adventitious roots from CS, SH, and SU. The expression levels of most unigenes in MVA biosynthesis were not significantly different among cultivars (Table 3). By contrast, the expression levels of unigenes involved in the triterpene biosynthesis pathway differed significantly ( Table 3). Unigene Pg_Root127352_c0_seq41 (encoding FPPS) was highly expressed in SU. Unigene Pg_Root104873_c0_seq1 (encoding SQE) and unigene Pg_Root126438_c1_seq1 (encoding DDS) were expressed at extremely low levels in CS compared to SH and SU. The expression level of each gene was correlated with the relative amounts of intermediates and total ginsenosides in adventitious roots of the three cultivars (Figure 4). Unigene Pg_Root123943_c0_seq1 (encoding PPDS) was the most highly expressed in SU whereas unigene Pg_Root91292_c0_seq1 (encoding PPTS) was highly expressed in all cultivars ( Table 3). Unigene Pg_Root120424_c0_seq15 (encoding β-amyrin synthase; β-AS) showed different expression among cultivars but its expression level was low compared to that of the DDS unigene (Pg_Root126438_c1_seq1).
The integration of gene expression with metabolite data by O2PLS model resulted in co-variation of 62% of transcripts with 80% of the metabolites. A pairwise Pearson correlations for those transcripts and metabolites showed a strong correlation (r > 0.7) with more than seven metabolites between unigene Pg_Root104873_c0_seq1 (encoding SQE) and Pg_Root126438_c1_seq1 (encoding DDS) (Supplementary Figure S1).

Ginsenoside Contents Are Altered after MeJA Treatment in Adventitious Roots from CS, SH, and SU
Methyl jasmonate treatment upregulates ginsenoside biosynthesis in P. ginseng (Yendo et al., 2010). We investigated the effects of MeJA treatment on ginsenoside biosynthesis in different ginseng cultivars to understand the underlying genetic factors ( Table 4). LC-UV-DAD analysis revealed that the three representative ginseng cultivars had different responses to MeJA treatment. In all three cultivars, the contents of PPD-type ginsenosides under MeJA treatment were higher than those of PPT-type ginsenosides. After 7 days of MeJA treatment, PPDtype ginsenoside levels increased 4.14-(SH), 13.15-(SU), and 31.44-(CS) fold, whereas PPT-type ginsenoside levels increased 1.34-(SH), 4.07-(SU), and 4.00-(CS) fold (Table 4).
Ginsenoside levels responded most strongly to MeJA treatment in CS, which had the lowest ginsenoside levels of the three cultivars under normal conditions, as well as the lowest expression of ginsenoside biosynthesis genes. In MeJA-treated CS adventitious roots, levels of Rd, a PPD-type ginsenoside, increased 90.32-fold, while the levels of other PPD-type ginsenosides, such as Rb 1 , Rc, and Rb 2 , increased 17.38-, 30.43-, and 32.92-fold, respectively. SH, the most ginsenoside-rich cultivar under normal conditions, was less responsive to MeJA treatment compared to CS and SU (Table 4).

Gene Expression and Quantity of Intermediates in Ginsenoside Biosynthesis under MeJA Treatment
To identify the genes and intermediates associated with changes in ginsenoside contents in the three cultivars following MeJA treatment, we analyzed the expression levels of genes and the contents of intermediates related to the triterpene biosynthetic pathway in MeJA-treated CS, SH, and SU adventitious roots (Figure 5).   Table 3 were compared among cultivars.
(B) Ginsenoside biosynthesis pathway and quantity of intermediates among three cultivars before (white bars) and after MeJA treatment (black bars).
Relative abundances of intermediates are represented in the bar graphs; significant differences are indicated by asterisks: * P < 0.05; * * P < 0.01 (t-test). (C) Expression (fold change) of each gene before and after MeJA treatment in adventitious roots of the three ginseng cultivars. Fold changes (Y -axis) in gene expression levels at 12 and 24 h (X-axis) compared to 0 h (analyzed by RT-PCR) are shown in the line graphs.
Gene expression analysis via RT-PCR revealed differences in the expression patterns of genes involved in the MVA pathway and the triterpene biosynthetic pathway (Figure 5 and Supplementary Table S4). Among genes in the MVA pathway, only AACT expression significantly increased after MeJA treatment (Supplementary Table S4). By contrast, most genes corresponding to the triterpene biosynthetic pathway were significantly induced by MeJA treatment, with different levels of increase detected among cultivars, indicating that the triterpene biosynthetic pathway could be responsible for the variations in ginsenoside accumulation among different genotypes of P. ginseng. In particular, the expression of DDS dramatically increased (17.62-to 71.83-fold) at 12 h post-elicitation in all cultivars, and SQE (7.9-to 15.97-fold) and SS (5.16-to 10.32-fold) were also significantly induced. Notably, the changes in expression of all three genes were highest in the cultivar with the lowest ginsenoside content, CS (Figure 5).
GC-MS analysis revealed that the relative abundance of intermediate compounds, from squalene to dammarenediol II, following MeJA treatment increased the most in CS ( Figure 5). Comparative analysis of adventitious roots among the three cultivars with or without MeJA treatment suggested that squalene to dammarenediol II might be the major limiting step that affects the variation in ginsenoside content among ginseng cultivars. Among these, DDS likely plays the most critical role in ginsenoside biosynthesis and accumulation.

Metabolic Profiles in Ginseng Adventitious Roots
Ginsenoside profiles and accumulation patterns can differ based on genotype Kim et al., 2009) and growth conditions (Matkowski, 2008). In this study, we focused on the genetic component of ginsenoside accumulation by using five genetically fixed inbred cultivars, CP, CS, GO, SH, and SU (Kwon et al., 2001;Ahn et al., 2008;Lee et al., 2008) and isolating ginsenosides from in vitro cultured adventitious roots grown under a controlled environment. We quantified the variation among cultivars in the contents of ginsenosides and intermediates in the ginsenoside biosynthesis pathway as well as the relevant differences in expression of the underlying genes.
We characterized 21 ginsenosides, 15 of which were PPD/PPT-type ginsenosides. Among these PPD/PPT-type ginsenosides, six metabolites were identified as malonylginsenosides (Table 1). Malonylation, a modification process that occurs widely in plants, takes place on the sugar moieties of phytohormones, xenobiotic compounds, and several secondary metabolites (Muth et al., 2008;Kochkin et al., 2013). Malonylation plays a role in improving the solubility of target compounds and facilitating their transport into storage organelles such as vacuoles (Taguchi et al., 2010). In the current study, malonyl-conjugated PPD-type ginsenosides accounted for five of the six identified malonyl-ginsenosides, which is consistent with a previous report that malonylation occurs more often in PPD-type ginsenosides than in PPT-type ginsenosides of Panax (Kochkin et al., 2013). Thus, some PPD-type ginsenosides that accumulate in adventitious roots are malonylated and, therefore, might be stored in organelles such as vacuoles.
The ginsenosides present in adventitious roots had different accumulation patterns in different cultivars. PCA loading plots revealed that that five ginseng cultivars could be divided into three groups according to their secondary metabolite profiles, although the PCA plot based on ginsenoside profile did not correlate precisely with that based on primary metabolite profiles (Lee Y.S. et al., 2016). Notably, however, the CS cultivar exhibited distinct metabolic accumulation in both the PCA loading plots based on primary metabolite profiles (Lee Y.S. et al., 2016) and those based on secondary metabolite profiles (Figure 3). These results indicate that the ginsenoside pathway is likely distinct in CS adventitious roots compared to other ginseng cultivars.

Dammarane-Type Ginsenoside Synthesis Rapidly Increases in Response to MeJA Treatment
Jasmonic acid and MeJA are signaling compounds related to plant defense (Derksen et al., 2013). They also have been used as important plant-derived elicitors for elevated production of secondary metabolites in various plant species (Ali et al., 2006;Han et al., 2006;Lee et al., 2013;Wang J. et al., 2015). Oxidosqualene, a precursor for dammarenediol synthase, is present at a branch point in the triterpene biosynthesis pathway. This compound is cyclized by various oxidosqualene cyclases, such as lanosterol synthase (LSS), cycloartenol synthase (CAS), β-AS, and DDS to produce phytosterols and ginsenosides (Thimmappa et al., 2014). LSS and CAS catalyze the conversion of oxidosqualene to lanosterol and cycloartenol, respectively, which is the first step for entry into phytosterol biosynthesis (Suzuki et al., 2006;Ohyama et al., 2009). The expression levels of LSS and CAS were constant or decreased in response to MeJA treatment in Centella asiatica (Bonfill et al., 2011) and also in P. ginseng (Han et al., 2006), suggesting that the phytosterols in ginseng play a fundamental role in plant growth and development rather than in defense mechanisms involving MeJA. In addition, the expression of β-AS, encoding an enzyme that produces β-amyrin from oxidosqualene during the first committed step of oleanane-type ginsenoside biosynthesis (Suzuki et al., 2006;Tansakul et al., 2006), is not significantly altered under MeJA treatment (Han et al., 2006). By contrast, we found that both the contents of dammarane-type ginsenosides (especially PPD-type) and the expression of DDS significantly increased after MeJA treatment ( Figure 5 and Table 4). This result suggests that the biological activities of dammaranetype ginsenosides are more closely involved in defense in ginseng than are those of oleanane-type ginsenosides. The differential regulation of dammarane-type ginsenoside biosynthesis in the ginseng cultivars we examined might influence the properties of these cultivars in terms of plant defense responses.

Dynamic Alteration of Ginsenoside Biosynthesis via MeJA Treatment
In the absence of MeJA treatment, the ginsenoside contents were considerably lower in adventitious roots of CS as compared to the other cultivars (Table 4 and Figure 4), which correlated well with the lower expression in CS of genes that function in the pathway from squalene to dammarenediol II (Table 3). In addition, the ginsenoside contents significantly increased in CS adventitious roots upon MeJA treatment; PPD-type ginsenoside levels were increased 17-to 90-fold, and the expression levels of genes such as SQS, SQE, and DDS, as well as the contents of intermediates such as squalene, 2,3oxidosqualene, and dammarenediol II, were highly elevated ( Figure 5). Remarkably, the expression level of DDS was upregulated 71-fold in CS in response to MeJA treatment and the relative dammarenediol II level in CS adventitious roots became comparable to those in SH and SU (Figure 5). In the case of CS, our data indicate that the low ginsenoside contents in adventitious roots are likely related to the low expression levels of SQE and DDS. The factors that underlie the low expression of these genes in CS are yet not known, but might be directly related to differences in function of the genes, such as unknown transcription factors of DDS.
The dramatic increase in DDS expression under MeJA treatment coincided with the rapid increase in ginsenoside contents, especially in CS adventitious roots. Indeed, the content of dammarane-type ginsenosides is decreased by 84.5% in the roots of transgenic ginseng plants with downregulated DDS expression (Han et al., 2006). These findings suggest that the route from squalene to dammarenediol II is a major factor responsible for the diversity of the flux through the ginsenoside biosynthetic pathway in various ginseng cultivars and that DDS might play a key role in the upregulation of ginsenoside biosynthesis.

CONCLUSION
We investigated the diversity of the ginsenoside biosynthetic pathway among adventitious roots from five ginseng cultivars using integrated transcriptomic and metabolomic analysis. We elucidated the variation in the ginsenoside biosynthetic pathway by comparing the expression levels of genes, the relative abundance of intermediate compounds, and the ginsenoside contents in adventitious roots of five genetically fixed inbred ginseng cultivars grown in a controlled environment. We also examined the intrinsic variation in ginsenoside contents and gene expression levels in the cultivars before and after MeJA treatment. Major genes in the triterpene biosynthesis pathway, including SQE and DDS, were found to play an important role in regulating ginsenoside biosynthesis, as well as in the variation of ginsenoside contents. The genes responsible for the metabolic diversity of these ginseng cultivars are good candidates to be utilized for the genetic improvement of this important medicinal plant.

AUTHOR CONTRIBUTIONS
T-JY designed the research and organized the manuscript. YSL and H-SP conducted adventitious root culture and transcriptome analysis. D-KL and SWK conducted metabolome analysis. YSL, H-SP, D-KL, MJ, N-HK, HJK, S-CL, and YJK conducted analysis of integrated bioinformatics data. YSL, H-SP, and T-JY wrote and revised the manuscript. All authors approved the final manuscript.

ACKNOWLEDGMENT
This work was carried out with the support of the "Cooperative Research Program for Agriculture Science & Technology Development (PJ01103001), " Rural Development Administration, Republic of Korea.