For aseptic seedling culture, tuberous roots of R. glutinosa Libosch “Wen 85-5” were surface sterilized with 0.1% mercuric chloride and cultured on hormone-free MS agar medium (Murashige and Skoog, 1962). The MS propagation medium contained 30 g/L sucrose and 9 g/L agar. The seedlings were maintained at 26°C under a 14 h light/10 h dark photoperiod in a growth chamber. Tuberous roots of two R. glutinosa cultivars Wen 85-5 and QH were planted on April 20, 2015. The fields were maintained with locally standard production conditions. The tuberous roots (TR), stems (S), tender leaves (TL), unfolded leaves (UL), and senescing leaves (SL) of the Wen 85-5 plants were collected from six plants grown in the field for 4 months. The TR of QH plants were collected at the same time. All the tissues from the six plants were separately pooled for each determination. The experiments were repeated in triplicate.

The culture of Agrobacterium rhizogenes strain ACCC10060 was initiated from glycerol stock and grown overnight at 28°C with shaking (180 rpm) in liquid YEB medium to the mid-log phase (OD₆₀₀ = 0.8). Excised leaves of Wen 85-5 were dipped into A. rhizogenes ACCC10060 liquid inoculation medium for 8 min to infect leaf explants of 25-day-old seedlings, blotted dry on sterile filter paper, and incubated in the dark at 26°C on solid MS agar medium containing 100 μmol/L acetosyringone (AS). After co-cultivation for 2 days, the leaf explants were transferred to liquid MS agar medium containing 500 mg/L cefotaxime (Cef) and 100 μmol/L AS for hairy root induction. Numerous hairy roots were observed emerging from the wound sites after ~2–4 weeks. When the hairy roots had grown to a length of 2–3 cm, they were separated from the explants and cultured in dark at 26°C on MS agar medium with a gradual decrease of Cef at a 10 day interval to obtain bacteria-free culture. The hairy root segments were checked for bacterial contamination by culturing them in MS medium without Cef. The rapidly growing hairy root culture line was transferred to 50 mL of MS liquid medium in 100 mL flasks. The hairy root cultures were maintained at 26°C on a gyratory shaker at 120 rpm in the dark.

For confirmation, isolation of total DNA from bacterium-free R. glutinosa hairy roots and non-transformed (NT) roots was conducted using a modified CTAB method (Cuc et al., 2008). Polymerase chain reactions (PCR) were performed for rolB and rolC using specific primers recognizing rolB (5′-GCTCTTGCAGTGCTAGATTT-3′; 5′-GAAGGTGCAAGCTACCTCTC-3′) and rolC (5′-CTCCTGACATCAAACTCGTC-3′; 5′-TGCTTCGAGTTATGGGTACA-3′). The amplification system and procedure were the same as in Tong et al. (2012), except for the primer annealing temperature, which was 58°C.

A solution of 10 mmol/L SA (Sigma-Aldrich) was prepared by dissolving 0.0414 g SA powder in 30 mL sterile distilled H₂O. MeJA (Sigma, China) was dissolved in DMSO to a concentration of 100 mmol/L. AgNO₃ (Ag⁺, 0.3397 g) was dissolved in 20 mL distilled H₂O with a concentration of 100 mmol/L. Put (Sigma, China) was dissolved in sterile distilled H₂O to a concentration of 10 mmol/L. All the solutions were filter sterilized through 0.22 μm filters and added to cultures to the desired final concentrations. Experiments were conducted in a complete random block design with three replicates of each treatment and three concentration levels (high, medium, and low concentration) for each elicitor, generating a total of 39 experimental units (10 treatments, 13 × 3 replicates). The final concentrations were 10, 25, and 40 μmol/L for SA and 5, 15, and 25 μmol/L for MeJA, Ag⁺, and Put, respectively. These four stock solutions of elicitors were stored at 4°C prior to use. For elicitor treatments, 0.1 g (fresh weight) hairy roots were incubated in 50 ml liquid MS medium in 100 ml Erlenmeyer flasks. Sterilized solutions were individually added to 20-day-old hairy root culture media at the designed concentrations to investigate their effects on hairy root growth and acteoside content. Water was added as a negative control. After 10 days of treatment with different elicitors, hairy roots were harvested for determination of dry weight and acteoside content. After 0, 12, and 24 h treatment with SA, hairy roots were collected for determination of acteoside content and RNA-seq analysis. The 0, 3, 9, 12, and 24 h SA-treated hairy roots were used for expression analysis. All treatments were performed in triplicate, and the results were averaged.

For RNA-seq analysis, pooled hairy roots from three replicate cultures were used to prepare an RNA sample for each treatment. Frozen hairy root tissue (100 mg) was ground to a powder in liquid nitrogen. Total RNA was isolated with TRIzol® Reagent (Invitrogen) as described by the manufacturer's protocol, then treated with RNase-free DNaseI (Invitrogen). The RNA quality and quantity were determined with a Nanodrop™ 2000 spectrophotometer (Thermo Fisher Scientific, USA) and a Bioanalyzer 2100 (Agilent, USA). Library preparation and Illumina sequencing were performed using an Illumina HiSeq™ 2000 from BGI-Tech (Shenzhen, China; Project ID: F15FTSCCKF2309). Transcripts from tuberous roots and leaves transcriptomes of R. glutinosa (Project ID: F13FTSCCKF1467) were used as references for read mapping and gene annotation.

The raw reads were subjected to quality control (QC) analysis to identify high-quality sequencing data and clean reads. Filtered clean reads from each sample were then separately aligned to the reference transcripts of R. glutinosa using the Bowtie2 software (Langmead et al., 2009) and used to estimate the abundance of gene transcripts using the RSEM method (Li and Dewey, 2011), measured as fragments per kilobase of transcript per million fragments sequenced (FPKM) (Trapnell et al., 2010). Differentially expressed genes were identified based on the method described by Audic and Claverie (1997). We used a false discovery rate (FDR) of ≤0.001 and an absolute value of FPKM fold-change of ≥2 as the thresholds to evaluate the significance of differentially expressed transcripts (DETs). The expression profiles of DETs from different samples were analyzed by hierarchical clustering, and a heat map of expression values was generated using the T-MeV 4.9.0 software (Howe et al., 2011). Annotation analyses of DETs were performed using WEGO software (Ye et al., 2006) for GO term functional classification, and pathway enrichment analysis of DETs was performed based on the Kyoto Encyclopedia of Genes and Genomes (KEGG) database (Kanehisa et al., 2008).

Total RNA was extracted using TRIzol reagent (Invitrogen) according to the manufacturer's protocol. Total RNA samples were treated extensively with RNase-free DNase I (Invitrogen) to remove any contaminating genomic DNA. cDNA was obtained from 1 μg of total RNA using a PrimeScript™ II 1st Strand cDNA Synthesis Kit (TaKaRa Bio, Dalian). Quantitative real-time PCR (qRT-PCR) assays were performed using SYBR® Premix Ex Taq™ II (Tli RNaseH Plus) (Takara Bio, Dalian) on a Bio-Rad iQ5 Real-Time PCR System (Bio-Rad, USA) as described by Wang et al. (2015). The RgTIP41 gene (GenBank accession number KT306007) was used as the reference gene to calculate relative expression levels based on the 2^−ΔΔCt method (Schmittgen and Livak, 2008). Data from three biological replicates were analyzed using analysis of variance (ANOVA) followed by Student's t-test (p < 0.05). The primers used are listed in Table S1.

Dried and powdered R. glutinosa hairy roots (0.8 g) were extracted with methanol for 1.5 h at room temperature, and the weight loss was made up with methanol. After filtration, the 20 mL extract was vacuum evaporated to dryness and the residue was dissolved in 5 mL acetonitrile-acetic acid solution (16:84, v/v), and then was determined by HPLC. Targeted analysis of acteoside was performed on an Agilent 1200 HPLC with a C18 column (4.6 × 250 mm, 5 μm) at 30°C. The mobile phase was acetonitrile-acetic acid solution (16:84, v/v) and was run at 1 mL/min. The injection volume was 20 μL for each sample, and detector wavelength was 334 nm. The reference standard of acteoside was purchased from Chengdu Must Bio-Technology Co., Ltd (Sichuan, China). Acteoside concentration was calculated by interpolating the peak area with a calibration curve obtained from standard purified compounds in the range of 0.038–0.56 mg/mL and expressed as mg/g of hairy root dry weight or mg/L of liquid culture.

The hairy root lines of R. glutinosa were established by Agrobacterium rhizogenes infection, and its state was confirmed by PCR for the amplification of two segments of (423 and 626 bp) corresponding to the bacterial rolB and rolC genes, respectively (Figure S1). After 30 days of growth in liquid MS medium, the biomass of hairy roots had essentially reached a maximum (Figure S2), and the roots were used for acteoside determination. To investigate the effects of abiotic elicitors on acteoside accumulation, 20-day-old R. glutinosa transgenic hairy roots were treated with SA, MeJA, silver nitrate (Ag⁺) and Put and then cultivated continuously for 10 days (Figure 1A). Relative to control untreated hairy roots, SA or MeJA both induced a significant increase in acteoside content at an appropriate concentration (Figure 1B). However, SA appeared to be more effective in triggering acteoside accumulation than MeJA. An approximate concentration of 25 μmol/L SA was the most efficient and acteoside content under this treatment (11.66 ± 0.19 mg/g dw) was 2.28-fold that of the control. Because the final biomass of hairy roots after SA elicitation was substantially the same and did not differ from that of control hairy roots (Figure 1C), the total amount of acteoside extracted from 1 L of SA hairy root liquid culture (53.87 ± 9.91 mg/L) was significantly higher than the amount in untreated hairy roots (20.75 ± 2.15 mg/L) (Figure 1D). However, Ag⁺ inhibited acteoside accumulation, and the content of acteoside was significantly lower with increasing concentrations of Ag⁺ (Figure 1B). When 25 μmol/L Put was added, the biomass was increased by 33% relative to the control, while no significant increase of content and a final yield of acteoside was abserved (Figures 1C,D). It can be concluded that the best elicitor for further enhancement of acteoside accumulation in transgenic R. glutinosa hairy root line Wen 85-5 is SA at a concentration of 25 μmol/L.

Elicitors can trigger physiological responses and secondary metabolites accumulation in a short time (Sá et al., 1992). To study the short-term effect of SA treatment on acteoside accumulation, hairy roots were collected at 0, 12, and 24 h after SA treated. The result showed that acteoside accumulation were significantly improved in 12 and 24 h SA treated hairy roots compared to untreated hairy roots (Figure S3).

As the genome sequence of R. glutinosa is not yet available and we were most interested in transcribed genes, we conducted transcriptomic analysis of the induced R. glutinosa hairy root cultures. To investigate the transcriptome response to SA induction in the hairy roots of R. glutinosa, the Illumina HiSeq 2000 platform was used to perform high throughput tag-seq analysis on R. glutinosa hairy root libraries constructed at 3 time points before and during the 24 h SA treatment period. The samples were collected at 0, 12, and 24 h after SA treatment. The major characteristics of these libraries are summarized in Table 1. Approximately 13 million sequence tags per library were obtained. Prior to mapping these tag sequences to the reference transcriptome sequences, adaptor, low quality and single copy tags were removed, resulting in a total of more than 12.7 million clean sequence tags per library. The RNA-seq reads data of three samples have been deposited in NCBI's Sequence Read Archive (SRA) under accession number SRP103641, including SRR5438036, SRR5438037, and SRR5438042.

Using our assembled R. glutinosa unigenes as the reference transcriptome, 47.44–49.61% of the unique clean reads (6.03–6.34 million) were uniquely mapped, representing 57.37–58.39% of the total mapped reads. In the corresponding individual samples, 48,246, 48,377, and 49,009 reference unigenes were identified, indicating that our sequencing depth was sufficient to approach saturation.

To identify the genes with significant change in expression level during SA treatment, DETs between 0 h and the other two libraries (12 and 24 h) were identified. Our digital expression analysis identified 2,715–4,018 DETs with at least 2-fold difference in expression levels and a FDR < 0.001 during 24 h SA treatment (Figure 2A). The differential expression patterns among libraries revealed that the largest differences occurred between 0 and 24 h, and 2,401 up-regulated and 1,617 down-regulated DETs were identified at 24 h; whereas the smallest differences existed between 12 and 24 h, for which only 1,470 up-regulated and 1,245 down-regulated DETs were identified (Figure 2A). Furthermore, we could identify 272 common genes through a Venn diagram of the three comparison periods in SA treatment (Figure 2B). In detail, 90 common up-regulated genes were found among the three comparisons, including alcohol dehydrogenase gene (ALDH, CL6693.Contig36), 4-coumarate-CoA ligase gene (4CL, CL929.Contig4), and cinnamoyl-CoA reductase gene (CCR, Unigene7) (Figure S4A; Table S2), etc. In contrast, only 54 common down-regulated genes were found in the same comparisons (Figure S4B; Table S2).

To further validate the results of mRNA sequencing, 10 transcripts were randomly selected along with their specific primers for qRT-PCR analysis (Table S1). All 10 transcripts showed similar expression patterns to the in silico differential analysis results from DET sequencing. Expression levels measured by qRT-PCR were strongly correlated with those of DETs identified by mRNA-seq (r = 0.8182; Figure S5). In addition, the production of expected fragment sizes using designed primers also supported the reliability of the de novo assembly.

Gene ontology functional classification analyses were performed to classify the functions of the DETs during SA treatment. Based on sequence homology, all DETs could be categorized into 38 functional groups at time points 12 and 24 h (Figure 3). In the biological process (BP) category, the largest groups were metabolic processes, cellular process, and single-organism process. In the cellular component (CC) category, the greatest numbers of genes were found in the cell part and cell terms. In the molecular function (MF) category, most of the DETs were mapped into the catalytic activity and binding groups.

We conducted a KEGG-pathway-based analysis to obtain a better understanding of the biological functions of these DETs. Among these DETs with a KEGG pathway annotation, 3,504 DETs were identified in the 0 h-vs.-12 h library and mapped onto 120 KEGG pathways (Table S3), while 3,537 DETs were identified in the 0 h-vs.-24 h library and mapped onto 121 KEGG pathways (Table S4). To determine whether genes involved in secondary metabolites were enriched, the KEGG pathway database was searched using the DETs to reveal the top 40 significantly enriched pathways. The following KEGG pathways were enriched after the 12 h SA treatment: 127 DETs were enriched in “phenylpropanoid biosynthesis,” 439 DETs were related to “biosynthesis of secondary metabolites,” 73 DETs were related to “flavonoid biosynthesis,” and 41 DETs were related to “phenylalanine metabolism” (Figure 4, Table S3). All 14 DETs related to “phenylalanine, tyrosine and tryptophan biosynthesis” were up-regulated after SA treatment (Figure S6), suggest that SA may stimulate accumulation of the two precursors in phenylpropanoid biosynthesis. Furthermore, 18 DETs were found to be associated with the “tyrosine metabolism.” Among these, more than 88% of the DETs were up-regulated (Figure S7). Similar results were observed for the 12 h SA treatment sample (Figure 4, Table S4). These results indicate that phenylalanine and tyrosine biosynthesis pathways were up-regulated by SA treatment, the increased precursors maybe a cause of improving accumulation of acteoside.

Acteoside biosynthesis begins with the generation of phenylalanine and tyrosine precursors by the shikimate pathway (Alipieva et al., 2014). The hydroxytyrosol moiety of acteoside is synthesized from tyrosine through either tyramine and/or dopamine, whereas its caffeoyl moiety is synthesized from phenylalanine via a cinnamate pathway (Ellis, 1983; Saimaru and Orihara, 2010). The BAHD acyltransferase superfamily encodes proteins catalyzing the acyl-transfer from coenzyme A-activated acids to varying acceptor molecules (D'Auria, 2006). Therefore, we inferred that shikimate O-hydroxycinnamoyltransferase (HCT) maybe involved in acteoside biosynthesis. Using the annotated R. glutinosa transcriptome assembly, we identified 219 unigenes with sequence lengths of more than 200 bps, annotated as known enzymes involved in the pathway for acteoside biosynthesis (Figure 5A; Table S5). All of the genes encoding enzymes involved in the biosynthesis of the acteoside were present in our R. glutinosa transcriptome (Table 2). In most cases, more than one unigene was annotated as the same enzyme. Such unigenes may represent different fragments of a single transcript, different members of a gene family, or both.

To identify which of these genes are important for acteoside biosynthesis, we performed hierarchical clustering analysis. The 219 unigenes were placed into nine clusters based on their expression patterns in hairy roots treated with SA (Figure 5B). In cluster 6 and 9, gene expression significantly increased and reached the highest expression level at the 12 or 24 h time point after SA treatment (Figure 5B; Figure S8), indicating that the expression of cluster 6 and 9 genes coincides with acteoside accumulation. In clusters 5, 7, and 8, the expression level of unigenes was up-regulated only at the 12 or 24 h time points. In contrast, the genes in clusters 1, 2, and 4 were down-regulated at at least one time point. Oddly, the FPKM-values of genes in cluster 3 were always zero at the three time points, indicating these genes did not expressed in hairy roots. We noted that similar expression trends exhibited for some genes under both 12 and 24 h SA treatment were identical to the above-mentioned enzyme families commonly expressed in both genotypes, but the transcript expressed was different.

To identify genes displaying significant expression changes during SA treatment, DETs were analyzed by comparing 12 and 24 h libraries with the control library. A total of 54 unigenes were found to be differentially expressed at at least one of the time points after SA treatment (Table S6). Surprisingly, all DETs were significantly up-regulated genes, and none were significantly down-regulated. Forty-one and thirty DETs exhibited a significant increase in expression levels at 12 and 24 h time points, respectively. There were 17 up-regulated DETs in both the 12/0 h and 24/0 h comparisons. Among the 54 identified acteoside biosynthesis genes, 14 of them encode 5 enzymes predominantly expressed in the early stage (12 h) of SA treatment. They include five phenylalanine ammonia-lyase (PAL) genes, three copper-containing amine oxidase (CuAO) genes, and two each of the cinnamate-4-hydroxylase (C4H), coumaroylquinate (coumaroylshikimate) 3′-monooxygenase (C3H) and tyrosine decarboxylase (TyDC) genes. In contrast, most of the ALDH genes (8/10) had higher expression levels at the late stage (24 h) after SA treatment. Sixteen DETs encoding four enzymes involved in acteoside biosynthesis were significantly up-regulated at both 12 and 24 h time points. They include four 4CL genes, three polyphenol oxidase (PPO) genes, five UDP-glucosyl transferase (UGT) genes, and four Shikimate HCT genes.

To gain a more comprehensive understanding of the modes of action of SA, MeJA, Ag⁺, and Put in affecting the accumulation of acteoside in R. glutinosa hairy roots, the expression levels of genes encoding enzymes belonging to the biosynthesis pathway of acteoside were examined by qRT-PCR during the 24 h period after elicitation. RgTIP41 gene was stably expressed in roots, stems, leaves and flower tissues of R. glutinosa, and has been used as internal reference gene for research the expression characteristics of R. glutinosa coding genes (Sun et al., 2014; Wang et al., 2015). Gene expression levels were normalized using RgTIP41 as the reference gene as internal standard, and the transcripts of all associated genes in the control were set to 1. The results of gene transcript analysis are shown in Figure 6, indicating that all the selected genes investigated were up-regulated by SA treatment, although the time and degree of response differed. Under SA treatment, the expression levels of CL1389.Contig1 (PAL), CL2673.Contig1 (C4H) and Unigene9725 (C3H) were observed to gradually increase and reached peak values at 12 h with transcript levels 1.6-, 3.6-, and 7.0-fold higher, respectively, than those of the control. The relative up-regulated levels of CL1583.Contig3 (TyDC), Unigene35591 (ALDH), Unigene11579 (PPO), CL5641.Contig1 (UGT), and CL592.Contig1 (UGT) after 9 h SA treatment were the most significant, reaching 3.9-, 4.3-, 222.6-, 32.5-, and 23.6-fold, respectively. CL929.Contig4 (4CL), Unigene9959 (CuAO) and Unigene19512 (HCT) showed a similar expression trend: the largest difference between treatments and control occurred at 12 h after treatment, reaching 5.4-, 6.2-, and 3.8-fold, respectively. In addition, Unigene11579 (PPO) and CL5641.Contig1 (UGT) exhibited significantly higher transcriptional levels than other genes, indicating that they might be more sensitive and critical than other genes for acteoside biosynthesis in hairy roots following SA elicitation.

Under MeJA treatment, most of the investigated genes (10 genes) were up-regulated at at least one time point (Figure 6). The largest relative difference in the expression level of CL2673.Contig1 (C4H), Unigene9725 (C3H), CL929.Contig4 (4CL), CL1583.Contig3 (TyDC), Unigene35591 (ALDH), Unigene11579 (PPO), Unigene9959 (CuAO), and Unigene19512 (HCT) between treatments and control occurred at 3 h after treatment, reaching 3.8-, 4.9-, 4.6-, 8.9-, 7.0-, 10.6-, 5.7-, and 4.5-fold, respectively. Furthermore, CL592.Contig1 (UGT) was up-regulated to ~81.9- and 79.3-fold after being treated for 9 and 12 h, respectively. The expression profiles of these selected genes suggested that they may promote acteoside accumulation in hairy roots of R. glutinosa treated with MeJA. However, after Ag⁺ and Put treatments, the expression levels of the major selected genes showed only slight up-regulation at some time points or were even down-regulated, indicating that Ag⁺ and Put did not induce the transcription of these genes to increase acteoside accumulation in hairy root of R. glutinosa.

To illuminate the relative expression of enzyme-encoding genes with acteoside accumulation, we detected relative expression levels of the 11 candidate genes in the tuberous root (TR), stem (S), tender leaf (TL), unfolded leaf (UL), and senescing leaf (SL) of R. glutinosa by qRT-PCR. This analysis revealed that CL1389.Contig1 (PAL), CL2673.Contig1 (C4H), Unigene9725 (C3H), CL929.Contig4 (4CL), CL1583.Contig3 (TyDC), Unigene11579 (PPO), Unigene9959 (CuAO), CL5641.Contig1 (UGT), CL592.Contig1 (UGT), and Unigene19512 (HCT) showed higher expression levels in the leaf (including TL, UL, and SL) than in the tuberous root (Figure 7), consistent with the acteoside content in the tuberous root, which was significantly lower than that of the leaf for R. glutinosa cultivars Wen 85-5 and QH (Figure S9A).

The secondary metabolite content in many medicinal plants has been found to vary significantly with genotype (Alagna et al., 2012; Song and Li, 2015), reflecting genetic control by the plant. QH is a high-acteoside (HA) R. glutinosa cultivar with a higher content of acteoside (up to 1.61 mg/g dw) in tuberous roots at three harvest times (Figure S9B). In contrast, Wen 85-5 is a low-acteoside (LA) cultivar with a lower content of acteoside (lower than 0.69 mg/g) in tuberous roots (Figure S9B). To detect the expression levels of the 11 candidate genes in Wen 85-5 and QH tuberous roots, qRT-PCR was performed. The results showed that most of candidate genes, including CL1389.Contig1 (PAL), CL2673.Contig1 (C4H), CL929.Contig4 (4CL), CL1583.Contig3 (TyDC), Unigene9959 (CuAO), CL5641.Contig1 (UGT), CL592.Contig1 (UGT), and Unigene19512 (HCT) had higher expression levels in QH than in Wen 85-5 (Figure 8). The identification of candidate genes highly expressed in high acteoside cultivars may provide an insight into the relationship between genotype and acteoside biosynthesis.

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.