Hybrid sequencing of the Gynostemma pentaphyllum transcriptome provides new insights into gypenoside biosynthesis

Background Gypenosides are a group of triterpene saponins from Gynostemma pentaphyllum that are the same as or very similar to ginsenosides from the Panax species. Several enzymes involved in ginsenoside biosynthesis have been characterized, which provide important clues for elucidating the gypenoside biosynthetic pathway. We suppose that gypenosides and ginsenosides may have a similar biosynthetic mechanism and that the corresponding enzymes in the two pathways may have considerable similarity in their sequences. To further understand gypenoside biosynthesis, we sequenced the G. pentaphyllum transcriptome with a hybrid sequencing-based strategy and then determined the candidate genes involved in this pathway using phylogenetic tree construction and gene expression analysis. Results Following the PacBio standard analysis pipeline, 66,046 polished consensus sequences were obtained, while Illumina data were assembled into 140,601 unigenes with Trinity software. Then, these output sequences from the two analytical routes were merged. After removing redundant data with CD-HIT software, a total of 140,157 final unigenes were obtained. After functional annotation, five 2,3-oxidosqualene cyclase genes, 145 cytochrome P450 genes and 254 UDP-glycosyltransferase genes were selected for the screening of genes involved in gypenoside biosynthesis. Using phylogenetic analysis, several genes were divided into the same subfamilies or closely related evolutionary branches with characterized enzymes involved in ginsenoside biosynthesis. Using real-time PCR technology, their expression patterns were investigated in different tissues and at different times after methyl jasmonate induction. Since the genes in the same biosynthetic pathway are generally coexpressed, we speculated that GpOSC1, GpCYP89, and GpUGT35 were the leading candidates for gypenoside biosynthesis. In addition, six GpWRKYs and one GpbHLH might play a possible role in regulating gypenoside biosynthesis. Conclusions We developed a hybrid sequencing strategy to obtain longer length transcriptomes with increased accuracy, which will greatly contribute to downstream gene screening and characterization, thus improving our ability to elucidate secondary metabolite biosynthetic pathways. With this strategy, we found several candidate genes that may be involved in gypenoside biosynthesis, which laid an important foundation for the elucidation of this biosynthetic pathway, thus greatly contributing to further research in metabolic regulation, synthetic biology and molecular breeding in this species. Electronic supplementary material The online version of this article (10.1186/s12864-019-6000-y) contains supplementary material, which is available to authorized users.


Background
Gynostemma pentaphyllum (Thunb.) Makino, belonging to the Cucurbitaceae family, is a kind of perennial creeping herbaceous plant widely distributed in China, India, Japan, Korea and Southeast Asia [1,2]. G. pentaphyllum, also known as "Southern Ginseng", is a traditional Chinese medicinal herb and is reported to have an adaptogenic nature that enhances the "yin" and "yang" properties of the human body [3]. The main active components of G. pentaphyllum are gypenosides (triterpene saponins). To date, more than 200 different gypenosides have been isolated from G. pentaphyllum [4][5][6][7]. Modern pharmacological studies have shown that gypenosides contribute to antitumor, hypoglycemic, hypolipidemic, cardiovascular and cerebrovascular protection and immunoprotection [8][9][10]. Previous studies have shown that 25% of all gypenosides are similar to ginsenosides, especially Gyp III, IV, VIII and XII, which are exactly the same as Gin-Rb1, -Rb3, -Rd and -F2 [11]. The content of triterpene saponins in G. pentaphyllum is almost five times higher than that of Panax ginseng. In addition, G. pentaphyllum is much easier to culture and grows faster than P. ginseng. Therefore, G. pentaphyllum is a good alternative resource for gypenoside production and has attracted great interest worldwide [3,12].
RNA sequencing (RNA-Seq) based next-generation sequencing (NGS) technology is a powerful and economical way to discover novel genes on a large scale. This technology has been widely used to uncover genes involved in the biosynthesis of secondary metabolites. Although NGS has advantages regarding sequencing depth and cost, the short read length may restrict correct assembly and annotation [32,33]. More recently, isoform sequencing (Iso-Seq) based on Pacific Biosciences (PacBio) SMRT technology has been developed. Compared to RNA-Seq, Iso-Seq has longer reads and an assembly-free analysis pipeline, which provide more fulllength transcripts and direct evidence of the structural variation of isoforms [34,35]. However, the higher cost and relatively low sequencing quality of Iso-Seq limit its wider application in gene mining.
Here, we provide a hybrid sequencing strategy that can produce high-quality long transcripts for gene mining and identification. Briefly, we used SMRT sequencing technology and NGS to generate comprehensive insight into full-length sequences from root, stem and leaf tissues. Moreover, we estimated the candidate enzymes involved in the biosynthesis pathway, OSCs, CYP450, and UGT through a methyl jasmonate (MeJA)-treated experiment and tissue-specific expression analyses. Finally, one OSC, one CYP450, one UGT, six WRKYs and one bHLH were selected as the candidates most likely to be involved in gypenoside biosynthesis. Our results will provide a valuable resource for investigating novel genes in the biosynthesis of gypenosides.

Hybrid sequencing, assembly and annotation
To obtain as many high-quality unigenes as possible, a hybrid sequencing strategy that combined SMRT and NGS technology was used to analyze the transcriptome of G. pentaphyllum (Fig. 2). A total of 13.17 Gb raw data were generated from eight SMRT cells by Iso-Seq. Following the PacBio standard analysis pipeline, 268,927 reads of insert were generated and then classified into 99,739 full-length non-chimeric reads and 142,079 nonfull-length reads. The full-length non-chimeric sequences were clustered by ICE and then polished by Quiver to obtain 66,046 polished consensus sequences. In addition, 85.82 Gb of Illumina raw data were obtained by sequencing the transcriptomes from roots, stems, leaves and MeJA-treated leaves (Additional file 1: Table S1). To obtain more transcripts for gene discovery, all raw data from different samples were combined. After removal of adapter sequences and low-quality reads, 576,532,682 clean reads were generated. All clean reads were then assembled with Trinity software to generate 140,601 Trinity unigenes. The polished consensus sequences and Trinity unigenes were merged. Finally, after removing redundant data with CD-HIT software, 140,157 final unigenes were obtained, ranging from 200 bp to 13,760 bp in size, with an average length of 750 bp. The transcript length distribution generated using these two platforms is shown in Fig. 3a. Trinity unigenes were mostly < 500 bp, and the average length was 650 bp. PacBio polished consensus sequences and PacBio unigenes were mostly distributed in the 800-2500 bp range. Their average lengths were also longer than Trinity unigenes, at 1564 bp and 2075 bp, respectively. These results demonstrate that PacBio sequencing reads are an invaluable resource to generate longer full-length transcripts without assembly, an element of critical importance for genomic studies on species without a reference genome assembly. All unigenes were functionally annotated by searching against the SwissProt, KEGG, Pfam and GO databases. The results indicated that 69,101 (49.3%) unigenes showed annotation hits in at least one of four databases, and 16,007 unigenes showed annotation hits in all databases (Fig. 3b). In total, 126,493 unigenes were assigned GO terms, which were classified into three major categories (molecular function, cellular component and biological process; Additional file 2: Figure S1). In the biological process category, the major subgroups were "metabolic process" (GO: 0008152) and "cellular process" (GO: 0009987). In the cellular component category, unigenes categorized in the "cell" (9416, 7.44% of the total), "membrane" (9354, 7.39% of the total) and "cell part" (9258, 7.31%) were highly represented. In the molecular function category, the major categories were "catalytic activity" (GO:0003824) and "binding" (GO: 0005488). Furthermore, based on KEGG annotation, 41, 549 unigenes were categorized into four major categories with different functions (Cellular Processes, Environmental Information Processing, Genetic Information Processing and Metabolism). Unigenes related to energy metabolism represented the largest subcategory, followed by genes related to amino acid metabolism and carbohydrate metabolism (Fig. 3c). More importantly, 614 unigenes were annotated as related to terpenoid metabolism and polyketides, accounting for 2.7% of the metabolism category, and 190 unigenes were related to terpenoid backbone biosynthesis. These findings are helpful for mining the candidate genes involved in gypenoside biosynthesis.

Candidate OSCs involved in gypenoside biosynthesis
OSCs catalyze the first committed step in triterpene biosynthesis, namely, the cyclization of the universal triterpene precursor 2,3-oxidosqualene and therefore define sterol and triterpene skeletal diversity. More than 100 diverse triterpene skeletons are currently known in plants.
In this study, we obtained 5 full-length candidate OSCs after data mining and manual curation (GpOSC1-5). Alignment analysis showed that all candidates contained six repeats of QW (QXXXGXW) motifs and one DCTAE motif (Additional file 3: Figure S2). The former can enhance the structure of the protein and stabilize the carbocation intermediates during cyclization, while the latter is presumed to be responsible for the initiation of the cyclization reaction [36,37]. Additionally, GpOSC2 and GpOSC4 have one MWCHCR motif, within which the histidine residue has been presumed to be important for stabilization of the protosteryl cation intermediate [38]. The protosteryl cation intermediate is catalyzed by the enzymes cucurbitadienol synthase (CbQ), CAS and LS.
To further predict the function of these OSCs, we analyzed the phylogenetic relationship between the G. pentaphyllum OSCs and characterized OSCs in other plants (Fig. 4). In the phylogenetic tree, LS, CAS and CbQ group together, while DS, BAS and lupeol synthase (LUS) form another group, which is in concordance with previous studies. The two groups of OSCs have different In general, sterol skeletons are synthesized through the CBC cyclization pathway, and triterpene skeletons are synthesized through the CCC cyclization pathway. As shown in Fig. 4, GpOSC2 belonged to the CAS clade, which showed 88% identity with the CAS from Cucurbita pepo, demonstrating that GpOSC2 might have a similar catalytic activity with CAS. GpOSC3 belonged to the LS clade and that it had 77% identity with the LS from Luffa aegyptiaca. GpOSC4 is predicted as a CbQ, which can catalyze 2,3-oxidosqualene to cyclize to cucurbitacins. Cucurbitacins are widely distributed in plants in the Cucurbitaceae family. GpOSC1 and GpOSC5 phylogenetically belonged to the BAS clade and had 70 and 64% identity with the identified BAS of P. quinquefolius, respectively. Simultaneously, GpOSC1 had 57% identity with DS from P. ginseng. However, no β-amyrin-type saponins were isolated from G. pentaphyllum, and the two enzymes need further functional characterization.
In general, the genes encoding enzymes in the same pathway are coexpressed. Therefore, with the four known genes involved in triterpene biosynthesis as references, we analyzed the expression levels of the 5 GpOSCs in different tissues using real-time PCR (Fig. 5a). 3-hydroxy-3-methyl glutaryl coenzyme A reductase (HMGR), FPS, SE and SS were found to be tightly coexpressed; they were all expressed highest in leaves, followed by stems and then roots. The results agree with the gypenoside contents in these tissues. GpOSC5 was expressed highest in stems; other GpOSCs had similar expression patterns to those of the four upstream genes, although GpOSC3 had much lower expression levels than GpOSC1, GpOSC2 and GpOSC4. In particular, the expression of GpOSC1 in leaves was 551.84-fold greater than that in roots and 33.38-fold greater than that in stems.
In plants, MeJA is an important signaling molecule in response to abiotic or biotic stresses, which can induce the biosynthesis of many secondary metabolites for plant defense [39]. MeJA can induce the expression of genes involved in the ginsenoside synthetic pathway, such as FPS, SS, SE and DS, thereby increasing the content of ginsenosides in P. ginseng [40,41]. This study showed that HMGR, FPS, SE and SS in the gypenoside biosynthetic pathway can also be induced by MeJA in G. pentaphyllum (Fig. 5b). More interestingly, the expression of all the genes was depressed at 6 h after MeJA treatment and then increased gradually. Among them, SS, SE and HMGR increased between 6 h and 18 h, reaching a maximum expression of 3.54-fold, 1.64-fold and 4.98fold greater than that at 0 h, respectively. FPS sharply increased between 6 and 12 h, peaking at 12 h with a 4.20fold increase. Among the four GpOSCs, only GpOSC2 was not induced by MeJA. GpOSC1 had similar expression patterns with the four known genes in the gypenoside biosynthetic pathway. In particular, GpOSC1 and HMGR showed the highest similarity in expression pattern with the Spearman correlation coefficient of 1, suggesting their consistency in response to MeJA induction (Fig. 5c). According to coexpression and phylogenetic analyses, we speculated that GpOSC1 most likely functions as a DS.

Candidate CYP450s involved in gypenoside biosynthesis
CYP450s that catalyze structural modification are crucial for the diversification and functionalization of triterpenes. Based on the annotation results, 145 CYP450s with lengths from 300 to 581 amino acids were identified in G. pentaphyllum (GpCYP1-145). They were classified by alignment with the CYP450 database using standard sequence similarity cutoffs, specifically 40, 55 and 97% for family, subfamily and allelic variants, respectively [42]. Thus, the 145 GpCYP450s were classified into 40 families and 39 subfamilies (Additional file 4: Figure S3). However, seven CYP450s (GpCYP3, 50, 51, 83, 96, 97, 114) did not belong to any old family. We considered that they belonged to novel P450 families. In addition, 68 CYP450s were found to belong to novel P450 subfamilies.
To date, all characterized CYP450s involved in ginsenoside biosynthesis belonged to the CYP716A subfamily (CYP716 family, 85 Clan). We found that two GpCYP450s (GpCYP12 and GpCYP89) from G. pentaphyllum belong to the CYP716 family. Therefore, phylogenetic analysis was carried out with the 85 Clan GpCYP450s and three characterized CYP450s involved in ginsenoside biosynthesis from P. ginseng (Fig. 6). Considering that almost all gypenosides are protopanaxadiol-type saponins, GpCYP450s with high identity with protopanaxadiol synthase (CYP716A47) are most likely involved in gypenoside biosynthesis. GpCYP12 and GpCYP89 had 49 and 48% identity to CYP716A47, respectively. However, they had higher identity (approximately 57%) with β-amyrin-28 oxidase (CYP716A52v2) of P. ginseng, although we have not found any β-amyrin-type triterpene saponins from G. pentaphyllum until now.
The expression patterns of these two GpCYP450 genes in different tissues and MeJA-treated leaves were then further analyzed with real-time PCR. As shown in Fig. 5a, GpCYP89 had a similar tissue-specific expression pattern to that of the upstream genes. The expression of GpCYP89 in leaves was more than 11.53-fold higher than that of stems, and 265.41-fold higher than that of roots. Unsurprisingly, GpCYP89 can be induced by MeJA because most CYP450s are involved in secondary metabolism in plants (Fig. 5c). GpCYP89 had similar expression patterns to those of SS, SE and HMGR, which were downregulated before 6 h and then continuously increased and reached the highest expression at 18 h; this highest expression was 1.82-fold greater than that at 0 h, with Spearman correlation coefficients of 0.9 with HMGR. According to phylogenetic and expression analyses, GpCYP89 from the CYP716 family may have functioned as protopanaxadiol synthase in G. pentaphyllum.

Candidate UGTs involved in gypenoside biosynthesis
Glycosylation is the last step in gypenoside biosynthesis. In this study, we obtained a total of 254 unigenes annotated as UGTs. We selected 68 UGTs with peptide lengths greater than 352 amino acids (GpUGT1-68) and divided them into 20 UGT families based on the criteria of amino acid sequence identity cutoffs, briefly, 40% for family, and 60% for subfamily [43]. Nine GpUGTs were noted as novel UGT families (Fig. 7).
To date, five UGTs from P. ginseng have been shown to be involved in ginsenoside biosynthesis, which belong to the UGT71, UGT74 and UGT94 families. Phylogenetic analysis showed that a total of 13 GpUGTs clustered into the above three families (Fig. 7). Therefore, they are candidate UGTs for gypenoside biosynthesis. Among them, GpUGT35 was clustered into the UGT94 family with the highest identity (50%) to PgUGT94Q2. GpUGT24 was clustered into the same branch as PgUGT74AE2 with 42% identity. GpUGT8 showed the highest identity (46%) with UGTPg1 (UGT71 family) from P. ginseng.
According to phylogenetic analysis, 13 GpUGTs were chosen for analysis of their expression patterns by realtime PCR. In the tissue-specific pattern assay, five GpUGTs (GpUGT2, GpUGT3, GpUGT35, GpUGT46 and GpUGT62) showed strong similarity to FPS, SS, SE and HMGR (Fig. 5a). In particular, the expression of GpUGT35 in leaves was 28.00-fold higher than that in stems and 453.00-fold higher than that in roots. In the MeJA treatment experiment, the expression of all five GpUGTs can be regulated by MeJA (Fig. 5c). However, only GpUGT35 had a similar expression pattern to SS, SE and HMGR, which was continuously increased by MeJA between 6 h-18 h and then reached the highest at 18 h, which was 2.17-fold compared to that at 0 h. The Spearman correlation coefficient was 1 with HMGR. According to the phylogenetic and expression analysis, GpUGT35 from the UGT94 family was regarded as a candidate enzyme responsible for gypenoside biosynthesis and will be the subject of further study.

Candidate TFs possibly regulating gypenoside biosynthesis
In this study, we selected 1362 TFs distributed in 64 families, including ARF, B3, MYB, bHLH, bZIP, and WRKY. As shown in Fig. 8, the largest family is AP2/ ERF-ERF, accounting for 9.9% of the total, followed by bHLH and C2H2, accounting for 7.5 and 6.6%, respectively. To date, only two TFs have been found to be involved in the regulation of the triterpenoid pathway, PqWRKY1 from P. quinquefolius and PnbHLH from P. notoginseng [30,31].
WRKY proteins are one of the largest families of TFs in plants. The WRKY domain can specifically bind to the target gene promoter W-box sequence to regulate the expression of target genes [44]. PqWRKY1 is reported as a positive regulator related to gypenoside biosynthesis in American ginseng [30]. A total of 58 putative WRKY proteins were identified by transcriptome-wide identification from G. pentaphyllum, among which 38 WRKYs (GpWRKY1-38) had complete WRKY domains (Additional file 5: Figure S4). According to Eulgem's method, 38 GpWRKYs were categorized into Groups I, II and III, and Group II can be further classified into subgroups IIa, IIb, IIc, IId and IIe. (Additional file 6: Table S2) [45]. PqWRKY1 belonged to IIc, and 13 GpWRKYs from G. pentaphyllum were clustered into this subgroup. The expression of 13 GpWRKYs was analyzed with real-time PCR. Among them, six GpWRKYs (GpWRKY7, GpWRKY11, GpWRKY14, GpWRKY17, GpWRKY26 and GpWRKY33) showed the highest expression in leaves, compared with roots and stems (Fig. 5a). Furthermore, the expression of all of these 6 GpWRKYs was increased by MeJA treatment and reached the highest level (> 4.56-fold) at 24 h compared with that at 0 h (Fig. 5d).
bHLH proteins are one of the largest families of TFs in both plants and animals, and play important roles in plant metabolism [46]. PnbHLH1 from P. notoginseng is regarded as another positive regulator involved in ginsenoside biosynthesis [31]. A total of 89 TFs with complete bHLH domains (Additional file 7: Figure S5) were grouped into 21 subgroups by phylogenetic tree construction [47,48], with the XII subgroup having the largest number of members (13 sequences) and the VIIIa, II, and IIIf subgroups having the fewest (1 sequence)  Table S3). PnbHLH1 belonged to IVa and four GpbHLHs (GpbHLH32, GpbHLH45, GpbHLH48, GpbHLH68) clustered with PnbHLH1 together. The expression of four GpbHLHs was then analyzed with real-time PCR. As shown in Fig. 5a, only GpbHLH68 was highly expressed in leaves. Moreover, GpbHLH68 expression was increased by MeJA treatment and reached the highest level (~2-fold) at 18 h compared with that at 0 h, suggesting its possible positive role in regulating gypenoside biosynthesis (Fig. 5d). Based on the above phylogenetic and expression analysis, six GpWRKYs and one GpbHLH were regarded as possible positive regulators of gypenoside biosynthesis and will be the subject of further study.

Discussion
Transcriptome sequencing is a powerful and economical tool for obtaining genetic information on a large scale in organisms without available genome sequences. Undoubtedly, low-cost NGS has provided deep insights into the characterization and quantification of the transcriptome; however, the relatively short reads generated by NGS have prohibited the accuracy of transcript reconstruction after assembly with software. In most cases, this assembly step may suffer from a misassembly of the sequences, lowering the accuracy of gene annotation and expression quantification in the following steps. The third generation sequencing technology, SMRT sequencing developed by PacBio, offers an alternative way to overcome this limitation. SMRT sequencing with a read length up to 20 kb renders the PacBio platform very effective in the sequencing of full-length cDNAs that exhibit long-transcript isoforms, avoiding the mistakes introduced by assembly steps [49]. However, the sequencing accuracy of SMRT technology is much lower than that of NGS. An alternative way to overcome these limitations is to integrate NGS short reads and PacBio long reads via hybrid sequencing. In hybrid sequencing, NGS short reads are mapped to and correct SMRT reads and could improve the accuracy of gene isoform identification and abundance estimation.
Natural products are important sources for drug discovery. Modern pharmacological studies have shown that ginsenosides have great potential in the prevention and treatment of multiple cancers [3]. Ginsenosides are generally produced by plants from the Panax genus, such as P. ginseng, P. quinquefolius and P. notoginseng. However, these plants grow slowly and cannot be continuously cropped in the same fields [12]. Thus, the ginsenosides isolated from the Panax plants are not able to meet the need for drug production and development. G. pentaphyllum can produce some triterpenoids named gypenosides, which have structures that are the same as or similar to those of ginsenosides [11]. Compared to Panax plants, G. pentaphyllum grows much faster, can be harvested four times a year and has a higher triterpenoid content. Therefore, this plant is a promising alternative resource for ginsenoside production. In addition, because of its short life cycle, easy tissue culture and genetic transformation, it is also an emerging medicinal model plant for triterpenoid biosynthesis.
Most ginsenosides and gypenosides have dammaranetype aglycones, which are not widely distributed triterpenoids, mainly produced by the Panax and Gynostemma genera [4][5][6][7]. Interestingly, even though the two genera belong to different evolutionary branches that are very far from each other, they can produce the same triterpenoids. DS is the signature enzyme of the dammarane biosynthetic pathway [16]. We found five OSCs in the G. pentaphyllum transcriptomes. Phylogenetic analysis showed that no GpOSCs clustered with DSs from the Panax genus, while two GpOSCs belonged to the BAS group, the group nearest to the DS group in the phylogenetic tree, among which GpOSC1 showed the highest similarity with the characterized DS (Fig. 4). Since the genes in the same biosynthetic pathway are generally coexpressed, we compared the expression patterns of all GpOSCs with four upstream genes encoding HMGR, FPS, SE and SS and found that the expression patterns of GpOSC1 were the most similar to those of the upstream genes (Fig. 5). Combined phylogenetic and expression analyses were used to screen two tailoring enzymes, GpCYP89 and GpUGT35, which may be involved in gypenoside biosynthesis. GpCYP89 belongs to the CYP716 family. To date, approximately ten enzymes from the CYP716 family have been reported to be involved in triterpenoid oxidation reactions, such as CYP716A86 and CYP716A83 from Centella asiatica, CYP716A111 from Aquilegia coerulea and CYP716A141 from Platycodon grandiflorus, suggesting that CYP716 enzymes may play an important role in triterpenoid biosynthesis in eudicots [50]. GpUGT35 belongs to the UGT94 family. It was reported that, in this family, SgUGT94-289-2 from Siraitia grosvenorii can glycosylate mogroside M2-E and mogroside M3 at the C-24 position to yield mogroside M3x and mogroside Sia [51,52]. Further study is needed to elucidate the roles of GpCYP89 and GpUGT35 in gypenoside biosynthesis. TFs play an important role in response to biotic and abiotic stress, and may further regulate the biosynthesis of secondary metabolites [53]. In our study, we annotated 1362 TFs and further identified six GpWRKYs and one GpbHLH as candidates for the regulation of gypenoside biosynthesis. The six GpWRKYs all come from subgroup IIc. In this subgroup, it was reported that some WRKYs can regulate secondary metabolite biosynthesis. For example, CjWRKY1 in Coptis japonica is a positive regulator of berberine biosynthesis through increasing the level of transcripts of pathway genes [54]. The GpbHLH candidate belongs to subgroup IVa. In this subgroup, CrBIS1 in Catharanthus roseus is reported to increase the content of iridoids and MIAs in C. roseus by transactivating the expression of genes encoding the enzymes that catalyze geranyl diphosphate to iridoid loganic acid [55]. Further analysis is needed to verify the predicted roles of these WRKYs and bHLH in the regulation of gypenoside biosynthesis.

Conclusions
G. pentaphyllum is not only a promising gypenoside or ginsenoside provider but also an emerging medicinal model plant for triterpenoid biosynthesis. Here, we first combined PacBio sequencing and Illumina sequencing (hybrid sequencing) to mine the uncharacterized genes encoding enzymes involved in gypenoside biosynthesis in G. pentaphyllum. Highly accurate full-length transcripts produced by hybrid sequencing can significantly improve the accuracy of gene annotation and gene expression quantification, which will be helpful for gene discovery related to gypenoside biosynthesis. Hybrid sequencing first provided abundant genetic information resources for G. pentaphyllum, laying an important foundation for research on genetic breeding, metabolic regulation, synthetic biology and other aspects of G. pentaphyllum. In addition, combined phylogenetic and coexpression analyses identified one OSC, one CYP450, one UGT, six WRKYs and one bHLH from the G. pentaphyllum transcriptomes to be the lead candidates involved in gypenoside biosynthesis. These results will provide new insights into gypenoside biosynthesis.

Plant materials and sample preparation
Plant materials of G. pentaphyllum (two years old) were grown in the Institute of Medicinal Plant Development (Additional file 9: Figure S6). Roots, stems, leaves, and MeJA-treated leaves of G. pentaphyllum were collected. The induced leaves were sprayed with 200 μm MeJA and collected at different treatment times (0 h, 6 h, 12 h, 18 h and 24 h). After collection, all samples were immediately frozen in liquid nitrogen and stored at − 80°C until RNA was extracted.

RNA extraction
Total RNA was extracted using the RNAprep Pure Kit for plant (Tiangen Biotech, China) and quantified by Qubit (Invitrogen TM Life Technologies, USA). The RNA integrity was evaluated with an Agilent 2100 Bioanalyzer (Agilent Technologies, USA).

PacBio sequencing
The total RNA from different tissues was equally mixed together. PolyA RNAs were isolated from total RNA using Dynal oligo (dT) 25 beads (Life Technologies, USA) and subjected to Iso-Seq library construction. Briefly, cDNA was synthesized from the polyA RNA using the SMARTer ® PCR cDNA Synthesis Kit (Clontech, Mountain View, CA, USA). Size selection was carried out on a BluePippin (Sage Science, USA), and 0.8-2 kb, 2-3 kb, 3-6 kb fractions were collected. For each fraction, a SMRTbell template library was prepared and sequenced using the PacBio RSII platform. A total of eight SMRT cells were carried out in this study.

Illumina sequencing
Illumina sequencing libraries were constructed with the root, stem and leaves from G. pentaphyllum with three replicates. In addition, MeJA-treated leaves at different times (0 h, 6 h, 12 h, 24 h) were collected to construct Illumina sequencing libraries. These 13 transcriptome libraries were paired-end sequenced using the Illumina NextSeq 500 system.

Bioinformatics pipeline for unigene generation
The high-quality reads of insert (circular consensus sequences) were generated from SMRT raw data and classified into full-length non-chimeric reads and non-fulllength reads with SMRT Analysis software v2.3.0 [56]. ICE and Quiver were used to cluster and polish the fulllength non-chimeric sequences to generate the polished consensus sequences. Illumina clean reads were obtained after the adapters, and low-quality sequences were filtered out from the raw data. These clean data were assembled into Trinity unigenes with the Trinity assembler [57]. Finally, the polished consensus sequences and Trinity assembled unigenes were merged together, and redundant data were removed with CD-HIT software (-T 12 -M 45000 -c 0.85) to obtain final unigenes.

Phylogenetic analysis
Amino acid sequences were aligned using the ClustalW program, and evolutionary distances were computed using the pairwise deletion method, and a neighbor-joining (NJ) tree was constructed with MEGA6 [64]. Bootstrap values obtained after 1000 replications are given on the branches. Values less than 50% are not shown.

Qualitative analysis by real-time PCR analysis and heatmap visualization
RNA samples were isolated from roots, stems, leaves and MeJA-treated leaves (0 h, 6 h, 12 h, 18 h and 24 h). Reverse transcription was performed using the GoScript™ Reverse Transcription System Kit (Promega, USA). Primers were designed using Primer3 version 4.1.0 (http://primer3.ut.ee). The primers used in this study are listed in Additional file 10: Table S4. RT-PCR analyses were then conducted using a Bio-Rad CFX96 RT-PCR system. The reaction mixture (20 μL) contained 10 μL 2 × SYBR Premix Ex Taq (Takara, Tokyo, Japan), 0.5 μL each of the forward and reverse primers and 1 μL of template cDNA. PCR amplification was performed under the following conditions: 95°C for 30 s; 40 cycles of 95°C for 5 s, 60°C for 30 s and 72°C for 15 s; 95°C for 10 s. The heatmap was plotted with the pheatmap package in R [65]. Here, the upstream gene SE was used as a reference and the 2 -ΔΔCt values were plotted with scale = "row" enabled.

WRKY and bHLH domain analysis
The domains of WRKY and bHLH were analyzed by MEME (http://meme-suite.org). Then, the alignment of the conserved parts was trimmed out and submitted to Skylign (http://skylign.org) to generate their corresponding sequence logos.