Transcriptomic and targeted metabolome analyses revealed the regulatory mechanisms of the synthesis of bioactive compounds in Citrus grandis ‘tomentosa’

Materials

The ‘Huanglong’ variety of C. grandis ‘tomentosa’ was sampled from the planting base of Xiangxiu Huajuhong Industrial Co., Ltd (Linchen Town, Huazhou City, 22°49′11″N, 110°36′11″E). The fruits (F, 30-day-old), petals (P, just blooming), and tender leaves (L, 20-day-old) were sampled on March 2, 2018. The plant materials, including 10 fruits, 20 leaves, and 50 flowers, were obtained each of nine 10-year-old individuals, , and the samples of three individuals were obtained as biological repeats. Three biological replicates were collected from the tissues, frozen rapidly with liquid nitrogen, and stored at −80 °C for further analysis.

Extraction and determination of the flavonoids and coumarin contents of C. grandis

The extraction and determination of the flavonoid and coumarin contents of C. grandis were performed using a method developed by Han et al. (2014). Samples of the different tissues were dried to a constant weight in an oven at 60 °C (based on the traditional method), crushed, and sieved through a 40 µm mesh. The powder was used for extraction. For the extraction, 0.25 g of flower, fruit, and leaf powder samples, respectively, were put into a 25 ml extraction tube. Then, 25 ml of 60% ethanol solution was added and the mixture was well shaken in an ultrasonic water bath at 40 °C for 30 min. Thereafter, the mixture was centrifuged at 6,000 rpm for 10 min. Qualitative and quantitative analysis of four types of flavonoids in the supernatants were performed using high-performance liquid chromatography (HPLC), using a SPD-20a HPLC system (Shimadzu Corporation, Kyoto, Japan) coupled with a thermo ODS chromatographic column (250 mm × 4.6 mm, 5 µm) and a Lc-20at UV detector. Mobile phase A comprised pure methanol, whereas mobile phase B consisted of 1% glacial acetic acid aqueous solution. The detection conditions were as follows: gradient elution (0–10 min, 35% a; 10–30 min, 35% a–65% a; 30–35 min, 65% a; 35–42 min, 65% a–35% a); a flow rate of 1.0 ml/min; a column temperature of 40 °C; and an injection volume of 10 µL. Naringin and naringenin were detected at 283 nm, while rhoifolin, apigenin, and coumarin were detected at 340 nm. Quantitative analysis of the bioactive compounds was performed based on the peak area of the standard curve of each compound. Naringin standard (no.110722-200810) was provided by the China Institute for the control of pharmaceutical and biological products, while naringin (no.120824), rhoifolin (no.120711), and apigenin (no.120906) standards were provided by Sichuan Vicki Biotechnology Co., Ltd. All standards used for HPLC were of high purity (≥98%). The naringin, rhoifolin, naringenin, apigenin, and coumarin contents of ECG were recorded as the average values of three biological repetitions.

Extraction and determination of the volatile compounds in C. grandis

Approximately 0.1 g of flower, fruit, and leaf powder samples were extracted using 10 mL of absolute alcohol at 30 °C for 30 min, followed by centrifugation. The supernatant was dehydrated using anhydrous sodium sulfate at 4 °C for 24 h. The mixture was then centrifuged at 8000 rpm for 10 min and then filtered through a 0.22 µm filter membrane. The volatile organic compound profile of the supernatant was determined by on-line gel permeation chromatography-gas chromatography-mass spectrometry (Han et al., 2018). Principal Components Analysis (PCA) and heatmap analysis and Venn diagram construction were performed using the omicshare (http://www.omicshare.com/tools) online tools.

Extraction and high-throughput sequencing of RNA from C. grandis

Total RNA was extracted from the tissues using the RNeasy plant Mini Kit (Qiagen, Valencia, CA, USA) according to the manufacturer’s instructions. RNA was treated with RNase-free DNase I (Takara Bio, Ostu, Japan) to remove residual DNA. RNA quantity and integrity were assessed using electrophoresis (1% agarose gel) and an Agilent 2100 Bioanalyzer (Agilent Technologies, Santa Clara, CA, USA), respectively. The mRNA was enriched and purified using oligo magnetic beads and then broken into 200 NT short fragments using a fragmentation buffer. Thereafter, first strand cDNA was synthesized using mRNA as the template. The second strand was synthesized immediately. The QiaQuick PCR extraction kit (Qiagen, Venlo, Holland) was used to purify the cDNA fragment, repair base ends, and connect the Illumina sequencing adapters. The size of the ligated product was determined using agarose gel electrophoresis. Finally, the product was amplified using PCR (Wang et al., 2016) and sequenced on an Illumina Hiseq 4000 (Illumina, San Diego, CA, USA) platform. The sequencing length was double ended 150 bp. All sequence data were deposited at the Genome Sequence Archive of the Beijing Institute of Genomics Data Centerunder accessions CRA006910.

Transcriptome data analysis: annotation and differential expression analysis

First, the raw reads were processed to remove joint sequence, low-quality reads (the base number of an alkali with mass value Q ≤ 20 accounts for more than 50% of the whole read) and reads with N ratio greater than 10%. The read was compared to the ribosomal database to remove rRNAs and obtain high-quality reads using bowtie v.2.2.8 software (Langmead et al., 2009). Thereafter, the clear reads were aligned against C. grandis genome version 1 (http://citrus.hzau.edu.cn/download.php) using Tophat2 v.2.1.1 software (Kim et al., 2013; Wang et al., 2017). The genomic reads were assembled by comparing them using Cufflinks (Trapnell et al., 2010), and then the assembly results of multiple samples were combined using cuffmerge. Transcripts with assembly errors were filtered to generate a unique annotation file for subsequent analysis. Genes with size ≥ 200 bp and exon ≥ 2 were mapped and functionally annotated to NCBI RefSeq, Kyoto Encyclopedia of Genes and Genomes (KEGG), and Swissprot databases.

The expressed values of all genes were calculated and normalized to fragments per kilobase of transcript per million fragments mapped (FPKM). Differential expression analysis of genes between two samples was performed using the DESeq package in R software (R Core Team, 2021). Genes with false discovery rate <0.05 and fold change —log2— ≥ 1 were considered differentially expressed genes (DEGs) between samples. Pearson’s correlation analysis was performed to determine the relationship between two samples, which was visually display in the form of a heatmap. Transcriptome power analysis using functions in RNASeqPower (Corley et al., 2017)

Gene ontology and KEGG pathway enrichment analysis

The identified DEGs were aligned against gene ontology (GO) and KEGG databases for functional annotation and pathway analysis, using KOBAS (http://kobas.cbi.pku.edu.cn/) and omicshare (http://www.omicshare.com/tools) online tools, respectively.

Quantitative reverse transcription PCR

To verify the RNA-seq results, nine genes were selected for quantitative reverse transcription PCR (qRT-PCR) analysis using Actin as the internal reference gene (Tian, Fan & Zeng, 2021). The reactions were performed on Bio-Rad CFX96 Touch (Bio-Rad, CA, USA) using 5 µL of 2 × ChamQ Universal SYBR qPCR Master Mix (Vazyme Biotech Co., Ltd., Nanjing, China) containing 1.5 µL of cDNA template, 0.4 µL of each primer (10 µmol/ µL), and 2.7 µL of nuclease-free water. PCR conditions and Melting curve analyses are described by Huang et al. (2017). Relative gene expression levels were normalized to that of Actin and determined using the 2^−ΔΔCT method (Livak & Schmittgen, 2001). All reactions were performed in triplicate.

Data analysis

All data are presented as the mean of several values. Differences in the bioactive contents between the plant tissues were determined using Duncan’s new compound polar difference method and SigmaPlot v.11 software (Systat Software Inc., CA, USA). Means were considered statistically significant at p < 0.05.

Flavonoid and coumarin contents of different tissues of C. grandis

The naringin (418.32 mg/g), naringenin (2.73 mg/g), and coumarin (0.42 mg/g) contents of the tender fruits were significantly higher than those of the tender leaves and petals, whereas the tender leaves had significantly higher levels of rhoifolin (57.70 mg/g) and apigenin (2.41 mg/g). Specifically, the naringin content of the tender fruits was more than six times those of the tender leaves and petals, whereas the rhoifolin content of the young leaves was more than 10 times those of the other two tissues (Table 1, Table S1).

Volatile oil content of different tissues of C. grandis

A total of 49 volatile compounds were identified in the tender leaves, fruits, and petals of ECG, among which 15 were monoterpenes, 16 were sesquiterpenes, one was triterpene, and 17 were fatty acids (Table S2). Additionally, 26 volatile compounds were common to the petals, fruitlets, and tender leaves, three were found only in the leaves, two only in the fruits and five only in the petals (Fig. 1A). PCA results showed that the three biological replicates from different tissues were well clustered, indicating good reproducibility. The good separation of different organizations indicates that there were differences in the compositions of volatile oils among different organizations (Fig. 1B). Moreover, among 38 volatile compounds identified the petals, the main compounds were β-pinene (7.53%), β-myrcene (15.32%), γ-terpinene (8.23%), nerolidol (20.57%), and farnesol (14.85%). Among 36 volatile compounds identified in the fruit, the main compounds were β-myrcene (24.35%), γ-Terpinene (20.70%), and germacrene D (18.46%). Furthermore, a total of 40 components were identified in the leaves, among which the main components were β-pinene (8.27%), β-myrcene (10.31%), γ-terpinene (8.76%), ficusin (8.52%), and squalene (9.25%) (Table S2). Furthermore, heat map analysis results showed that the volatile oil contents in different tissues could be divided into four categories. The first category mainly exists in the petals, such as 2,4-Di tert butylphenol and Nerolidyl acetate. The second category mainly exists in flowers, such as Myrcene. The third category mainly exists in leaves, such as a-Ocimene, and the fourth category is distributed in petals and leaves, such as 3-Thujene (Fig. 1C).

Sequence analysis of transcriptome in different tissues

The amount of transcriptome sequencing data of per sample was 6 G, and the saturation analysis showed that per sample had reached saturation (Fig. 2). A total of 55.02 G HQ clean data were obtained from nine samples (Table 2), following the removal of low-quality reads. The Q30 values (i.e., the sequencing accuracy was 99.9%) were >95.0%. The percentages of high-quality reads were >96%, indicating that the sequencing quality was high. Subsequently, the clean reads were aligned against the pomelo genome (version 1), of which an average of 86.83% of reads from nine libraries were successfully mapped to the grapefruit genome. The statistical power of this experimental design, calculated in RNASeqPower, is 0.76% (Table S3). A total of 31,215 genes were obtained, among which 21,672 were known genes in the genome and 1,092 were novel genes (Table 2).

GO analysis of the novel genes showed that a total of 180 new transcripts were annotated in 15 terms in “biological processes,” 11 terms in “cellular components,” and six terms in “molecular functions” (Fig. 3). Specifically, 47.2, 41.7, and 37.2% of the genes were annotated in “metallic process,” “cellular process,” and “single organization process,” respectively, in ‘biological processes.’ Additionally, 31.1, 31.1, and 21.7% of the genes were annotated in “cell part,” “cell,” and “organelle,” respectively, in ‘cellular components.’ Moreover, 48.9% and 28.3% of the genes were annotated in “catalytic activity” and “binding” in ‘molecular functions.’

Analysis of sequence repeatability of the transcriptome in different tissues

Principal component analysis (PCA) was performed to visualize the gene expression profiles of samples from the different tissues of ECG. The results show that the three biological replicates from each tissue clustered together, indicating that the sequencing repeatability was high. Additionally, the gene expression profiles of the leaves formed a distinct cluster in the second principal component PC2, whereas the gene expression profiles of the fruitlets formed distinct clusters in the first principal component PC1, indicating that there are significant differences in gene expression in different tissues (Fig. 4A). Furthermore, Pearson’s correlation analysis was performed to determine the relationships between the expression profiles of different tissues, using the FPKM values of the genes (R² > 0.8 indicates a significant correlation between the two samples). The results show that the Pearson correlation coefficient among three replicates of the same tissue was greater than 0.95 (Fig. 4B), indicating that the sequencing had high repeatability and could be used for subsequent differential expression analysis.

The expression profiles of twelve genes were amplified using RT-qPCR to confirm the accuracy and reproducibility of the RNA-seq results, using specific primers (Table S4). The results of the RT-qPCR and RNA-seq show similar expression profiles for the twelve genes (Table S5). Additionally, the Pearson’s correlation coefficient between the RNA-Seq and RT-qPCR data was 0.9462 (P < 0.0001), indicating the reliability of the RNA-Seq data (Fig. 5).

Analysis of gene differential expression in different tissues of C. grandis

Differential expression analysis was performed to identify DEGs between the three tissues. A total of 11,514 DEGs were identified in a pairwise comparison of three tissues, among which 5,605, 10,466, and 9,942 DEGs were identified in the L vs F, L vs P, and P vs F comparison groups, respectively. Additionally, 1,081 and 3,804 genes, 2,820 and 7,646 genes, and 6,660 and 3,282 genes were upregulated and downregulated in the L vs F, L vs P, and P vs F comparison groups, respectively (Fig. 6A). Among the differential genes, 2,420 genes were differentially expressed in the three tissues, 405 were differentially expressed only between tender fruits and leaves, 1,281 genes were differentially expressed only between petals and leaves, and 1,021 genes were differentially expressed only between tender fruits and petals (Fig. 6B).

Functional analysis of differentially expressed genes

KEGG pathway analysis of the DEGs showed that a total of 2,566 DEGs were significantly enriched in 43 KEGG pathways, among which “metabolic pathways” (1,228 DEGs, 47.86%) and “biosynthesis of secondary metabolites” (787 DEGs, 30.67%) were the most significantly enriched pathways (Fig. 7, Table S6). In addition, the pathways involved in flavonoid synthesis, such as “Flavonoid biosynthesis” (46 DEGs, 1.80%), “Flavone and Flavonol biosynthesis” (six DEGs, 0.23%), as well as pathways involved in the synthesis of volatile oil compounds, such as “Dieterpenoid biosynthesis” (32 DEGs, 1.24%), “Monoterpenoid biosynthesis” (19 DEGs, 0.74%), and “Fatty acid metabolism” (43 DEGs, 1.68%), were enriched significantly (Fig. 7, Table S6). The results indicate that the DEGs are mainly involved in the biosynthesis of flavonoids and volatile oils.

Gene expression of flavonoids and coumarin synthesis pathway

The expression profiles of genes related to naringin, rhoifolin, naringenin, and apigenin synthesis were examined in the different tissues of ECG (Fig. 8, Table S7). The results show that five phenylalanine enzyme (PAL) genes, which are involved in phenylpropane biosynthesis, were differentially expressed in the three different tissues, with the petals having the highest expression of three of the genes. Additionally, five 4CL genes regulating cinnamic acid production of cinnamoyl-COA were differentially expressed in the three tissues, among which two were upregulated in the fruitlets. Moreover, two CYP73A genes involved in regulating the synthesis of p-coumarinyl-COA were significantly upregulated in the fruitlets. Furthermore, five chalcone synthase (CHS) genes involved in regulating the production of grapefruit ligand chalcone were differentially expressed in the three tissues. One CHI1 gene involved in regulating naringenin synthesis was significantly upregulated in the fruitlets, which is evidenced by the higher naringenin content of the fruitlets. Moreover, naringenin is a precursor of apigenin, which in turn is a precursor of rhoifolin. Additionally, the fruitlets had higher expression levels of one 1,2Rhat gene (Cg1g023820), which is involved in the synthesis of naringin. Furthermore, 15 BGLU genes, which are involved in coumarin synthesis, were differentially expressed between the three tissues, with the fruitlets having significantly higher expression of two BGLU genes (Fig. 8, Table S7).

Gene expression of the terpenoids synthesis pathway

The biosynthesis of terpenoids comes from the main chain formed by the cytosol and plastid pathways, and further forms a series of substances such as monoterpenes, sesquiterpenes, and diterpenes. This process involves many key enzymes. A total of 29 DEGs of encoding key enzymes are involved in the regulation of terpenoid backbone biosynthesis, including acetyl CoA thiolase (AACT), 3-hydroxy-3-methylglutaryl CoA synthase (HMGS), hydroxymethylglutaryl CoA reductase, mevalonate kinase, mevalonate 5-diphosphate decarboxylase, 1-deoxy-d-xylose-5-phosphate synthase, cytidine 4-diphosphate-2-c-methyl-d-erythritose kinase (CMK), 2-c-methyl-d-erythritol-2,4-cyclophosphate synthase, 1-hydroxy-2-methyl-2-(E)-butenyl-4-diphosphate synthase (HDS), and 1-hydroxy-2-methyl-2-(E)-butenyl-4-pyrophosphate reductase (HDR) (Table S8, Fig. 9). Further analysis showed that 19 encoding key enzyme genes were differentially expressed in the monoterpene biosynthesis pathway, while 26 encoding key enzyme genes were differentially expressed in the sesquiterpene and triterpene biosynthesis pathway (Table S8, Fig. 9).