Comparative transcriptome analyses of different Salvia miltiorrhiza varieties during the accumulation of tanshinones

Plant material and RNA sample collection

The selected four varieties of S. miltiorrhiza were collected from Sichuan (SC), Yunnan (YN), Shaanxi (SX), and Henan (HN) (Table 1). All the varieties were asexually propagated through root cutting and cultivated in the experimental field of Chengdu University of Traditional Chinese Medicine for six years. Roots were sampled at 2 days post-anthesis (S1, the early stage of tanshinones accumulation) and 60 days post-anthesis (S2, the late stage of tanshinones accumulation). Each sample had three biological replicates and was frozen in liquid nitrogen as soon as sampled and then stored at −80 °C for subsequent tanshinones contents determination and transcriptome analysis.

Determination of tanshinones contents in S. miltiorrhiza roots

The contents of tanshinones in four varieties of S. miltiorrhiza roots at two periods were determined using ultra-performance liquid chromatography. Tanshinones content was determined on an Agilent 1,290 Infinity UPLC system equipped with a diode array detector (DAD, G4212A) using a ZORBAX Eclipse Plus C18 Rapid Resolution HD HPLC column (2.1 mm × 50 mm, 1.8 um, Agilent). The mobile phase consisted of 0.02% (vol/vol) phosphoric acid solution (A) and acetonitrile (B) at a flow rate of 0.4 mL ∙ min⁻¹. The UPLC chromatogram was monitored at 270 nm and the column temperature was set at 25 °C. The six bioactive compositions including tanshinone IIB (TNIIB), dihydrotanshinone I (DTNI), tanshinone I (TNI), cryptotanshinone (CTN), tanshinone IIA (TNIIA) and miltrinoe (MTN), were compared with the authorized standard (Chengdu Alpha Biotechnology Co., Ltd., China). Total tanshinone in this study was comprised of TNIIB, DTNI, TN1, CTN, TNIIA, and MTN.

RNA extraction, library preparation, and sequencing

Total RNA was extracted using Trizol reagent kit (Invitrogen, Carlsbad, CA, USA) according to the manufacturer’s protocol. RNA’s integrity was confirmed by using the 2,100 Bioanalyzer (Agilent Technologies, Palo Alto, CA, USA) and checked using RNase-free agarose gel electrophoresis. After total RNA was extracted, mRNA was enriched by Oligo(dT) beads. Then the enriched mRNA was fragmented into short fragments using fragmentation buffer and reverse transcripted into cDNA with random primers by using NEBNext Ultra II RNA Library Prep Kit for Illumina (New England Biolabs, USA). Second strand cDNA was synthesized by DNA polymerase I, RNase H, dNTP and buffer. Then the cDNA fragments were purified with QiaQuick PCR purification kit (Qiagen, Germany), end repaired, poly(A) added, and ligated to Illumina sequencing adapters. The ligation products were size selected by agarose gel electrophoresis, PCR amplified, and sequenced using Illumina HiSeq2500 by Gene Denovo Biotechnology Co., Ltd (Guangzhou, China), and 150 bp paired-end reads were generated.

Data processing, assembly, and annotation

Fastp (Chen et al., 2018b) was used to remove sequencing junctions, primers, and low-quality sequences to obtain high-quality clean reads. All downstream analyses were based on high-quality clean data as determined by Q30, meaning that the error rate was less than 0.1%. High-throughput sequencing data were de novo assembled using Trinity (Grabherr et al., 2011).

The unigenes were annotated using BLASTX program against the NCBI non-redundant protein (Nr) database (http://www.ncbi.nlm.nih.gov), the Swiss-Prot protein database (http://www.expasy.ch/sprot), the Kyoto Encyclopedia of Genes and Genomes (KEGG) database (http://www.genome.jp/kegg). When annotating the unigenes with the BLAST program (E-value threshold was set <10⁻⁵), the target sequence with the lowest E-value was selected as the best aligning results. To further assign functions to each unigene, Gene Ontology (GO) annotation was performed by the Blast2GO program according to the categories of Cellular component, Molecular Function, and Biological Process (http://www.geneontology.org/).

Protein sequences obtained from CDS prediction of unigenes were aligned by Hmmscan (Thiriet-Rupert et al., 2016) to search the domain of the plant transcription factors and these unigenes were classified according to plant TFdb3 to predicted transcription factors (Jin et al., 2014).

Gene expression levels calculation and identification of differentially expressed genes

The assembled unigenes were quantified using RSEM (Li & Dewey, 2011) software, and the results were normalized to Fragments Per Kilobase of transcript sequence per Million base pairs sequenced (FPKM). For each unigene, the threshold of FPKM >1 was considered as an expressed gene, and the results could be directly used to compare the differences in gene expression.

DESeq2 (Love, Huber & Anders, 2014) was used to analysis of significant differential expression in two-group comparisons. An adjusted P-value < 0.05 and |log2(FoldChange)| >1 was set as the threshold for a significant differentially expression level. DEGs among four S. miltiorrhiza varieties at two stages were identified by comparing each pair of varieties (i.e., HN-SC, HN-SX, HN-YN, SX-SC, SX-YN, and YN-SC). Tanshinone accumulation-related DEGs were identified by comparing the genes expression of each variety during accumulation stages (i.e., SC1-SC2, YN1-YN2, SX1-SX2, and HN1-HN2).

Statistical analysis

R (version 3.5.1) was used for the basic statistical analysis. To evaluate differences of tanshinones content among the four varieties during the post-anthesis stage, Student’s t-test was carried out. The two-way ANOVA method was used to analyze the effect of different stages, varieties and interaction between stages and varieties on tanshinones accumulation. Significant differences were shown to be statistically significant (P < 0.01). Spearman’s rank correlation coefficient was used to calculate the correlation between gene expression and tanshinone content, and correlation coefficient > 0.70 was set as a significance level.

Validation of RNA-Seq data by quantitative real-time RT-PCR

The validation was performed with the same RNA samples used for RNA-seq analysis. The three biological replicates of the four varieties of S. miltiorrhiza (SC, YN, SX, and YN) at two accumulation stages were analyzed. cDNA synthesis was done by using 400 ng of the total RNA samples with RT Easy™II Kit (Foregene CO., LTD., Chengdu, China). For quantitative PCR, reverse transcribed cDNA products were used as templates. Differentially expressed genes related to tanshinones biosynthesis and candidate transcription factors related to tanshinones accumulation were selected for validation. The primers employed in the qRT-PCR experiments were designed according to the assembled sequences. Primer sequences were reported in Table S1. Actin was used as an endogenous control for the normalization of expression levels of genes (Jiang et al., 2020). The reaction was carried out on a CFX Opus Real-Time PCR System using the Real-Time PCR Easy™-SYBR Green I (Foregene CO., LTD., Chengdu, China) with a total reaction volume of 20 µL, 0.4 μM of the primer, and 80 ng of cDNA. Relative gene expression levels were calculated using the 2^−ΔΔCt method (Yuan et al., 2006). To ensure reproducibility and reliability, three biological replications and three technical replications were implemented for each sample.

Determination of tanshinones contents of different S. miltiorrhiza varieties during tanshinones accumulation

The tanshinones contents of the four S. miltiorrhiza varieties (SC, YN, SX, and HN) were measured at different tanshinones accumulation stages (Fig. 1). The results showed that the accumulation of tanshinones in the four S. miltiorrhiza varieties were different at the two stages under the same environmental conditions (Fig. 1). At S1 stage, the tanshinones contents of each variety were low. The contents of tanshinones in all samples was significantly different between S1 and S2 stage (P < 0.01, Table S2), and the content of total tanshinone at S2 was about 6.48 ~ 20.77 times higher than that at S1. By comparing the tanshinones content in different varieties of S. miltiorrhiza, it was observed that at S2 the content of CTN, TNIIA, and TNI increased by 7.71, 7.88, and 24.01 times respectively. While, compared with at S1, the content of DTNI, TNIIB and MTN increased only by 3.44, 5.00 and 5.42 times at S2. This phenomenon is consistent with the research previous results of Yang et al. (2013), Chang et al. (2019), Contreras et al. (2019), and Zhan et al. (2019) in differential tanshinones contents accumulating in varied degrees. The two-way ANOVA of tanshinone showed that both stage and variety had significant effects (P < 0.01) on tanshinones accumulation, and the interaction between them was significant (P < 0.01), indicating that the combination of varieties and growth stages had an impact on the tanshinones accumulation (Table S3).

Sequencing, assembly, and functional annotation of genes in four varieties of S. miltiorrhiza

To explore the expression of genes in the four varieties of S. miltiorrhiza during tanshinones accumulation, the transcriptome at tanshinones accumulation stages of the four varieties was analyzed by RNA-Seq. A total of 24 cDNA libraries were constructed from root samples with three biological replicates for each stage and each S. miltiorrhiza variety after removing the adapter reads. The average GC content of the transcriptome was 41.17%. Each sample yielded more than 5.27 Gb clean data (Table S4). The Q20 and Q30 percentages were more than 92.80% and 75.22%, respectively. The transcriptome sequencing data from the four S. miltiorrhiza varieties were filtered and assembled by Trinity resulting in a total of 70,357 unigenes. Overall, the unigenes obtained showed a total assembled bases of 75,359,646, a mean length of 1,071 base pairs (bp) and a mean N50 of 1,911 bp.

These assembled unigenes were searched against the NCBI Nr, SwissProt, KEGG, and KOG databases, and it was found that approximately 56.88% of the unigenes could be annotated in at least one database. The Venn map revealed that 19,454 unigenes were found to be common annotated among four databases (Fig. S1). By sifting out the genes with low expression levels (FPKM < 1), 52,216 genes were obtained (Table S5). KEGG database can link genomic information with functional information, and with KEGG annotation, the functional genes needed for different research can be deeply mined. In this study, a total of 112 genes related to tanshinones biosynthesis were identified, including the terpenoid backbone biosynthesis pathway (Ko00900), diterpenoid biosynthesis pathway (Ko00904) (Table S6).

Differential gene expression regulation in four S. miltiorrhiza varieties

To understand the different regulation of gene expression in different varieties of S. miltiorrhiza during tanshinones accumulation stages, we identified the DEGs between each pair of varieties (i.e., HN-SC, HN-SX, HN-YN, SX-SC, SX-YN, and YN-SC) at S1 and S2 (Fig. 2). Among the pairwise comparison groups, a total of 11,368 and 17,639 DEGs were detected in S1 and S2 stages respectively (Fig. S3). Only one IREG3 (Unigene 0009899) gene was differentially expressed in all the pairwise comparison groups at S1 (Fig. S3A), while there was no common differentially expressed gene at S2 (Fig. S3B). Additionally, the number of DEGs in each pairwise comparison during tanshinones accumulation revealed that the dynamic regulation of gene expression was varied among the four varieties (Table 2). Except for YN-HN and YN-SC, the number of DEGs at S1 was more than that at S2 in the other four pairwise comparison groups (Table 2).

KEGG enrichment analysis is helpful to further understand the potential function of the genes. In stage S1, DEGs in six pairwise comparison groups were mapped to 100–128 KEGG pathways (Table S7). Of these KEGG pathways, 88 biological processes related to primary metabolism, secondary metabolism, and signal transduction processes were identified in all pairwise comparison groups (HN1-SC1, HN1-SX1, HN1-YN1, SX1-SC1, SX1-YN1 and YN1-SC1) (Table S7). In stage S2, DEGs in six pairwise comparison groups were mapped to 73–136 KEGG pathways, and 51 pathways related to secondary metabolism and signal transduction were identified in all pairwise comparison groups (Table S8). Of these KEGG pathways, 48 metabolic processes including plant-pathogen interaction, MAPK signaling pathway, galactose metabolism, pyruvate metabolism, amino sugar and nucleotide sugar metabolism, and oxidative phosphorylation were the identified in all pairwise comparison at the two stages (Table S7). At stage S1, DEGs were specifically enriched in some processes of primary metabolism (photosynthesis and citrate cycle) and terpenoid backbone biosynthesis. At stage S2, DEGs were mainly related to metabolic processes (thiamine metabolism and propanoate metabolism) and plant circadian regulation (Table S8).

DEGs related to the changes of tanshinones content during the accumulation of tanshinones

To comprehensively understand the regulation of S. miltiorrhiza gene expression during the accumulation of tanshinones period, the gene expression of the four varieties in two stages (i.e., SC1-SC2, YN1-YN2, SX1-SX2, and HN1-HN2) were compared. A total of 2016 DEGs containing 18 genes involved in the diterpenoid biosynthesis were obtained (Table S9). Among them, 10 DEGs including HSP TF (Unigene0006797), RBG (Unigene0014447), and GAP3 (Unigene0014163) etc. were differentially expressed in all the four pairwise comparison groups (Fig. 3A, Table S9). Compared with stage S1, the expression of these TFs at stage S2 was up-regulated (Figs. 3B, 3C). The KEGG enrichment analysis of 2016 DEGs showed that photosynthesis, secondary metabolites biosynthesis, oxidative phosphorylation, diterpenoid biosynthesis, terpenoid backbone biosynthesis, sesquiterpene and triterpene biosynthesis were the most differentially regulated biological processes during the tanshinones accumulation stages among four varieties of S. miltiorrhiza (Fig. 3D).

To identify the tanshinones accumulation-related genes, the correlation between 2016 DEGs with the significant accumulated components (DTNI, TNIIA, and TNI) was analyzed, respectively. Using this method, 277 DEGs correlated with all the three components were obtained, containing one dehydrogenase gene and three CYP450 genes (Table S10, Fig. S4). To understand which categories of these genes were overrepresented, these DEGs were further analyzed by a KEGG enrichment analysis. The KEGG pathway analysis mapped to 63 categories, including photosynthesis, carbon fixation in photosynthetic organisms, galactose metabolism, oxidative phosphorylation, biosynthesis of secondary metabolites, ubiquinone and other terpenoid-quinone biosynthesis, metabolic pathways, MAPK signaling pathway, and plant-pathogen interaction (Table S10). The top 20 KEGG pathways with the highest representation are shown in Fig. 3E. These results suggest that the accumulation of tanshinones is not only related to the process of tanshinones biosynthesis, but also related to other secondary metabolic processes, signal transduction processes, and even primary metabolic processes.

Transcription factors (TFs) are a class of proteins that specifically bind to gene promoters and regulate their expression at different levels. In our transcriptome analysis, a total of 1,420 TFs was identified, which could be classified to 55 families (Table S11). During the tanshinones accumulation stage, a total of 124 differentially expressed TFs were identified in the four varieties (i.e., SC1-SC2; YN1-YN2; SX1-SX2; HN1-HN2) (Table S12). Among the genes associated with tanshinone accumulation, there are 24 TFs, including 6 NAC, 4 AP2/ERF, 3 WRKY, 2 MYB and 1 MYB-related transcription factor (Table S12). These TFs were up-regulated in all varieties of S. miltiorrhiza at stage S2, and significantly up-regulated in YN2 with high tanshinones concentrations (Fig. 3F).

Experimental validation of differential expressed genes by qRT-PCR

To verify the reliability of the RNA-seq results, 10 DEGs at S1 and 10 DEGs at S2 were selected for qRT-PCR analysis (Table S1). These gens included differentially expressed genes related to tanshinones biosynthesis and candidate transcription factors related to tanshinones accumulation (Table S1). The expression levels of these genes were consistent with those determined using RNA-Seq both in stage S1 (R² = 0.8066) and S2 (R² = 0.8484) (Figs. 4, S5), indicating that the dataset obtained by RNA-seq was reliable for gene expression analysis.