Integrated transcriptome and metabolome analyses revealed regulatory mechanisms of flavonoid biosynthesis in Radix Ardisia

Determination of total saponin, flavonoid, and phenolic content

Ardisia Crenata Sims (ACS) and Ardisia Crispa (Thunb.) A.DC. (AC) plants were collected in Guiyang County, Guizhou Province, China (N.109.437569, E.19.19680). The samples were authenticated by Professor Sheng Hua Wei and stored at the Laboratory of Traditional Chinese Medicine and Ethnic Medicine, Guizhou University of Traditional Chinese Medicine, Guiyang, Guizhou, China. The specimen accession number was GUTCM0059.

Total saponin content (TSC) was measured according to the method described by Medina-Meza et al. (2016) with minor modifications. 1 mL of ACS or AC extract was added to a 10 mL test tube, and then 0.2 mL 5% Vanillin acetic acid solution and 0.8 mL HCLO₄ were added after 6 min, vortexed vigorously, and heated at 60 °C for 10 min; after the mixture was cooled in ice water, five mL glacial acetic acid was added. The absorbance was measured at 545 nm. Glacial acetic acid was used as the blank. The TSC was obtained by comparison with the standard curve of oleanolic acid; the analysis was done in triplicate.

Total flavonoid content (TFC) was determined according to the method described by Xu et al. (2014) with some modifications. 1 mL Ardisia Crenata Sims (ACS) or Ardisia Crispa (Thunb.) DC (AC) extract was added to a 10 mL test tube, and then 2.4 mL 70% ethyl alcohol and 5% NaNO2 (0.4 mL) were added after 6 min. Then, 0.4 mL 10% Al (NO₃)₃ was added. After 6 min, 4 mL 4% NaOH and 70% ethyl alcohol were added at room temperature for 15 min. It was then measured against ethyl alcohol as the blank at 510 nm, a calibration curve of rutin was plotted to calculate the TFC, and the analysis was done in triplicate.

The analysis of total phenolic content (TPC) was carried out according to the Folin-Ciocalteau spectrophotometric method, with some modifications (Vieira de Morais et al., 2021). The 1 mL gallic acid (standard solution) or ACS and AC extract and 3 mL Folin–Ciocalteau solution (10% in water) were pipetted into a 25 mL test tube. After 5 min, 6 mL 10% sodium carbonate aqueous solution was added to each tube. A control was prepared by replacing the sample with distilled water. The absorbance was measured at 750 nm after 50 min. The TPC was calculated by linear regression using gallic acid as the standard; all samples were analyzed in triplicate.

RNA extraction, library preparation and sequencing

In this study, 12 root samples produced from Ardisia Crenata Sims (ACS) and Ardisia Crispa (Thunb.) DC (AC) (six biological replicates) were tested using RNA-seq. Total RNA was extracted with the RNAprep Pure Plant Kit (Tiangen, China) and detected on 1% agarose gel. The total RNA concentration was determined using an Agilent 2100 Bioanalyser, and the purity of the samples was measured using a NanoDrop™ ultraviolet spectrophotometer. The purity, concentration, and integrity of the total RNA samples were assessed prior to further analysis. After the isolation and fragmentation of the total RNA, mRNA was enriched using poly-T oligo-attached magnetic beads. The enriched mRNA was fragmented into short fragments using M-MuLV Reverse Transcriptase and reverse transcribed into cDNA using random hexamer primer. Second strand cDNA was synthesized using DNA Polymerase I and RNase H. Then, 3 µl USER Enzyme (NEB, USA) was used with size-selected, adaptor-ligated cDNA at 37 °C for 15 min followed by 5 min at 95 °C before PCR. Finally, PCR products were purified (AMPure XP system) and library quality was assessed on the Agilent Bioanalyzer 2100 system. After cluster generation, the library preparations were sequenced using an Illumina HiSeqTM 4000 platform by Biomarker Technologies (Beijing, China). In this work, twelve cDNA libraries produced from ACS and AC were sequenced using the Illumina HiSeq™ 4000 platform and 150bp paired-end reads were generated. Clean reads were obtained by removing lower quality reads as well as raw reads containing an adapter or ploy-N. Transcriptome assemblies for the twelve libraries were performed separately using Trinity (Garber et al., 2011) with min_kmer_cov set to two by default and all other parameters set to default.

DEGs related to secondary metabolism pathways

Gene expression levels were calculated using RSEM (Li & Dewey, 2011) for each sample. The differentially expressed genes (DEGs) were identified by comparing raw readcounts using the DESeq 2 (Love, Huber & Anders, 2014). Genes with FC ≥ 1.0 and p-value <0.05 were identified as DEGs. We used the OmicShare tools (https://www.omicshare.com/tools) and KOBAS 2.0 (Xie et al., 2011) software to test the statistical enrichment of DEGs in KEGG pathways, and the significance of KEGG terms was assessed using the Bonferroni corrected Fisher exact test (P < 0.05). We then selected the pathways associated with flavonoid biosynthesis for a more detailed analysis.

qRT-PCR verification

In order to verify the accuracy of the RNA-seq, a qRT-PCR analysis was used to assess the quality of the RNA-seq data. Total RNA was extracted from root tissue with TRIzol according to the manufacturer’s instructions (TIANGEN, China), and the isolated RNA was reverse-transcribed into cDNA with the RevertAid First Strand cDNA Synthesis Kit (Thermo Fisher, China). The expression of randomly selected genes was monitored by qRT-PCR using the SYBR Green qPCR Mastermix (TIANGEN, China) real-time PCR system by Bio-Rad CFX 96™ following the manufacturer’s instructions. Detailed information about the primer sequences for qRT-PCR was provided by the primer3 platform (http://frodo.wi.mit.edu/primer3/; Table S6). The qRT-PCR was used on three biological replicates with the expression levels of the genes determined using the 2^−−ΔΔCT method. Expression levels were normalized against the GAPDH. Data are represented as mean values  ± standard deviation, and the GraphPad Prism 6 software was employed to draw the histogram.

Metabolite profiling using UHPLC-QTOF-MS

In this study, 12 samples from ACS and AC (six biological replicates, respectively) were scanned for metabolite determination by UHPLC-QTOF-MS, and the samples were consistent with RNA-seq.25 mg of root tissue sample after liquid nitrogen grinding was weighed and placed in a 1.5 ml EP tube, and 500 µL extract solution (acetonitrile: methanol: water = 2: 2: 1) containing isotopically-labelled internal standard mixture was added to each EP tube. After a 30s vortex, the samples were placed in the TissueLyser at 35 Hz for 4 min and sonicated for 5 min in an ice-water bath. The samples were then incubated for 1 h at −40 °C and centrifuged at 12,000 rpm for 15 min at 4 °C. The supernatant was transferred to a fresh EP tube and dried at room temperature. Then, the samples were reconstituted in 300 µL 50% acetonitrile by sonication for 10 min, then centrifuged at 13,000 rpm for 15 min at 4 °C, and 75 µL of supernatant was transferred to a fresh glass vial for LC-MS analysis.

Metabolomics data analysis

The UHPLC-MS/MS analyses were applied to the 1290 Infinity series UHPLC System (Agilent Technologies) coupled with a TripleTOF 6600 mass spectrometer (AB Sciex) from Biomarker Technologies (Beijing, China). The MS raw data were converted to the mzXML format by ProteoWizard, and processed by R package XCMS (version 3.2). This process included peak deconvolution, alignment and integration. Minfrac and cut off were set as 0.5 and 0.3, respectively. An in-house MS/MS database was applied for metabolites identification. The principal component analysis (PCA) and partial least squares discriminant analysis (PLS-DA) were performed using R. We employed a univariate analysis to calculate statistical significance; metabolites with VIP ≥ 1, fold change ≥ 2.0 or ≤ 0.50, and a p-value <0.05, were considered differentially expressed metabolites (DEM). Volcano plots were used to filter the metabolites of interest based on the log2 (FC) and log10 (p-value) of the metabolites. The KEGG database was used to annotate the differentially expressed metabolites with a p-value less than 0.05 set as the threshold. The enrichment factor represented the ratio between the proportion of differentially expressed metabolites in the pathway and all the metabolites in the pathway; the greater the value, the greater the degree of enrichment.

Combined transcriptome and metabolome analyses

To uncover the regulatory mechanism of flavonoid biosynthesis, a correlation analysis was performed using Pearson’s correlation coefficient to calculate the correlation coefficient (CC). The p-value of DEG and DEM, and the —CC—>0.80 and CCP <0.05 were used for the cluster analysis. The correlation matrix was imported into the Cytoscape software to visualize the DEG and DEM network. The canonical correlation analysis (CCA) is statistical technique for studying associations between two sets of variables (Hotelling, 1936), and applied to the integration of data originating from different “omics” technologies (Le Cao, Gonzlez & Djean, 2009).

Total saponin, flavonoid, and phenolic content in Radix Ardisia from different origins

The total saponin content (TSC), total flavonoid content (TFC), and total phenol content (TPC) in Radix Ardisia from different origins were measured using the colorimetric method. The TSC in ACS was markedly higher than in ACSV (P < 0.05), however, there was no difference between the TSC in ACS and AC. The TFC and TPC in ACS were dramatically higher than in AC (P < 0.05) (Fig. 1). The differences in gene expression at varying origins in Radix Ardisia may lead to differences in the accumulation of secondary metabolites. We performed RNA-seq to examine the differences in gene expression between ACS and AC, and explore the regulation of flavonoid biosynthesis.

RNA sequencing and assembly results

We compared transcriptional profiles, using RNA-Seq. The root samples from ACS and AC were each sequenced in six replicates. Twelve RNA-Seq libraries of root tissue generated approximately 89.73 GB of clean data. The high-quality reads totaled 52,249 Unigenes, with an average length of 1440 bp; N50 was 2336 bp, the alignment rate was more than 67.52%, and the Q30 base percentages were greater than or equal to 94.85% (Table 1). The results indicated that the assembly quality was reliable, and satisfactory for further analysis. The raw sequencing data were deposited in the NCBI Bioproject database under accession number PRJNA739135. The Unigenes data were deposited in the Transcriptome Shotgun Assembly(TSA) Database under accession number GJZC00000000, and SUBID was SUB11534240.

Identification of differentially expressed genes

A total of 52,249 Unigenes were identified, which were expressed in at least one sample (Table S1). DESeq 2 was used to analyze the differentially expressed genes (DEGs) in Ardisia Crispa (Thunb.) DC. from origins ACS and AC. There were 17,685 DEGs between ACS vs. AC: 8,854 were upregulated and 8,831 were downregulated (Table S2; Fig. 2). To confirm the reliability of the gene expression, 20 genes were randomly selected for quantitative real-time polymerase chain reaction (qRT-PCR) detection. The qRT-PCR showed that the tendency of gene expression was similar to the RNA-Seq results (Figs. 3A and 3B). The results showed that the RNA-Seq results were reliable in this study.

DEG functional enrichment analysis

A KEGG pathway analysis showed that 139 pathways were represented in the transcriptome dataset (Table S3). Interestingly, these DEGs in ACS and AC were mainly enriched in the metabolism, genetic information processing, environment information processing, cellular process, and organismal systems, involved in Inositol phosphate metabolism, Other glycan degradation, Homologous recombination, Mismatch repair, Phosphatidylinositol signaling system, Glycosylphosphatidylinositol (GPI)-anchor biosynthesis, Autophagy - other, Ubiquitin mediated proteolysis, Basal transcription factors, Galactose metabolism, mRNA surveillance pathway, Glycosaminoglycan degradation, Butanoate metabolism, Fructose and mannose metabolism, Phenylalanine, tyrosine and tryptophan biosynthesis, Circadian rhythm - plant, Starch and sucrose metabolism, Nicotinate and nicotinamide metabolism, and Taurine and hypotaurine metabolism (P < 0.05) (Fig. 4; Table 2).

DEGs related to flavonoid biosynthesis

DEGs involved in the flavonoid biosynthetic pathway included those integral to phenylpropanoid biosynthesis (ko00940, 52 genes), phenylalanine metabolism (ko00360, 26 genes), flavonol and flavone biosynthesis (ko00944, three gene), flavonoid biosynthesis (ko00941, 32 genes) and anthocyanin biosynthesis (ko00942, one gene). Table S4 showed that 100 DEGs (49 upregulated and 61 downregulated) were differentially expressed in ACS and AC. Furthermore, eight chalcone synthase (CHS, Unigene_171144, Unigene_186235, Unigene_108709, Unigene_182412, Unigene_168887, Unigene_030909, Unigene_001289, Unigene_029241), one chalcone isomerase (CHI, Unigene_026449), six dihydroflavonol 4-reductase (DFR, Unigene_030972, Unigene_181720, Unigene_106630, Unigene_173665, Unigene_092658, Unigene_180983), two flavonoid 3′5′-hydroxylase (F3′5′H, Unigene_094347, Unigene_094346), one flavanone 3-hydroxylase (F3H, Unigene_026466), two flavonol synthase (FLS, Unigene_013848, Unigene_015774), one leucoanthocyanidin reductase (LAR, Unigene_174994), one phenylalanine ammonia-lyase (PAL, Unigene_177796), four 4-coumarate: CoA ligase (4CL, 008420, Unigene_171229, Unigene_025326, Unigene_007828), and two cinnamate 4-hydroxylase (C4H, Unigene_181127, Unigene_013550, not significant) were identified in the second developmental stage (Fig. 5).

Metabolic differences in ACS and AC

RNA-Seq results demonstrated significant differences in metabolism between ACS and AC. Therefore, we investigated the changes in the metabolic compositions of ACS and AC. In this study, we used 12 samples to survey the differences in the metabolic constituents of ACS and AC with six biological replicates. First, the principal component (PC1) in ESI+ mode (34% of the total variables) and PC1 in ESI − mode (47.20%) were clearly separated between the ACS and AC groups (Figs. S1A and S1B). The VIP values of the first two principal components of the multivariate orthogonal partial least square discriminant (OPLS-DA) (Figs. S1C and S1D), and VIP ≥ 1, fold change ≥2.0 or ≤ 0. 50, and p-value <0.05 to screen for differentially expressed metabolites (DEM). A total of 943 and 1,250 DEMs were identified between the ACS and AC groups in the ESI+ (ESI −) mode, respectively (Table 3; Fig. 6).

Metabolic pathway analysis

A DEM analysis and KEGG pathway annotation of metabolites revealed that 124 DEMs were enriched in a variety of functional pathways (Table S5), including: flavonoid biosynthesis, phenylalanine metabolism, zeatin biosynthesis, indole alkaloid biosynthesis, isoflavonoid biosynthesis, ateroid biosynthesis, flavone and flavonol biosynthesis, fatty acid biosynthesis, phenylpropanoid biosynthesis, phenylalanine metabolism, tyrosine metabolism, pyrimidine metabolism, glutathione metabolism, and glucosinolate biosynthesis (Figs. 7A and 7B). In this study, 11 metabolites were authenticated including L-Phenylalanine, p-Hydroxycinnamaldehyde, L-Tyrosine, Ferulic acid, Sinapic acid, Sinapyl alcohol, 4-Hydroxycinnamic acid, Naringenin, Luteolin, Catechin, and Quercetin, all of which participated in flavonoid biosynthesis and the upstream phenylpropane metabolic pathway (Fig. 7C). Compared with ACS, the expression levels of six metabolites (L-Phenylalanine, p-Hydroxycinnamaldehyde, L-Tyrosine, Ferulic acid, Sinapic acid, Sinapyl alcohol) were significantly increased in AC, while the expression levels of five metabolites (4-Hydroxycinnamic acid, Naringenin, Luteolin, Catechin, Quercetin) were significantly decreased in AC.

Correlation analysis between RNA-seq and metabolites uncovers the regulatory pathway of flavonoid biosynthesis

To explore the correlation between gene expression and metabolites, we performed correlation analyses of the metabolites related to flavonoid biosynthesis and the transcripts. The profiles of the metabolites and gene expression in ACS and AC were compared using a correlation coefficient and a canonical correlation analysis (CCA). These correlations and the CCA analysis results are shown in Figs. S2 and S3. The results indicated that metabolites such as Succinate, Salicylic acid, 2-Hydroxyphenylacetic acid, p-Hydroxycinnamaldehyde, L-Tyrosine, Phenylpyruvate, N-Acetyl-L-phenylalanine, Phenylacetic acid, Sinapic acid, Sinapyl alcohol, 4-Hydroxycinnamic acid, Alpha-N-Phenylacetyl-L-glutamine, Naringenin, 2-Phenylacetamide, Luteolin, Catechin, and Quercetin were more likely to be regulated in the flavonoid biosynthesis pathway. The correlation analysis between flavonoid-related genes and metabolites showed that 48 Key Unigenes (4CL, AOMT, CHS, CHI, DFR, LAR…) were significantly correlated with metabolites (Fig. 8; P < 0.05, —r— >0.8; Table 4). These results could provide insight into the relationship between genetically regulated metabolites and the metabolic impact on gene expression (Fig. 9).