Comparative transcriptomics characterized the distinct biosynthetic abilities of terpenoid and paeoniflorin biosynthesis in herbaceous peony strains

Plant materials for RNA-Seq and HPLC

The P. lactiflora Pu–Bang, Xian–Tiao, and Guan–Shang strains were conserved and cultivated under field conditions at Bozhou University, Bozhou, China. The shoots and roots of 3-year old plants were isolated. To avoid circadian effects, we harvested all the tissues in the afternoon of the same day. The samples for RNA-seq with three biological replicates were frozen in liquid nitrogen immediately after harvesting. The isolated samples and purified RNA were stored at −80 °C. To measure paeoniflorin level by HPLC, roots of 3-year old plants were dried after removing its barks.

RNA extraction, library construction, and Illumina sequencing

The total RNA of individual samples was extracted and purified with the RNeasy® Plant Mini Kit (QIAGEN, Hilden, Germany). RNA concentration was measured using a Nanodrop 2100 spectrophotometer. RNA Integrity values were checked using an Agilent Bioanalyzer. The samples with a RIN score >8.5 were used for library construction (Liu et al., 2014). The sequencing libraries were generated using a NEB Next Ultra RNA Library Prep Kit for Illumina (New England Biosystems, Waltham, MA, USA), following the manufacturer’s recommendations. Library sequencing was performed on a Hiseq X10 system with 150-cycle paired-end sequencing protocol (Illumina, San Diego, CA, USA).

Bioinformatics analysis of RNA-seq datasets

RNA-seq datasets were checked using FastQC (Brown, Pirrung & McCue, 2017). Referring the methods used in recent studies with modification (Bedre et al., 2016; Lu et al., 2018), we assembled the transcriptome using Trinity (Haas et al., 2013). We aligned read sequences using HISAT2 (Kim, Langmead & Salzberg, 2015). Read count of each gene was called using HTseq-count (Anders, Pyl & Huber, 2015). Fragments per kilobase of exon per million fragments mapped of assembled transcripts (FPKM) were calculated and normalized using DESeq2 with global normalization parameters (Anders, Pyl & Huber, 2015; Love, Huber & Anders, 2014; Quinn & Chang, 2016; Zhang et al., 2014). The sequences of the assembled unigenes were annotated by Trinotate (Haas et al., 2013). Coding regions of unigenes were predicted using Transdecoder (Haas et al., 2013). BLAST v2.7.1 was performed to determine the sequence homology (e-value cutoff of 1e⁻⁵) to UniProt/Swiss-Prot, HMMR v3.1b2, EggNOG v4.5.1, and metabolic pathways were analyzed using KEGG database (Kanehisa et al., 2016). For the unigene with annotated homolog, we used the criterion of FPKM more than 0.8 to filter out the lowly expressed genes; for those without homolog, we used FPKM more than 1 as the cut-off criterion. Differential expression analysis was carried out using DESeq2 (Anders & Huber, 2010). Genes with normalized fold-change greater than 2, significance P-value less than 0.05, and Benjamini–Hochberg false discovery rate less than 0.1 were considered to be differentially expressed genes.

Quantitative detection of paeoniflorin

We referred to the previous method in order to measure paeoniflorin in samples (Yuan et al., 2013). The dried samples (0.50 g) were weighed and extracted with 50 mL of 50% aqueous methanol with ultrasonication for 30 min. The extracted samples were diluted with 50 mL 50% aqueous methanol and filtered with a 0.45-μm Millipore filter membrane (Millipore, MA, USA) at 25 °C. We used the Agilent 1200 LC Series (Agilent Technologies, Palo Alto, CA, USA) High Performance Liquid Chromatography (HPLC) system to measure the paeoniflorin abundance. The wavelength was set at 230 nm with a flow rate of 1.0 mL/min at a temperature of 25 °C. Standard compounds were purchased from the National Institutes for Food and Drug Control and the linearity of the standard compounds was checked at seven concentration solutions.

Quantitative RT-PCR

A total of 1 to 2 µg RNA samples were treated by DNase I (RNeasy plant mini kit) and were reverse transcribed with oligo (dT) primer and SuperScript III (Invitrogen, Carlsbad, CS, USA). cDNA samples were analyzed using quantitative PCR with SYBR Premix Ex Taq (Takara, Shiga, Japan) and a Biorad CFX96 real-time PCR system. The P. lactiflora Actin transcript sequence (JN105299) was used as endogenous reference genes to design primers to normalize the expression levels among samples (Qi et al., 2018; Yuan et al., 2013). The qRT-PCR reactions were carried out with two-step cycles (5 s 95 °C denaturation, 30 s 60 °C annealing and extension) and 45 cycles of amplification (Fig. S1A). The melting curves of primers were checked to ensure the primer efficiencies (Fig. S1B–S1L). We used three technical replicates to produce the average expression levels of the genes relative to that of the reference gene using the 2^−ΔΔCT method (Livak & Schmittgen, 2001). The primers are listed in Table S1.

RNA-Seq and de novo assembly of the P. lactiflora transcriptome

The P. lactiflora Pu–Bang (PB) and Xian–Tiao (XT) accessions are the most widely used herbaceous stains for medical uses due to their high levels of paeoniflorin, which is derived mainly from roots of 3-year old plants; whereas the Guan–Shang (GS) accession contains less paeoniflorin and is usually cultivated for ornamental uses. The morphological features of the 1- and 3-year old plants of the three strains were shown in Fig. 2, respectively. We determined the paeoniflorin levels of the isolated shoot and root samples of 3-year old plants using High Performance Liquid Chromatography (HPLC). The results confirmed the accumulated levels of paeoniflorin in PB and XT roots compared with those of GS (Fig. 3A).

To systematically identify genes and explore the gene expression network underlying paeoniflorin biosynthesis in P. lactiflora, we purified the RNA samples with three biological replicates derived from the shoots and roots of 3-year old PB, XT, and GS plants, and carried out the RNA-sequencing analysis using the Illumina paired-end 150 bp protocol. After filtering out the low quality reads, we obtained 775.73 million reads in total (Table S2). Using Trinity (Haas et al., 2013), we performed de novo transcriptome assembly and obtained 36,264 unigenes encoding 72,910 transcripts with 986 nt contig N50 length and 42.7% average GC content (Table 1). We also checked the completeness of de novo assembled Unigenes by aligning the sequences with Swiss-Prot database (The UniProt Consortium, 2017). The results showed that 7,131 (50%) out of the 14,165 aligned sequences have more than 40% coverage of known transcripts (Tables S3 and S4).

Functional annotation of expressed genes in the three P. lactiflora accessions

Due to embryo abortion and vegetative propagation, the P. lactiflora accessions have been undergoing reproductive isolation with a long cultivation history in the Bozhou region and were thought to be genetically distinct from each other (Zhou et al., 2002). However, the genomic evidence supporting this point is still lacking. We measured the expression levels of unigenes by calculating normalized Fragments Per Kilobase of exon per million fragments Mapped (FPKM). The unigenes with an FPKM value higher than 1 in at least one were used to perform hierarchical clustering analysis based on the Pearson correlation efficiency. We analyzed the hierarchical structure of the gene expression levels on a genome-wide basis (Fig. S2A). Most of the biological replicate samples belonged to the same clusters. Principal component analysis results also showed a similar result confirming a high level of reproducibility of the biological replicate samples (Fig. S2B). However, the tissues and strains were distributed in different cluster clades. It was noted that all the clades derived from PB and XT were separated from those of GS, indicating that the strains for medicinal uses are genetically divergent from the strains used for ornamental purposes, possibly due to the selection and reproductive isolation among the strains.

We predicted the protein-coding potential for the unigenes using the Transdecoder and searched for the annotation for the unigenes by aligning the assembled transcripts and predicted peptide sequences to the protein sequences annotated by Swiss-Prot, TrEMBL, Pfam and KEGG databases using Trinotate, BLASTP and BLASTX (Boutet et al., 2016; Camon et al., 2003; El-Gebali et al., 2019; Haas et al., 2013; Kanehisa et al., 2017). In total, we identified 34,203 unigenes containing significant matches to the annotated genes/proteins in at least one database (Fig. 3B). Of them, 28,083 (82%) were reproducibly detected by at least 2 data resources. The annotation information, predicted protein sequences, and FASTA-formatted sequences of these genes were provided in Supplemental Materials that could serve as a reference annotation for future studies (Table S3).

Transcription factors (TFs) with DNA-binding domains are the major regulators controlling the activity and specificity of the gene transcription process. We predicted genes encoding TFs in our assembled unigene dataset and identified 521 TF-encoding unigenes belonging to 32 TF families (Fig. 3C). Of these, AP2/ERF, bHLH, FAR1, bZIP and HB are the most abundant TF genes. The detailed information of the TF genes was provided in Table S5.

Identification of tissue- and/or strain-differentially expressed unigenes

Next, we searched for the differentially expressed unigenes. For each strain, we compared shoots and roots and extracted lists of root/shoot specific upregulated genes (fold change of expression level >2 and false discovery rate <0.05). Our analysis identified 10,125 up-regulated unigenes in roots and 11,911 in shoots. There were 1,906–3,737 unigenes specifically up-regulated in the roots and/or shoots of each strain; whereas only 332 (3%) and 886 (7%) unigenes were up-regulated in the roots and shoots of all the three P. lactiflora strains, respectively (Fig. 4). This result showed that a large number of the tissue-specific genes were also strain-specific, which is consistent with the fact that the three P. lactiflora strains have been undergoing genetic separation and selection during the last several centuries.

We performed the correlation analysis between the paeoniflorin levels and the expression levels of all the genes in shoots and roots of the three strains (Table S3). Our analysis identified 3,952 genes with the Pearson correlation coefficients higher than 0.6. Among the 161 gene with Pearson correlation coefficients higher than 0.9, we found several genes encoding transcription factors, E3 ubiquitin-protein ligase and oxidation-reduction related process, such as GT-2, ABCF4, SKP2A, FREE, GOR, TLP, RMA1H1, HAT5, XBAT31 and DELLA, suggesting transcriptional and post-translational regulation may be highly associated with paeoniflorin accumulation.

Identification of genes in terpenoid and paeoniflorin biosynthesis pathway

To further dissect the regulation pathways of the unigenes, we analyzed the gene list enrichment using the Kyoto Encyclopedia of Genes and Genomes (KEGG) datasets and KOBAS3.0 with hypergeometric testing and Benjamini and Hochberg correction (Kanehisa et al., 2012; Xie et al., 2011). In total, we identified 71 significantly enriched KEGG pathways in P. lactiflora (Table S6). The list included several KEGG terms of secondary metabolite biosynthesis, such as: Terpenoid backbone biosynthesis (P-value < 0.0152), Glycerophospholipid metabolism (P-value < 0.0001), Inositol phosphate metabolism (P-value < 0.0001), Pyruvate metabolism (P-value < 0.0002), Seleno compound metabolism (P-value < 0.0298), Ascorbate and aldarate metabolism (P-value < 0.0172) and Butanoate metabolism (P-value < 0.0352). P. lactiflora are generally known to have abundant secondary metabolites (Li et al., 2016a; Liu et al., 2017; Ma et al., 2016). Our transcriptome and pathway enrichment results are consistent with the metabolism profiling studies.

The previous studies have identified 19 EST sequences in the terpenoid backbone biosynthesis pathways in P. lactiflora, including 7 with full-length cDNA sequences (Yuan et al., 2013). However, the genes in the terpenoid backbone biosynthesis pathway have not been globally profiled and the enzyme catalyzing the initiation step from GPP to monoterpenoid biosynthesis has not been identified in P. lactiflora yet. In our datasets, we identified 32 genes with full-length CDSs encoding the enzymes controlling the major catalytic reactions in MVA and MEP pathways (Table S7). Compared with the previous study (Yuan et al., 2013), the genes encoding 10 previously reported enzymes were also identified in our study (Fig. 1). Moreover, we identified the unigene (E_H33980_c1_g4, RLC1) encoding (−)-alpha-terpineol synthase (EC 4.2.3.111) that can catalyze the conversion from GPP to alpha-terpineol (Kulkarni et al., 2013), a monoterpene precursor of paeoniflorin. Our analysis identified the genes encoding the enzymes that almost completely catalyzed the reactions from the glycolysis products to the terpenoid backbone and monoterpenoid biosynthesis in P. lactiflora.

The high accumulation of paeoniflorin in XT and/or BP roots suggests some genes in the paeoniflorin biosynthesis pathway may be highly expressed in XT and/or PB roots. We analyzed the expression patterns of the aforementioned 33 unigenes (Fig. 5). In XT roots, the expression levels of AACT1, HMGS, MK, PMK, GGR, FPS1 were higher than those in GS roots; whereas the expression levels of PMK, MVD, GGR in PB roots were higher than those in GS roots. The other genes did not clearly reveal XT- and/or PB-root differential expression patterns. Note that XT and PB presented the distinct expression pattern of these genes, which was consistent with the fact that the two strains were genetically isolated with a long history (Zhou et al., 2002). To confirm the expression levels determined by the RNA-seq experiment, we used a quantitative Real-time PCR (qRT-PCR) assay to measure the expression levels of several paeoniflorin biosynthesis pathway genes and transcription factor genes (Fig. 6; Fig. S1). The expression levels of the transcription factor gene BHLH94 showed a highly correlation with paeoniflorin levels. The paeoniflorin biosynthesis pathway gene FPS1, AACT1, HMGS, IPI2 and MK as well as the transcription factor gene DREB1D and NAC098 presented root- and/or strain-preferential expression patterns. Our qRT-PCR verification results were highly consistent with the RNA-seq analyses (Fig. 6; Fig. S3).

Accession numbers

The fastq-formatted RNA-seq datasets and the sampling information (accession number CRA001881) are publically available on the Genome Sequence Archive database (BIG Data Center Members, 2018; Wang et al., 2017).