Transcriptional Responses and Gentiopicroside Biosynthesis in Methyl Jasmonate-Treated Gentiana macrophylla Seedlings

Gentiana macrophylla, a medicinal plant with significant pharmacological properties, contains the bioactive compound gentiopicroside. Methyl jasmonate (MeJA) is an effective elicitor for enhancing the production of such compounds. However, little is known about MeJA-mediated biosynthesis of gentiopicroside. We investigated this phenomenon as well as gene expression profiles to determine the molecular mechanisms for MeJA-mediated gentiopicroside biosynthesis and regulation in G. macrophylla. Our HPLC results showed that Gentiana macrophylla seedlings exposed to MeJA had significantly higher concentrations of gentiopicroside when compared with control plants. We used RNA sequencing to compare transcriptional profiles in seedlings treated for 5 d with either 0 μmol L-1 MeJA (C) or 250 μmol L-1 MeJA (M5) and detected differentially expressed genes (DEGs). In total, 77,482 unique sequences were obtained from approximately 34 million reads. Of these, 48,466 (57.46%) sequences were annotated based on BLASTs performed against public databases. We identified 5,206 DEGs between the C and M5 samples, including genes related to the α-lenolenic acid degradation pathway, JA signaling pathway, and gentiopicroside biosynthesis. Expression of numerous enzyme genes in the glycolysis pathway was significantly up-regulated. Many genes encoding transcription factors (e.g. ERF, bHLH, MYB, and WRKY) also responded to MeJA elicitation. Rapid acceleration of the glycolysis pathway that supplies precursors for IPP biosynthesis and up-regulates the expression of enzyme genes in that IPP pathway are probably most responsible for MeJA stimulation of gentiopicroside synthesis. Our qRT-PCR results showed that the expression profiles of 12 gentiopicroside biosynthesis genes were consistent with the RNA-Seq data. These results increase our understanding about how the gentiopicroside biosynthesis pathway in G. macrophylla responds to MeJA.

Gentiana macrophylla Pall (family Gentianaceae) is a perennial medicinal plant prescribed in China since ancient times to treat arthralgia, stroke, hemiplegia, pain, jaundice, infantile malnutrition, and osteoarthritis [8]. Its dried roots are officially listed in the Chinese Pharmacopoeia under the name Radix Gentianae Macrophyllae (Qin-jiao in Chinese) and are frequently used to dispel rheumatism and ease pain [9]. Gentiopicroside, an abundant and indicative ingredient in Qin-jiao, is the most important active component of total secoiridoid glycosides and has significant anti-inflammatory, analgesic, and antibacterial properties, as well as biological activity for treating osteoarthritis and strengthening gastric motility [10,11]. Although levels of gentiopicroside in Qin-jiao are influenced by soil elements and fertilization [12,13], the effect of MeJA in regulating those concentration has not been investigated.
Gentiopicroside is synthesized via the secoiridoid pathway. This pathway has been well studied and reviewed in Catharanthus roseus [14,15], and most steps have been identified. In the first step within higher plants, isopentenyl diphosphate (IPP), a precursor of terpenoids, is formed through either the plastidial 2-C-methyl-D-erythritol-4-phosphate (MEP) pathway or the cytosolic mevalonic acid (MVA) pathway. The allylic isomer of IPP, dimethylallyl diphosphate, reacts with one IPP in a head-to-tail fashion to form geranyl diphosphate, which is then catalyzed and converted into geraniol. The secoiridoid pathway starts with geraniol and proceeds through a series of reaction steps leading to the formation of secologanin [15]. Ultimately, secologanin is converted into gentiopicroside and other secoiridoids through several currently unknown steps [16].
In contrast to conventional methods, such as single gene cloning and DNA microarrays, that yield a limited amount of genetic information, RNA-seq is powerful tool for analyzing differential gene expression with high resolution at the whole-genome level [17,18]. In particular, transcriptome analysis can reveal relationships between plant gene expression and phenotype [19][20][21]. No previous regulatory mechanisms for gentiopicroside biosynthesis have been reported for G. macrophylla. However, de novo analysis using next-generation sequencing technologies can provide a robust platform for elucidating the mechanisms that might influence the accumulation of gentiopicroside in that species.
Because data are lacking for the means by which gentiopicroside production in G. macrophylla is modulated by MeJA, we monitored concentrations of that compound in MeJAtreated seedlings. Methyl jasmonate stimulates gentiopicroside biosynthesis. Therefore, RNAseq can be used to analyze differential gene expression over time in MeJA-treated plants versus the control. In general, applications of MeJA trigger profound transcriptional reprogramming in plant cells to manipulate the machinery that controls a wide range of metabolite biosynthesis via interplay of both positive and negative regulators [22]. Improving our understanding of the events between MeJA application and gentiopicroside accumulation will be useful for developing strategies to enhance production in G. macrophylla. Therefore, our objective here was to identify the relevant metabolic pathways for major MeJA-responsive genes and decipher the molecular mechanism by which MeJA stimulates yields.

Plant growth conditions and treatments
Seeds of Gentiana macrophylla were collected from the cultivation base of our lab, which is located in Taibai county, Shaanxi province, China. They were surface-sterilized with 2% sodium hypochlorite for 9 min, followed by five rinses with distilled water. After being kept in the dark for 2 d, they were germinated on culture dishes containing an MS solid medium (16-h photoperiod, 20˚C). One-month-old seedlings were transferred to either a fresh MS solid medium (control; 0 μmol L -1 MeJA) or an MS medium supplemented with 250 μmol L -1 MeJA. Whole seedling samples were collected after 1, 3, 5, 7, 9, and 11 d of treatment and dried to constant weight at 40˚C for the determination of gentiopicroside contents. The experiments were performed in three individual biological replications and every treatment contained more than 30 seedlings.

Measurements of gentiopicroside
Gentiopicroside concentrations were determined via High Performance Liquid Chromatography (HPLC). All sample solutions and stock solutions of gentiopicroside were prepared as we have described before [23]. Chromatographic separations were conducted with a C18 column (250 × 4.6 mm, 5 μm particle size; Agilent Technologies Inc., USA) on an Agilent 1260 Infinity LC system, using a solvent system comprising 70% ddH 2 O (A) and 30% methanol (B). The flow rate was adjusted to 0.8 mL min -1 and the detection wavelength was 245 nm. All separations were performed at 25˚C.

RNA isolation
Total RNA was isolated using TRIzol 1 Reagent (Invitrogen, USA) according to the manufacturer's protocol. Quality of the RNA was assessed on agarose gels and the concentration was determined with a NanoDrop ND2000 Spectrophotometer (NanoDrop Technologies Inc., USA).

cDNA library construction and sequencing
Two cDNA libraries-C (control) and M5 (MeJA treatment for 5 d)-were generated using mRNA-Seq Sample Prep Kits (Illumina, USA) according to the manufacturer's instructions. Magnetic beads containing poly-T molecules were used to isolate the poly(A) mRNA from 20 μg of total RNA. Following purification, the samples were fragmented into small pieces using divalent cations at 94˚C for 5 min, then converted into first-and second-strand cDNA with a SuperScript double-stranded cDNA synthesis kit (Invitrogen). The synthesized cDNA was subject to end repair and adenylation of the 3' ends and purified using a QIAquick PCR Purification Kit (QIAGEN, Germany). Afterward, Illumina paired-end adapters were ligated to the resulting cDNA fragments. Each cDNA library was constructed with an insert size of 200 bp. After quality was verified on an Agilent 2100 Bioanalyzer, deep-sequencing was performed with an Illumina HiSeq4000. In all, 150 bp paired-end reads were generated.

De novo assembly and gene annotation
Raw reads were filtered by the Illumina pipeline prior to assembly. We removed any reads that showed adapter contamination or for which more than 20% of the bases had quality values 10 or more than 5% were unknown nucleotides. The high quality clean reads were then randomly clipped into overlapping K-mers with default K = 25 for assembly with the Trinity, a short-reads assembling program [24]. The resulting sequences, termed unigenes, from each sample's assembly were further processed for sequence-splicing and removal of all redundancy to acquire non-redundant unigenes that were as long as possible. Finally, Blast X alignments were made (E-values <0.00001) between the unigenes and protein databases NR, SwissProt, Pfam, KEGG, KOG, and COG. Sequence directions were decided and functional annotations were assigned for the unigenes based on the best alignment results. Combining NR annotation with the Blast2GO program (v2. 5

Analysis of differentially expressed genes (DEGs)
Normalized expression levels for the unigenes were calculated per the FPKM method [28]. The significance of differential transcript abundance was valued according to the false discovery rate (FDR) [29]. Only those differentially expressed genes (DEGs) with FDRs 0.001 and absolute fold-changes !2 were reserved. For examining pathway enrichment, all DEGs were mapped to terms in the KEGG database to identify significantly over-represented metabolic pathways or signal transduction pathways. We primarily focused on differentially regulated pathways closely related to the biosynthesis of gentiopicroside.

Validation of gene expression by qRT-PCR
We selected 12 unigenes involved in gentiopicroside production for qRT-PCR experiments. Gene-specific primer pairs were designed by Primer Premier 5.0 software, and SAND1 served as the reference gene [30]. Total RNA was isolated from the C and M5 samples using TRIzol 1 Reagent (Invitrogen) according to the manufacturer's protocol. After treatment with DNase I (Tiangen, China), 1 μg of RNA was used in reverse-transcription with PrimeScript TM 1st Strand cDNA Synthesis Kits (TaKaRa, Japan). Quantitative reactions were performed on a LightCycler 1 96 real-time PCR detection system (Roche, Switzerland), using SYBR_ Premix Ex Taq™ (TaKaRa, Japan). Reaction conditions included an initial 95˚C for 10 min, then 40 cycles of 95˚C for 15 s, followed by 60˚C for 25 s. Relative expression levels for each unigene were compared between the two sample types and were calculated by the 2 -ΔΔCt method [31]. All data were expressed as means ± SD after normalization. Primer sequences used for qRT-PCR are listed in S8 Table. Availability of supporting data The data set supporting the results presented in this paper is available from the NCBI Sequence Read Archive repository (http://www.ncbi.nlm.nih.gov/sra/) under Accession Number SRP078971.

Effects of MeJA on gentiopicroside biosynthesis in Gentiana macrophylla seedlings
Induction of gentiopicroside biosynthesis in response to MeJA has not been previously described. Here, gentiopicroside production was stimulated in one-month-old G. macrophylla seedlings treated with 250 μmol L -1 MeJA (Fig 1). Our HPLC results showed that levels of gentiopicroside increased to 1.23, 1.76, 2.58, 2.16, and 2.07 mg g -1 , concentrations that were 19.3, 59.9, 123.3, 75.9, and 64.4% higher than those measured in corresponding control plants (C; 0 μmol L -1 MeJA) on Days 1, 3, 5, 7, and 9, respectively, of this experiment. This is the first report to demonstrate that MeJA treatment influences gentiopicroside biosynthesis in this species. We utilized G. macrophylla seedlings exposed to MeJA for 5 d for RNA-seq analysis because this span of time is associated with the greatest contrast in gentiopicroside concentrations between control (untreated) plants and those to which MeJA has been applied.

Sequencing and de novo assembly
Two cDNA libraries from G. macrophylla seedlings-C (control) and M5 (MeJA treatment for 5 d)-were sequenced by Illumina deep sequencing to obtain approximately 18.16 and 15.74 million high-quality clean reads for C and M5, respectively. Each read averaged 298 bp long. The Q30 value (percentage of sequences with a sequencing error rate <0.1%) was 94.93% for C and 94.28% for M5. Read lengths were 150 bp×2. Assembling all trimmed reads produced 14,952,699 contigs from the two libraries, which were then joined into unigenes based on the paired-end information. This generated 54,841 (C) and 62,179 (M5) unigenes (total of 77,482 unigenes, including overlaps between C and M5). These had an average length of 764 bp and an N50 of 1,169 bp (i.e., 50% of the assembled bases were incorporated into unigenes at least 1,169 nt long) ( Table 1). All unigenes were longer than 200 bp, 46.55% were more than 500 bp long, and 14.56% (12,559) were longer than 1,000 bp.

Functional annotation and classification
We performed function annotation of the generated unigenes by using BLASTX to search reference sequences against results from the NCBI non-redundant protein databases (NR),  Table).

Analysis of KEGG pathways and differentially expressed genes (DEGs)
We conducted a KEGG pathway-based analysis to obtain a better understanding of the biological functions of these unigenes. Among the 18,720 annotated transcripts mapped to 128 KEGG pathways, 20 of the most-represented pathways are shown in S2 Table. In all, 2,443 unigenes were assigned to 25 secondary-metabolite pathways ( Table 2). The Fragments Per Kb per Million reads (FPKM) method was used to calculate unigene expression levels and, ultimately, identify DEGs. When M5 and C treatments were compared, 3,805 unigenes were up-regulated and 1401 were down-regulated in response to MeJA applications. All DEGs are shown in S3 Table while the 20 unigenes most up-or down-regulated between M5 and C samples are presented in S4 Table. Our KEGG analysis enabled us to identify 38 significantly enriched metabolic pathways or signal transduction pathways for those DEGs (Fig 4). To examine the molecular basis for jasmonate (JA) stimulation of gentiopicroside production in G. macrophylla, we focused on glycolysis, the IPP biosynthesis pathway, α-linolenic acid metabolism, plant hormone signal transduction, and related transcription factors (TFs). DEGs involved in JA biosynthesis and the JA signaling pathway As an elicitor of the production of bioactive metabolites, JA triggers a transcriptional reprogramming of plant metabolism, resulting in a concerted upregulation of expression for genes that encode enzymes involved in specific, specialized metabolic pathways [32]. We investigated changes in the expression of genes closely associated with α-linolenic acid metabolism, which finally leads to JA biosynthesis, and were interested to find 22 DEGs (Table 3; Fig 5). Of those, 19 DEGs encoding putative lipoxygenase (LOX), allene oxide synthase (AOS), allene oxide cyclase (AOC), 12-oxophytodienoic acid reductase (OPR), OPC-8:0 CoA ligase 1 (OPCL1), acyl-CoA oxidase (ACOX), acetyl-CoA acyltransferase 1 (ACAA1), and jasmonate O-methyltransferase (JMT) showed increased transcript abundance (Table 3; Fig 5).

DEGs involved in gentiopicroside biosynthesis
Our RNA-seq data indicated that 243 unigenes, including 59 DEGs, are members of the gentiopicroside biosynthesis pathway. Among them, 114 unigenes, including 23 DEGs, are involved in the IPP pathway while 129 unigenes, including 36 DEGs, function in secoiridoid biosynthesis (Table 4). Although we did not find here that the latter pathway was significantly enriched, we did note that the former pathway is (Fig 6). Expression of 21 DEGs related to the IPP pathway was up-regulated in plants treated with MeJA. Among the 36 DEGs involved in iridoid biosynthesis, 18 were up-regulated and 18 were down-regulated. These results demonstrated that applying MeJA could significantly increase gentiopicroside biosynthesis by up- regulating the expression of genes related to the IPP pathway, but not to the secoiridoid biosynthesis pathway.

DEGs involved in glycolysis
Glucose is broken down by glycolysis to produce acetyl coenzyme A (acetyl-CoA) as the direct precursor in the MVP pathway, as well as glyceraldehyde-3-phosphate and pyruvate, which are precursors in the MEP pathway [33]. We found that the glycolysis pathway was significantly enriched in M5 samples when compared with C samples (Fig 7) and 52 of 57 related DEGs were up-regulated (S5 Table). For example, in that pathway, many putative genes encoding hexokinase, fructose-bisphosphate aldolase, triosephosphate isomerase, glyceraldehyde 3-phosphate dehydrogenase, phosphoglycerate kinase, gpmI, enolase, pyruvate kinase, phosphoenolpyruvate carboxykinase, pyruvate dehydrogenase E1 component alpha subunit, pyruvate decarboxylase, dihydrolipoamide acetyltransferase, L-lactate dehydrogenase, dihydrolipoamide dehydrogenase, and aldehyde dehydrogenase were significantly up-regulated by MeJA treatment (Fig 7, S5 Table).

DEGs associated with hormone signaling components
We observed that treatment with MeJA triggered the enrichment of signal transduction pathways for GA, ET, SA, and ABA. In addition to MYC2 and JAZs in the JA signaling pathway, transcripts of signaling components such as the gibberellin receptor GID1, DELLA, and phytochrome-interacting factor 4 for GA; ethylene receptor ETR, ethylene-insensitive protein 3, and ethylene-responsive transcription factor 1 (ERF1) for ET; pathogenesis-related protein 1 (PR1) for SA; and protein phosphatase 2C and serine/threonine-protein kinase for ABA were greatly accumulated in response to MeJA (S6 Table). We were interested to find that expression of the PR1 mRNA was increased by 576-fold.

DEGs associated with TFs
We identified 164 DEGs that responded to MeJA elicitation, including 131 that were up-regulated and 33 that were down-regulated (S7 Table). They were largely represented by TF families that influence secondary metabolism and stress responses, e.g., the ERF superfamily (32 members), bHLH superfamily (27 members), WRKY superfamily (19 members), and myeloblastosis (MYB) superfamily (15 members). Transcript levels of the five most up-regulated genes increased by 105.3-to 212.1-fold in M5 samples, with three of them putatively encoding members of the ERF superfamily. In contrast, transcript levels of the five most strongly down-

Validation of gene expression by qRT-PCR
We used qRT-PCR analysis to validate the important DEGs obtained from our assembled transcriptome as well as from the expression profiles revealed by RNA-Seq data. This examination entailed 12 unigenes involved in the biosynthesis of IPP (DXS, DXR, HDS, HDR, HDS, HMGR, and MVD) and secoiridoids (GES, G10H, and 8/10HGO) (Fig 8). Our results suggested that the assembled transcripts were reliable and that the designed primer pairs were suitable for subsequent expression experiments. Based on the 2 -ΔΔCt method [31], relative expression levels of the selected unigenes were calculated and compared between M5 and C samples. The expression patterns detected by qRT-PCR were consistent with those from the  RNA-Seq data. Overall, the qRT-PCR analysis confirmed that the unigenes obtained from the assembled transcriptome were trustworthy and the profiles were credible.

Discussion
Exogenous MeJA is believed to be a primary regulator of pathways for JA biosynthesis and signaling in plants. This compound has been studied extensively in Solanum lycopersicum, Arabidopsis thaliana, and Taxus chinensis [6,34,35]. Jasmonates are plant-specific signaling molecules that activate several defense mechanisms, inducing a massive reprogramming of gene expression [36]. The basic helix-loop-helix (bHLH) TF MYC2 is a central regulator in JA signaling cascades, including those leading to the biosynthesis of several classes of specialized metabolites [37]. In the absence of JA, MYC2 action is repressed by jasmonate ZIMdomain-containing (JAZ) proteins by forming repressor complexes with a group of other Table 4  . This JA elicitation leads to JAZ degradation by the SCFCOI1-ubiquitin-proteasome pathway, thereby releasing MYC2 from the repressor complex. In the present transcriptome pool, there are 7 putative unigenes encoding MYC2 and 18 putative genes encoding JAZs based on BLASTs performed against public databases (S1 Table). We found that the expression of all DEGs encoding putative MYC2 and JAZ were up-regulated in MeJA-treated seedlings ( Table 3). The JA signaling mechanism also oscillates through a negative feedback loop involving MYC2 and JAZ proteins, in which JAZ blocks MYC2 activity at the protein level while MYC2 transcriptionally induces JAZ expression [34]. Our results clearly confirmed that exogenous application of MeJA can regulate the pathways for JA biosynthesis and signaling in G. macrophylla. The derivation of gentiopicroside from secologanin entails 10 enzymatic conversions, starting from IPP (Fig 6). Catharanthus roseus has been used to investigate several genes that encode key enzymes in those pathways [15,40]. Briefly, IPP is produced via the plastidial MEP or cytosolic MVA pathway (Fig 6). However, the biosynthetic route from secologanin to gentiopicroside has been entirely unknown to date. We previously identified 114 putative unigenes involved in secoiridoid biosynthesis in our transcriptome library [41]. In the present database, 40 putative unigenes encoding iridoid synthase (IS), iridoid oxidase (IO), 7-deoxyloganetic acid glucosyltransferase (7-DLGT), and 7-deoxyloganic acid hydroxylase (7-DLH) were identified for the first time in G. macrophylla.

. Putative genes that encode enzymes involved in gentiopicroside biosynthesis and DEGs for which expression is altered in response to
The glycolysis pathway (Fig 7) not only plays a crucial role in energy generation but also provides carbon building blocks for the biosynthesis of gentiopicroside and other organic constituents of secoiridoid [42,43]. In particular, glucose is broken down by glycolysis to produce acetyl coenzyme A (acetyl-CoA) as the direct precursor in the MVP pathway, as well as glyceraldehyde-3-phosphate and pyruvate, which are precursors in the MEP pathway [33]. In the present study we found that the glycolysis pathway was significantly enriched in M5 samples when compared with C samples (Fig 7) and up-regulation of numerous enzyme genes led to an elevated rate of flux and replenished the precursors consumed in IPP biosynthesis. Research on the rubber tree (Hevea brasiliensis) has suggested that such rapid acceleration of the glycolysis pathway in supplying precursors for the production of IPP and natural rubber, rather than rubber biosynthesis per se, is responsible for ethylene (ET) stimulation of latex yields [44].
Here, we noted that the high expression of enzyme genes in both the glycolysis pathway and the IPP biosynthesis pathway contributed to greatly increased gentiopicroside concentrations in response to MeJA.
Hormone responses are generally the result of interactions and crosstalk among multiple pathways [45]. For example, JA signaling may function by interacting with other major plant hormones, such as ET, gibberellin (GA), alicylic acid (SA), abscisic acid (ABA), brassinosteroid, auxin, and cytokinin. Signaling crosstalk between SA and JA results in complementary action in mediating endophyte-induced accumulations of secondary metabolites [46,47]. Moreover, under stress conditions, ABA interacts with the SA/JA pathways [48]. Although some plant hormones, e.g., ABA, ET, and SA, induce the production of bioactive compounds [44,49], it is unclear how JA interacts and coordinates with other hormones with regard to this stimulation. Our results showed that treatment with MeJA triggered the enrichment of signal transduction pathways for GA, ET, SA, and ABA. One explanation for this crosstalk is that IPP is also the precursor of GA and ABA biosynthesis and that MeJA also elicits a rise in IPP production.
Transcription factors play important roles in controlling many biological processes in a cell or organism by modulating gene expression [50]. Many TFs help regulate the biosynthesis and accumulation of secondary metabolites [51], and also have crucial roles in crosstalk between hormone signalling pathways [52]. For example, CrWRKY1, a member of the WRKY family from Catharanthus roseus, can positively regulate TIA biosynthesis [53]. The phenylpropanoid pathway in different plant organs and tissues of higher plants is under the control of specific R2R3-MYB TFs and bHLH families [54]. We identified many putative transcription factors differentially expressed between C samples and M5 samples (S7 Table), including the ERF superfamily, bHLH superfamily, WRKY superfamily, and myeloblastosis (MYB) superfamily. Although it is still unknown how these genes function in parallel with MeJA applications to enhance the production of gentiopicroside, our data provide a valuable resource for discovering candidate genes related to the complex regulatory networks involved in that response.

Conclusion
Although MeJA is a ubiquitous and conserved elicitor of plant secondary metabolism, the degree to which metabolic pathways are stimulated is species-specific. We verified here that MeJA applications effectively enhance the production of gentiopicroside in Gentiana macrophylla. Our RNA-Seq analysis of MeJA-related transcriptional changes indicated that 5206 genes are differentially expressed. Transcriptome analysis revealed increased expression of genes in the α-lenolenic acid degradation pathway that produces abundant JA and quickly activates the holistic JA pathway in those seedlings. Rapid acceleration of the glycolysis pathway that supplies precursors for IPP biosynthesis and up-regulates the expression of enzyme genes in that IPP pathway are probably most responsible for MeJA stimulation of gentiopicroside synthesis. Furthermore, many genes encoding various TFs, e.g., ERF, bHLH, MYB, and WRKY, also respond in plants exposed to MeJA.
The results from this research improve our understanding of how MeJA applications alter the production of gentiopicroside in G. macrophylla seedlings. Our data will provide a massive genetic resource for further investigation of gentiopicroside biosynthesis and will lay the foundation for genetic engineering to boost yields of this compound in such a valuable medical plant.