Transcriptomic analyses provide new insights into green and purple color pigmentation in Rheum tanguticum medicinal plants

Plant materials, sample collection, and preparation of total RNA

Fresh tissues from both purple and green R. tanguticum (Voucher Number 51322420120805441) were collected from the Aba Tibetan and Qiang autonomous prefecture (33°68′N, 103°25′E) located in the Sichuan province of China by Yangfan Tang (Sichuan Academy of Chinese Medicine Sciences) and taxonomically identified by Qingmao Fang (Sichuan Academy of Chinese Medicine Sciences). The collection of plant material in this study complied with all relevant institutional, national, and international guidelines and legislation. The samples were flash-frozen in liquid nitrogen on-site. A total of 19 R. tanguticum samples were collected from different parts of R. tanguticum plants: five green leaf biological duplications, five green petiole biological duplications, two green rhizome biological duplications, four purple leaf biological duplications, and three purple petiole biological duplications. Total RNA was extracted with the CTAB-pBIOZOL (CAT#BSC55M1), according to the manufacturer’s instructions.

The statistical power of this experimental design, calculated in RnaSeqSampleSize (https://github.com/slzhao/RnaSeqSampleSize), was 0.7.

Quality evaluation of total RNAs, library preparation, and sequencing

The evaluations of the RNA samples were carried out using Qubit 2.0, Nanodrop, and Agilent 2100 (Simbolo et al., 2013) to ensure that the concentration, integrity, and purity were suitable for library preparations and RNA sequencing. The samples with an RNA integrity number (RIN) value over seven were moved to library preparation using the MGIEasy RNA kit (CAT# 1000006383). The constructed library was used for sequencing by the BGISEQ-500 (Huang et al., 2018) platform.

Pre-processing of raw reads and de novo transcriptome assembly

We used FastQC (version 0.11.3) (Andrews, 2010) to confirm the validation of the raw data and the low quality reads were filtered out using Trimmomatic (version 3) (Bolger, Lohse & Usadel, 2014). We then performed a second quality validation of just the clean reads, using FastQC, to ensure they were suitable for downstream analyses. All FastQC results were visualized using MultiQC (version: 1.9) (Ewels et al., 2016). Clean reads were mapped to the unigene using Bowtie2 (version 2.2.5) (Langmead & Salzberg, 2012). The Trinity (Haas et al., 2013) (v2.9.1) software was used to assemble the short k-mers into contigs based on the de Bruijin Graph algorithm. The sequencing depth of all used data was 678 ×.

Transcriptome annotation

We used the TransDecoder (v3.0.1; http://transdecoder.sourceforge.net) to predict the potential coding sequence of the unigenes. First, the longest open reading frame (ORF) was selected from the transcript sequences and then scanned with known protein sequences (Swiss-Prot 2020 database and Pfam-A.hmm 2020 database) using blastp and hmmscan. These results were then combined to predict the coding region of R. tanguticum.

The assembled unigene sequences were functionally annotated by aligning them with the Kyoto Encyclopedia of Genes and Genomes (KEGG) (Kanehisa & Goto, 2000), Gene Ontology (GO) (Ashburner et al., 2000; The Gene Ontology, 2021), clusters of orthologous groups for eukaryotic complete genomes (KOG) (Tatusov et al., 2000), SwissProt (Bairoch & Apweiler, 2000), Pfam (http://pfam.xfam.org/), and National Center for Biotechnology Information (NCBI) non- redundant (NR) protein databases, as well as the Nucleotide Sequence Database (NT) with BLAST ( E-value ≤ 1e–05).

Differential gene expression analysis

The RSEM pipeline was used to determine the FPKM (fragments per kilobase per million) value of different samples. The clean data were re-mapped to the assembled transcriptome using bowtie2 (version 2.2.5) (Langmead & Salzberg, 2012). Bowtie-build was used to make a reference library. The “–transcript-to-gene-map” parameter was then used to map the transcripts to the corresponding genes. Filtered sequencing reads were mapped to the reference by bowtie2 and then the “rsem-calculate-expression” option was used to quantify the expression level. A differential expression analysis was carried out with the DESeq2 (version 1.28.1) (Love, Huber & Anders, 2014) package in R (version 4.0.2). The screening criteria for differential expression genes were: adjusted pvalue <0.05, log2FoldChange >1 (up-regulated) or log2FoldChange <−1 (down-regulated). The information of the grouping and control samples from the differential expression gene analysis is summarized in Table S1. Based on the DEseq2 normalized data matrix, we calculated the Pearson correlation coefficients between different samples. The co-expression analysis of anthraquinone pathway genes and TFs was performed using the WGCNA R package.

Ethics approval

This study, including sample collection, was conducted according to the ethical clearance of the 10,000 Plant Genomes Project—10KP (NO. FT20026) by the institutional review board of BGI which permits the use of biological resources for scientific research purposes.

RNA data quality assessment and de novo assembly

In R. tanguticum, we constructed short-read RNA-seq libraries and a total of 221.09 Gb raw data and 213.94 Gb high-quality data were generated for quantification, and differential gene expression analyses. The percentage of clean reads ranged from 95% to 99%, the mapping rates were about 70% in R. tanguticum, and the Q20 value of the clean reads was around 90%. The clean data from a total of 19 samples were subjected to further transcriptome assembly. We retrieved 336,987 transcripts and 120,261 unigenes with contig N50 1,761 bp and contig N50 1,474 bp, respectively. The completeness of the transcriptome was evaluated using the Benchmarking Universal Single-Copy Ortholog (BUSCOs) resulting in 91.9% (S:25.5%, D:66.4%, F:5.7%, M:2.4%) coverage. The reads number, clean reads percentage, mapping rate, Q20 value, and GC content of the data were assessed and the results are summarized in Table S2.

The PCA analysis showed that the samples in the dataset clustered into four major groups corresponding to the tissue types and plant colors (Fig. 1A). PC1 and PC2 explained 44% and 33% of the total variance in gene expression of the R. tanguticum data set, respectively. The heatmap clustering of Pearson correlation coefficients from the comparison of all 19 samples with various tissue types revealed a strong correlation (0.81–1) between replicates (Fig. 1B). Taken together, these results suggest that our datasets are a reliable data resource for future studies.

Functional annotation

A total of 45,866 out of the 120,261 total unigenes (38.14%) showed significant similarities to known proteins in R. tanguticum (Table S3). Gene Ontology (GO) assignments were used to classify the unigene sequences based on functional annotation to determine their possible functions in R. tanguticum. There were 27,423 unigenes that could be categorized into functional groups under the “cellular component,” “molecular function,” and “biological process,” divisions in R. tanguticum (Fig. 2A). For the biological process group, genes involved in the “cellular process” (12,518) and “metabolic process” (10,908) were the most highly represented. For unigenes in the cellular component group, “cellular anatomical entity” (14,987) and “intracellular” (8,183) were the most highly represented categories, followed by “protein-containing complex” (2,519). For the molecular function group, “binding” (15,434), “catalytic activity” (12,960), and “transporter activity” (1,535) were the most represented categories.

A functional characterization analysis based on the KEGG database was performed on the unigenes generated in the present study. In summary, the most represented pathways in R. tanguticum were: “global and overview maps” (9,223), “carbohydrate metabolism” (4,170), “translation” (2,960), and “folding, sorting and degradation” (2,172) (Fig. 2B). The identified transcripts enriched in these diverse metabolic pathways will help us better understand the active ingredients in R. tanguticum.

Transcriptome expression and differentially expressed genes of R. tanguticum

The FPKM values of all unigenes in different samples were summarized in Table S4. To investigate the different colors between the leaves and petioles of R. tanguticum, we calculated the expression levels of the unigenes between these two tissues. Comparing green leaves (GL) and purple leaves (PL) showed 4,861 DEGs in total, with 2,426 up- and 2,435 down-regulated unigenes, respectively. Further, we did a KEGG enrichment analysis for these DEGs. Among the down-regulated genes, 51 genes were assigned to “photosynthesis,” 48 genes clustered in “carbon fixation in photosynthetic organisms,” and 22 genes were assigned to “photosynthesis –antenna proteins” (Fig. 3A). Interestingly, the unigenes enriched in secondary metabolism, such as the “flavone and flavonol biosynthesis” (eight) and “sesquiterpenoid biosynthesis” (19) pathways, were up-regulated (Figs. 3A; 3C). We also compared the green petioles and purple petioles and found 3,028 up- and 2,413 down-regulated unigenes. Notably, there were no up- or down-regulated unigenes associated with photosynthesis-related pathways. Three major up-regulated genes were enriched in “metabolic pathways” (694), “biosynthesis of secondary metabolism” (346), and “amino sugar and nucleotide sugar metabolism” (120). These results suggest that the expression profiles of the petiole are different than the expression profiles of the leaves (Figs. 3B; 3D).

We also performed a GO enrichment analysis. In this analysis, many down-regulated unigenes between the green and purple leaves were also found to be enriched in photosynthesis-related pathways, such as “photosystem II stabilization” (five), “photosystem II assembly” (six), “photosystem I reaction center” (10), and “photosynthesis, light harvesting” (23) (Fig. 4; Table S5). The GO enrichment analysis of up-regulated unigenes between GL and PL found that they were enriched in cell cycle-associated genes including “microtubule-based process” and “nucleic acids metabolisms process” (Table S6). The up-regulated genes in the petioles (GP vs. PP) were enriched in cell wall metabolism and development-associated genes (Table S7), while the down-regulated genes in the petioles were enriched in amino acids transport and cell wall catabolism (Table S8).

To compare tissue-specific expression patterns, we compared the expression profiles of the leaves and petioles of the same colors. The comparative transcriptomic study between GL and GP found 2,004 up-regulated and 1,690 down-regulated DEGs. The enrichment analysis results of these DEGs are summarized in Tables S9 and S10. By comparing the PL and PP, we found 2,687 up-regulated and 1,745 down-regulated DEGs. The results of the KEGG and GO enrichment analyses of the down- and up-regulated DEGs between PL and PP are summarized in Tables S11 and S12.

Comparing the expression of unigenes involved in chlorophyll and anthocyanin in R. tanguticum

Chlorophyll is the most abundant pigment on earth and is a key component of photosynthesis required for the absorption of sunlight (Hörtensteiner & Kräutler, 2011). We identified chlorophyll pathway genes (Table S13) and then constructed the expression maps of those genes. Interestingly, some enzyme genes (RtGSA, RtHEM, RtHEMC, RtHEME, RtHEMF, RtHEMG, RtCAO, RtNYC1/RtNOL, RtCLH, RtSGR) were highly expressed in green samples. In contrast, we found that genes involved in chlorophyll degradation like RtPPH (pheophytinase), RtPaO (pheophorbide a oxygenase) and RtRCCR (red chlorophyll catabolite reductase) exhibited higher expression levels in purple samples (Fig. 5). The purple color of R. tanguticum could be due to the combination of lower chlorophyll biosynthesis expression and higher RtPPH, RtPaO, and RtRCCR expression in the chlorophyll degradation pathway.

Anthocyanin is a major group of plant pigments that may appear red, purple, blue, or black in various tissues. The analysis of the R. tanguticum transcriptome data set revealed that 73 unigenes exerted a direct influence over 10 enzymes that are known to be involved in the anthocyanin pathway, and almost all of those 73 unigenes were demonstrated to be a multigene family. Anthocyanidin synthase (ANS) and Anthocyanidin 3-O-glucoside 2”’-O-xylosyltransferase (UFGT) catalyze the last two steps of anthocyanidin biosynthesis and are therefore key enzymes in the biosynthesis of blue or red pigments (Chen et al., 2011; Zhang et al., 2020). However, these two genes were more highly expressed in green samples compared to the purple tissues, which did not support the phenotypes we observed.

Expression profile of unigenes involved in the anthraquinone biosynthesis pathway

Anthraquinones are bioactive compounds in Rheum plants. Their biosynthetic pathway is thought to involve the shikimate, MVA, MEP and polyketide pathways (Leistner, 1971). To identify the potential candidate unigenes of the anthraquinones biosynthetic pathway, we identified homologous genes by aligning the transcriptome sequence of R. tanguticum with all the known enzyme sequences associated with the above-mentioned pathways. We then screened these results by combining the sequence similarity and functional annotation to obtain the most conceivable candidate unigenes involved in the anthraquinones biosynthesis pathway. Using this process, we identified 79 unigenes associated with the anthraquinones biosynthesis pathway in R. tanguticum (Table S14). In the shikimate pathway, most RtDAHPS, RtDHQS, RtMenE, and RtEPSP genes were highly expressed in green leaves. Interestingly, in the MEP pathway, most RtISPG genes were highly expressed in purple samples. Other catalytic genes (RtDXPS, RtDXR, RtISPD, RtISPF and RtISPH) had varying expression levels in the green and purple plants. In the MVA pathway, RtHMGR genes showed tissue-specific expression regardless of color types in both green and purple petioles. The R. tanguticum rhizome is considered to be the highest quality of all medical rhubarb in medicinal material markets. We found that RtHMGS, RtHMGR, RtMK, RtPMK and RtMPD were all highly expressed in green rhizome samples. In the PKS pathway, RtPKS was highly expressed in green rhizomes and RtPKC expressions were higher in green leaves.

The regulatory network of the anthraquinone biosynthetic pathway in R. tanguticum showed differential expression patterns (Fig. 6A). The co-expression analysis of anthraquinone and TFs showed that RtSMK, RtMK, RtDXPS, RtHMGR, RtHMGS, RtMenB, RtMPD, and RtSDH gene family members have expression patterns that closely resemble multiple TFs, including bHLH, WRKY, MADS, and MYB TFs (Fig. 6B). We cannot rule out the possibility that anthraquinone is synthesized in the leaves and petioles of R. tanguticum and then transported below ground for long-term storage. In future studies, the transcriptional network involved in the regulation of anthraquinone biosynthesis and transportation should be investigated.