PTHRCs were established as reported in our previous study (Wang et al., 2015). Briefly, sterile healthy P. tunicoides plants were maintained on MS media supplemented with 20 g/L sucrose (pH 5.8) at 25 °C and a 16-h photoperiod. Excised leaves of P. tunicoides were used to induce hairy roots by infection with A. tumefaciens strain ATCC15834. The hairy roots were routinely subcultured and maintained in 100 mL liquid MS media supplemented with 30 g/L sucrose (pH 5.8) in conical flasks (250 mL) in a shaker with a speed of 100 rpm for 35 days. The growth conditions for suspension cultivations were 25 ± 1°C in darkness.

SA (Sigma-Aldrich, United States) was dissolved in ultrapure water as a stock solution (10 mg/mL) and filter-sterilized. After 35 days of suspension cultivation, hairy roots at their maximum weight were used to perform an SA elicitation experiment. The subcultured media were renewed by 100 mL fresh liquid media supplemented with 30 g/L sucrose (pH 5.8), in which the SA stock solution was added to give the final concentrations of 5, 10, 20 mg/L SA. The same volume of ultrapure water was added as control (CK) cultures.

The treated (5, 10, 20 mg/L SA) and CK hairy root groups were sampled at 1, 3, 5, 7, and 9 days for saponin extraction and quantification. All experiments were repeated at least three times. The samples were vacuum freeze-dried for 12 h (Telstar, Spain), and ground into powders. The powders (1 g) were extracted using 80% ethanol with ultrasonication for 100 min. After extraction with 40 mL of n-butanol three times, the combined extracts were rotary evaporated at 40 °C (BUCHI, Switzerland), and redissolved in 10 mL of methanol.

Total saponin content was determined by vanillin-perchloric acid colorimetry according to a previous study (Xiang et al., 2001), and oleanolic acid was used as the reference standard. Levels of two representative saponin derivatives [quillaic acid (QA) and gypsogenin (Gyp)] were quantified by high-performance liquid chromatography (HPLC) based on a previous study (Kolodziej et al., 2018). Chromatographic analysis was carried out on a Shimadzu series LC-20 AD XR instrument, with an SPD-M20A diode array detector, on a reverse-phase Sunniest C18 (ChromaNik, Japan) analytical column (250 mm × 4.6 mm, 5 μm) at 25°C. The mobile phase consisted of acetonitrile (solvent A) and 0.1% (v/v) H₃PO₄ aqueous solution (solvent B) using the gradient elution as follows: 50% A 0–2 min, 50–100% A 2–35 min, 100% A 35–40 min. The injection volume was 10 μL. The detection wavelength was 210 nm, and the flow rate was 1 mL/min. Gyp and QA standards (Yuanye, China) were used as standard calibration curves.

For transcriptome sequencing, antioxidant enzyme and qRT-PCR analysis, the 35-day-old hairy roots were harvested from 5 mg/L SA treatment at 8 (SA_8 h) and 24 h (SA_24 h) after initiating treatments. Untreated hairy roots were collected before SA treatment, and indicated as controls (SA_0 h). Total RNA was extracted using TRIzol reagent (Invitrogen, United States) according to the manufacturer’s instructions. RNA purity was checked using the NanoPhotometer^® spectrophotometer (IMPLEN, United States). RNA concentration was assessed using a Qubit^® RNA Assay Kit in a Qubit^® 2.0 Fluorometer (Life Technologies, United States).

Equal amounts of RNAs from different hairy root groups (SA_0h, SA_8h, and SA_24h) were pooled to provide the total RNA of P. tunicoides. The mRNA was enriched using the Oligo d(T) magnetic beads, then reverse transcribed to cDNA using Clontech SMARTer PCR cDNA Synthesis Kit (Clontech, United States). Amplification of double-stranded cDNA was followed by size selection using the BluePippin system (Sage Science, United States), and fragments of 1–6 kb were retained. The full-length cDNA was generated, and the cDNA ends were repaired and ligated to sequencing adapters. SMRTbell template libraries were obtained and subsequently sequenced on the PacBio Sequel RS sequencing instrument.

Raw sequence data obtained from SMRT sequencing were processed using SMRTlink6.0 software (Chin et al., 2013). Circular consensus sequence (CCS) reads were generated from subread BAM files, and classified into full-length non chimera (FLNC) reads, full-length chimeric reads, non–full-length (NFL) reads, and short reads based on poly(A) signal, 5’ and 3’ adaptor check. NFL and full-length FASTA files were then fed into the cluster consensus isoforms by isoform-level clustering (ICE) (Gordon et al., 2015). Finally, the redundancies were removed using CD-HIT-EST to obtain unigenes (Huang et al., 2010).

RNA samples from three SA-treated hairy roots (SA_0h, SA_8h and SA_24h) were used for Illumina library construction and sequencing. Each group had three biological replicates. RNA-seq libraries were generated using NEBNext^® Ultra^TM RNA Library Prep Kit for Illumina^® (NEB, United States) following the manufacturer’s recommendations and index codes were added to attribute sequences of each sample. The libraries were constructed with insert sizes of 250–300 bp in length. Clustering of the index-coded samples was performed on a cBot Cluster Generation System using TruSeq PE Cluster Kit v3-cBot-HS (Illumina). After cluster generation, the libraries were sequenced on an Illumina Hiseq Xten platform from Novogene Experimental Department (Beijing, China) and the raw data were generated. The clean data were obtained by removing reads containing adapters, poly-N and low quality reads from raw data using Trimomatic (v0.36) (Bolger et al., 2014). The Q20, Q30 and GC-content of the clean data were calculated. The raw sequence data have been deposited in the Genome Sequence Archive¹ under accession number CRA002795.

The gene function was annotated, and classifications were based on NCBI non-redundant protein sequences (Nr)², NCBI nucleotide sequences (Nt; see text footnote 1), protein family (Pfam)³, eukaryotic/clusters of orthologous groups (KOG/COG)⁴, Swiss-Prot⁵, Kyoto Encyclopedia of Genes and Genomes (KEGG)⁶ and Gene Ontology (GO)⁷ databases with local BLAST programs (E-value < 1.0E^–5). Hmmscan was adopted for Pfam annotation⁸, and Blast2GO was used for GO annotation (Conesa et al., 2005). The ANGEL software was used to predict coding sequences (CDS)⁹. The iTAK software was applied to identify the TF family database (Zheng et al., 2016).

The gene expression levels were determined by RSEM software with bowtie2 (mismatch 0) (Li and Dewey, 2011). Illumina clean data were mapped onto the SMRT sequencing data, and the readcount for each gene was obtained from the mapping results. Then the unique readcounts for each transcript were normalized by calculating fragment per kilobase per million (FPKM). DEGs of different libraries were analyzed using the DESeq R package (1.10.1) (Anders and Huber, 2010). Genes with an adjusted P-value < 0.05 and | log2FoldChange| > 1 found by DESeq were assigned as differentially expressed. GO enrichment analysis of the DEGs was implemented by the GOseq R packages based Wallenius non-central hypergeometric distribution (Young et al., 2010). KEGG enrichment analysis of DEGs was done using KOBAS software (Xie et al., 2011). Heatmaps were drawn by TBtools software (Chen et al., 2020). A phylogenetic tree was generated by the neighbor-joining method in MEGA7 (1,000 bootstrap replicates) (Kumar et al., 2016) and iTOL v5¹⁰. The protein sequences were obtained from UniProt database¹¹.

Hairy root samples were fully ground into 0.2 M phosphate buffer solution (pH 7.4) at 4°C, then the homogenates were centrifuged at 12,000 rpm at 4°C for 10 min. The supernatants were collected for peroxidase (POD) and glutathione reductase (GR) activities measurement according to the instructions of the assay kits (Nanjing Jiancheng, China).

Total RNA was extracted using RNAprep Pure Plant Kit (Tiangen, China) according to the manufacturer’s protocol. cDNA was obtained from 1 μg of total RNA using a TransScript One-Step gDNA Removal and cDNA Synthesis SuperMix (TransGen Biotech, China). qRT-PCR assays were performed using SuperReal PreMix Plus kit (Tiangen, China) with SYBR Green method on an ABI 7500 Real-Time PCR System (Thermo Fisher Scientific, United States). The Ptβ-Actin was used as the reference gene to normalize the relative expression levels via 2^–ΔΔCt method (Li et al., 2016). All qRT-PCR experiments were performed in three biological replicates. The primers were designed using Primer Premier 5.0 (Lalitha, 2000), and listed in Supplementary Table 1.

To investigate the elicitation effects of SA on triterpenoid saponin accumulation, 35-day-old P. tunicoides hairy roots were treated with a series of concentrations (5, 10, 20 mg/L) of SA, then cultivated continuously for 9 days. In comparison with the CK, the total saponin levels of PTHRCs were significantly increased by all groups of SA elicitation (Figure 1A). The treatments with 5 mg/L SA were most efficient, triggering the total saponin content of PTHRCs to peak at 2.49-fold higher than that of controls after 1 day (6.96 ± 0.43 mg/g DW). Hence, SA acted as an effective elicitor to positively modulate triterpenoid saponin production in PTHRCs, and 5 mg/L SA might be the optimal concentration for further analysis.

Previous studies illuminated that the molecular structures of P. tunicoides saponins were dominated by two representative aglycones: QA and Gyp (Cheikh-Ali et al., 2019). Furthermore, the trends of QA and Gyp with 1 day of SA treatments were also considered in acidolysis extraction of PTHRCs using HPLC (Supplementary Figure 1). Compared with control, levels of QA (3.14 ± 0.21 mg/g DW to 5.23 ± 0.23 mg/g DW) and Gyp (0.77 ± 0.02 mg/g DW to 0.83 ± 0.03 mg/g DW) in SA-elicited PTHRCs (5, 10 or 20 mg/L SA for 1 day) were increased by 35.4–125.9% and 20.1–28.4%, respectively (Figure 1B). Taken together, SA promoted yields of the representative saponin aglycones in PTHRCs in a relatively short time.

To investigate how SA triggered the transcript pathway of triterpenoid saponin biosynthesis in PTHRCs, we conducted transcriptome analysis using combined Illumina- and PacBio SMRT-based sequencing. Nine RNA samples from the different time points of SA treatments (0, 8, and 24 h) were sequenced using the Illumina Hiseq Xten platform. A total of 487,691,344 raw reads ranging from 48 to 62 million for each sample and 475,534,604 clean reads (97.5% of raw reads) ranging from 47 to 61 million for each sample were obtained. The quality of clean reads was good with Q20 > 96%, Q30 > 91.8% (except SA_0h_3) and GC 42–44% (Supplementary Table 2). To obtain wide coverage of the P. tunicoides transcriptome, pooled samples from different time points of SA treatments (0, 8, and 24 h) were sequenced using the PacBio RS II platform. A total of 39.43 Gb raw reads and 15,034,708 subreads were yielded. Then, 430,117 CCS reads were obtained after processing raw sequencing data, including 325,905 (75.8%) FLNC reads, and 100,217 (23.3%) NFL reads (Figure 2A). The CCS reads contained a mean length of 3,172 bp with N50 of 3,530 bp, while FLNC reads showed an average length of 2,818 bp with N50 of 3,168 bp (Figures 2B,C). After correcting with Illumina reads and removing the redundant sequences via CD-Hit analysis, 35,262 non-redundant transcript isoforms were generated.

We found that the unigene number obtained from the PacBio transcriptome (16,375) was much lower than that of the Illumina transcriptome (182,277) (Table 1). However, PacBio unigenes showed a significantly longer average length and larger N50 values when compared with Illumina unigenes (Table 1). We found that nearly 26.4% of Illumina unigenes had a length of < 500 bp; whereas only 0.3% of PacBio unigenes had lengths < 500 bp, and approximately 96.9% genes from SMRT sequencing were over 1,000 bp (Table 1). In addition, 65.1% of unigenes from the PacBio platform had complete CDSs, while only 13.5% of unigenes in the Illumina transcriptome had full-length CDSs (Table 1).

Seven databases, including Nr, Nt, Pfam, KOG/COG, Swiss-Prot, KEGG and GO, were used to annotate the functions of all unigenes. A total of 15,969 unigenes (97.5%) were annotated from at least one database, while 8,154 annotated unigenes were found in all databases (Figure 3A). For GO analysis, “metabolic process” (5,779 unigenes) showed the most enrichment pathways in biological process. “Cell” and “Binding” represented the major groups in the category of cellular component and molecular function, respectively (Figure 3B). For KOG categorization, “general function prediction only” (2,055 unigenes) was the largest among 26 functional groups, followed by “signal transduction mechanisms” and “post-translational modification, protein turnover, chaperones” (Figure 3C). A total of 442 unigenes were classified in the “secondary metabolites biosynthesis, transport and catabolism” group. For KEGG classification, the metabolism group exhibited the most genes (Figure 3D). In this group, 168 unigenes were found to represent “xenobiotics biodegradation and metabolism.” Only 376 unigenes were related to secondary metabolism, including a large number of genes (22.9%) associated with terpenoid metabolism (Supplementary Table 3).

lncRNAs are key functional regulators in plant biological processes (Long et al., 2017). In this study, 569, 1,324, 3,001, and 2,223 candidate lncRNAs were identified by the Coding-Non-Coding-Index (CNCI), Coding Potential Calculator (CPC), Pfam, and PLEK databases, respectively. A total of 4,678 lncRNAs were identified from at least one database. The full lengths ranged from 1,500 to 4,000 bp, with an average length of 2,656 bp (Supplementary Figure 2).

To investigate transcript changes in PTHRCs after SA treatments, DEG profiles were determined using the FPKM mapped reads method. The DEG numbers of three comparisons (“SA_8h vs. SA_0h,” “SA_24h vs. SA_0h” and “SA_24h vs. SA_8h”) are shown in Figures 4A,B. Compared to controls, a total of 775 up-regulated and 811 down-regulated unigenes were, respectively, observed after 8 and 24 h of SA treatments in PTHRCs. Only 91 or 64 genes were up- or down-regulated in all three comparisons, respectively. The up-regulated genes involved in “SA_8h vs. SA_0h” or “SA_24h vs. SA_0h” were further assigned in KEGG pathway analysis (Figures 4C,D). We found a series of up-regulated genes related to secondary metabolism, including terpenoid, flavonoid, phenylpropanoid and glutathione. Notably, the pathway of “monoterpenoid biosynthesis” reached the peak of factors among all KEGG pathways in both comparisons.

To illuminate how the SA signaling network functions in P. tunicoides, DEGs involved in SA perception or transduction pathways were analyzed based on transcriptome data. In PTHRCs, exogenous application of SA increased the expression levels of two NPR genes, named PtNPR1 (22,821/f3p0/2,141) and PtNPR4 (18,761/f3p0/2,457), which are widely reported as SA-binding proteins in model plants (Innes, 2018; Figure 5). SA-binding protein 2 (SABP2) was also up-regulated after SA elicitation (Figure 5). Furthermore, a series of genes encoding reactive oxygen species (ROS) scavenging enzymes were found in the transcriptome, including PODs and GRs. SA continuously up-regulated four PtPODs and one PtGR, as well as promoting enzyme activities of POD and GR at 8 h elicitation (Figure 5 and Supplementary Figure 3).

A total of 1,129 TF unigenes from 29 families were found in the PTHRC transcriptome, and the three largest TF families were related to C3H (89), NAC (56) and WRKY (56) (Figure 6A). Compared to controls, 275 or 300 DEGs associated with TFs were induced after 8 h or 24 h of SA treatments, respectively. Among them, 42 DEGs were continuously up-regulated with SA elicitation in both 8 and 24 h time points, of which several TF families generally participate in stress response and secondary metabolism (Endt et al., 2002), such as NAC (9), WRKY (4), bZIP (2), HSF (2), and MYB (1) families (Figure 6B). In particular, SA increased gene expression of PtWRKY70 (31,348/f5p0/1,467), PtMYB4 (33,896/f3p0/1,108), and PtWRKY40 (33,257/f4p0/1,251), suggesting that these unigenes might participate in SA-elicited secondary metabolism.

In the PTHRC transcriptome, six unigenes were identified as putative genes in the MVA pathway, including one AACT (acetyl-CoA acyltransferase), one HMGS (3-hydroxy-3-methylglutaryl-CoA synthase), two HMGRs, one PMK (phosphomevalonate kinase) and one MVD (mevalonate 5-diphosphate decarboxylase). In addition, 17 putative unigenes were recognized in the MEP pathway, containing eight DXSs (1-deoxy-D-xylulose-5-phosphate synthase), two DXRs (1-deoxy-D-xylulose 5-phosphate reductoisomerase), one MDS (2-C-methyl-D-erythritol 2,4-cyclodiphosphate synthase), four HDSs (4-hydroxy-3-methylbut-2-en-1-yl diphosphate synthase) and two HDRs (4-hydroxy-3-methylbut-2-enyl diphosphate reductase). Furthermore, 31 unigenes were putatively related to formation of β-amyrin, including one GPPS (geranyl diphosphate synthase), two FPSs (farnesyl pyrophosphate synthase), one SS (squalene synthase), 22 SEs (squalene epoxidase) and six bASs (β-amyrin synthase). Information for all of these genes is present in Figure 7 and Supplementary Table 4. Among them, a total of 35 DEGs were found after SA treatments in PTHRCs. We observed that SA significantly promoted the expression levels of one AACT, one HMGS and three DXSs (Figure 7). These enzyme genes were all located at the early stage of MVA and MEP pathways. Moreover, unigenes including one GPPS, 16 SE and one bAS were rapidly up-regulated at 8 h of SA treatment. Five SEs showed at least 10-fold up-regulation compared to control (Figure 7), suggesting that these PtSE genes might be the key regulatory factors in response to SA elicitation in PTHRCs.

Cytochrome P450s (CYPs) function in site-specific oxidization of the oleanane skeleton (Banerjee and Hamberger, 2018). A total of 114 putative CYP unigenes were recognized in transcriptome data, which belong to diverse CYP subfamilies like CYP71A, CYP72A, CYP716A, CYP89A, and CYP94A. Given that CYP72A and CYP716A subfamilies are primarily involved in the biosynthesis of triterpenoid saponin diversification (Ghosh, 2017), we found 29 genes belonging to CYP72A subfamilies, while only two genes belonged to CYP716A families (Supplementary Table 4). In total, 21 CYP72A genes were significantly up-regulated after 8 h of SA applications, and PtCYP72A219 (22,782/f10p0/2,106) showed the largest increase compared to controls (206.8-folds) (Figure 7 and Supplementary Table 4). Phylogenetic trees were inferred for PtCYP72A219 and eight known CYP72As that participated in terpenoid metabolism, and PtCYP72A219 showed the highest similarity (54.7%) to MtCYP72A58 in Medicago truncatula (Figure 8).

To validate the RNA-Seq transcriptome results, 10 candidate DEGs involved in the SA signaling network and saponin biosynthesis in PTHRCs were selected for measurement of transcript levels by qRT-PCR analysis, including PtNPR1, PtNPR4, PtAACT, PtHMGS, PtDXS, PtSE, PtbAS, PtCYP72A219, PtNAC29, and PtWRKY70. The fold change of qRT-PCR and RNA-seq data were compared. As shown in Figure 9, the relative expression results of these 10 candidate DEGs were generally consistent with the RNA-seq values. In particular, the expression levels of PtSE, PtbAS, PtCYP72A219, PtNAC29, and PtWRKY70 increased significantly at the early stage of SA treatments.