Transcriptome analysis identifies key genes involved in carotenoid biosynthesis in the flesh of red pummelo (Citrus maxima)

Abstract

Red-fleshed pummelo (Citrus grandis), a mutant of ‘Guanxi’, is an important economic plant rich in carotenoids. The transcriptome of ‘Guanxi’ pummelo was compared with red-fleshed pummelo, with the aim of identifying the key genes responsible for red flesh through RNA sequencing. In total, 44.11 and 47.71 Gb clean reads were respectively obtained after filtering in red and white flesh. Altogether 46,425 unigenes, with an average length of 1208 bp, were assembled. Compared with that of white flesh, 3722 genes were up-regulated, whereas 1421 genes appeared down-regulated in the red-fleshed pummelo. A GO analysis assigned 11,983 differentially expressed genes to 53 subcategories, and 3582, 4722 and 3679 differentially expressed genes were involved in biological processes, cellular components and molecular functions, respectively. The mevalonate pathway involved in terpenoid backbone biosynthesis was significantly different between ‘Guanxi’ pummelo and its red-fleshed mutant, and the differences were closely associated with carotenoid biosynthesis. The expression levels of the carotenoid biosynthetic pathway genes were different between ‘Guanxi’ pummelo and its red-fleshed mutant, with the upstream genes having higher expression levels in red-fleshed pummelo than in ‘Guanxi’ pummelo, whereas the downstream gene β-CHX was downregulated. Two bHLH transcription factors (CL4778.Contig1 and CL4778.Contig2) were of the same branch as CubHLH1, regulating carotenoid biosynthesis in Citrus unshiu Marc., and Unigene CL4778.Contig2 had expression differences in two kinds of pummelo, which suggested CL4778.Contig2 should be a candidate gene controlling carotenoid accumulation in red pummelo. The study would be beneficial to identify the major gene regulating the red flesh trait.

Keywords:

• Pummelo
• red flesh
• carotenoid biosynthesis
• transcriptome analysis
• bHLH

Introduction

Pummelo (Citrus maxima) is an important species of the Citrus genus and a common germplasm for citrus breeding. It is widely commercially grown in some countries, including Japan, Thailand, Malaysia and China, partially because pummelo contains plentiful antioxidant compounds, such as carotenoids, flavonoids and limonoids [1,2]. In pummelo, the flesh colour of the fruit is a quality trait that correlates with nutritional value. Carotenoids are the major pigments responsible for the colour accumulation in Citrus fruit [3–5]. Carotenoids impart the yellow, orange or red colouration that is characteristic of many fruit and required for optimum seed dispersal [6,7].

In higher plants, carotenoids are synthesized mainly in cell plastids [8]. The carotenoid biosynthetic pathway plays various roles in photosynthesis, photoprotection, development, stress hormone responses, and the production of various volatiles and signalling apocarotenoids [9]. The molecular mechanism of carotenoid biosynthesis is clear, and schematic diagrams of the main pathway steps have been built [3,4,9,10]. The first step involves the condensation of geranylgeranyl pyrophosphate, the C20 precursor, to build the colourless phytoene. This step is catalysed by a rate-limiting enzyme, phytoene synthase (PSY). Second, phytoene undergoes four consecutive desaturation reactions, catalysed by phytoene desaturase (PDS) and ζ-carotene desaturase (ZDS), and two isomerization reactions that are regulated by carotene isomerase (CRTISO) and 15-cis-ζ-carotene isomerase, to form lycopene. In addition, the cyclization of lycopene, as a key reaction, is mediated by two enzymes encoded by lycopene β-cyclase (β-LCY) and ε-LCY. Thus, two pathways exist, which result in the formation of β- and α-carotene. Lutein is formed from α-carotene catalysed by β-carotene hydroxylase (β-CHX) and ε-CHX. β-Carotene is first converted into zeaxanthin through a hydroxylation reaction catalysed by β-CHX, and then violaxanthin is generated through an epoxydation reaction catalysed by zeaxanthin epoxidase (ZEP). Finally, the end product of the pathway, abscisic acid, is generated from violaxanthin through a reaction mediated by 9-cis-epoxycarotenoid dioxygenase (NCED). Recently, the molecular bases of carotenoid biosynthesis and accumulation have been analysed. The accumulation of lycopene during fruit ripening results from the downregulation of β-LCY activity [11]. Lycopene appears to be the determining factor in the red colouration of papaya flesh [12]. Structural genes encoding enzymes of the carotenoid biosynthetic pathway, including PDS [13–15], PSY [14,16,17], ZDS [14], β-LCY [14], ε-LCY [14, 17], ZEP [14], CRTISO [14] and HYb [14], have been identified and characterized. A comprehensive available transcriptome dataset for pummelo has been obtained, and the associated molecular markers have been identified [18]. Importantly, genes relating to secondary metabolites were mainly regulated by transcription factors [19]. The carotenoid biosynthetic pathway is mainly regulated by the transcription factors MYB [20,21] and bHLH [22–24].

Next-generation sequencing technology is a powerful and cost-effective tool for genome sequencing, genome re-sequencing, mRNA expression profiling and DNA methylation analyses. Recently, de novo transcriptome sequencing has been widely used for functional gene discovery in non-model organisms for which reference-genome information is lacking [25]. In addition, this approach has also allowed the analyses of gene expression profiles using the quantification of short cDNA reads [26–28]. In this manuscript, RNA sequencing (RNA-Seq) using next-generation sequencing was employed to determine the transcriptome differences in ‘Guanxi’ pummelo and its red-fleshed mutant to discover the key genes responsible for red flesh in pummelo.

Materials and methods

Plant materials

The materials used in this study were ‘Guanxi’ pummelo and its red-fleshed mutant, both of which are widely planted in Fujian Province, China. The fruit were collected at the ripe stage, with three biological replicates, and were immediately frozen in liquid nitrogen and kept at –80 °C until further use. Total RNA was extracted from the fruit employing the Tiangen RNAprep Pure Plant Kit (Tiangen, China) in accordance with the manufacturer’s protocol. Electrophoresis on a 1.0% (w/v) agarose gel demonstrated the quality of total RNA. The concentration of total RNA was determined by NanoDrop (Thermo Fisher Scientific, Wilmington, DE, USA).

Illumina sequencing

The cDNA libraries of the ‘Guanxi’ pummelo and red-fleshed pummelo fruit were prepared in line with the manufacturer’s instructions for mRNA-Seq sample preparation (Illumina Inc, San Diego, CA, USA). The cDNA library products were addressed by Illumina paired-end sequencing technology with read lengths of 150 bp, using an Illumina HiSeqXten-PE150 instrument (Novogene Bioinformatics Technology Co. Ltd., Beijing, China). The raw sequence reads were stored in the National Center for Biotechnology Information (NCBI) SRA database.

Sequence data filtering, de novo assembly and annotation

The raw data from the images were assembled with Illumina GA Pipeline 1.6 by removing low-quality reads (reads with unknown ‘N’ sequences, adaptor sequence fragments and empty reads). Next, the de novo assembly of the transcriptome into unique-sequence unigenes was processed with Trinity software, a short-read assembly program [29]. In short, Trinity integrates reads with certain lengths of overlap to form longer fragments or contigs in the first place. Then, the reads are mapped back to contigs, and with paired-end reads, it becomes possible to detect contigs from the same transcript as well as the distances between these contigs. Finally, Trinity connects the contigs and generates sequences that can be extended on neither end. Such sequences are defined as unigenes [29]. Later on, a Basic Local Alignment Search Tool (BLAST)x alignment (e-value < 0.00001) was carried out on the unigenes from protein databases, involving NCBI non-redundant protein (Nr) and nucleotide sequence (Nt), Swiss-Prot, Kyoto Encyclopedia of Genes and Genomes (KEGG), Cluster of Orthologous Groups of proteins (COG), Interpro and Gene Ontology (GO). Also, the best alignments were employed to determine the sequence directions of the unigenes. If the results from the different databases showed a conflict between each other, then the priority order of Nr, Swiss-Prot, KEGG and COG was followed to determine the sequence direction. When a unigene was not aligned in any of the above databases, the software ESTScan was used to predict its coding regions and determine its sequence direction [30]. In the final step, using Nr annotations, the Blast2GO program was used to obtain the GO and KEGG annotations for the unigenes [31]. After the GO annotation was obtained for each unigene, WEGO software was used to classify the unigenes in terms of function and to determine the distribution of gene functions in the species at the macro level [32]. A Venn diagram of the annotated proteins was constructed using the online software venny (http://bioinfogp.cnb.csic.es/tools/venny/index.html).

Analysis of differential gene expression

To map the differentially expressed reads, the sequenced raw data were filtered to remove low-quality reads. For read annotations, clean read sequences were mapped to our transcriptome reference database, allowing no more than a two-nucleotide mismatch [29,33]. Gene coverage is the percentage of a gene covered by reads. This value is equal to the ratio of the base number in a gene covered by unique mapping reads to the total base number of that gene. For the gene expression analysis, the number of expressed sequences was calculated and normalized to the number of Reads per Kb per Million reads.

To compare the differences in gene expression, the read frequency in each library was statistically analysed according to the method of Audic and Claverie [33]. We used a false discovery rate < 0.001 and a |log2Ratio| > 1 as the threshold for judging the statistical significances of the gene expression differences [34]. Next, the differentially expressed genes (DEGs) were subjected to GO and KEGG enrichment analyses [35]. The genes related to carotenoid biosynthesis in KEGG pathways were collected and aligned to the unigenes from a transcriptome mixture of the fruit from ‘Guanxi’ pummelo and its red-fleshed mutant, using the BLASTx algorithm with an E-value of < 1 e − 5. The selected genes were: 1-deoxy-d-xylulose 5-phosphate synthase, 1-deoxy-d-xylulose 5-phosphate reductoisomerase, geranylgeranyl pyrophosphate synthase (GGPS), PSY, PDS, ZDS, CRTISO, β-LCY, ε-LCY, β-CHX, ZEP and NCED.

Bioinformatics analysis

The sequence alignments were conducted using the Vector NTI 10 software package (Thermo Fisher Scientific, Wilmington, DE, USA). The phylogenetic tree was constructed with MEGA5 software (Tokyo Metropolitan University, Hachioji, Tokyo, Japan) using the neighbour-joining method [36].

Results and discussion

Sequencing and sequence assembly

RNA-Seq efficiently and conveniently enables rapid access to all transcripts from biological tissues. After filtering, 47,712,090 reads from the wild-type pummelo and 44,111,638 from the red-fleshed mutant remained, with a Q20 percentage of 96.69%. Trinity software was used to link the valid data to 46,425 unigenes, totalling 5,60,96,761 bp. The average length of each gene was 1208 bp (Supporting Information Table S1), which was longer than in another study [18]. The N50 length was 1922 bp (Supporting Information Table S1). Accurate and effective assemblies usually have longer mean and N50 lengths [37]. BLASTX software was employed to annotate the 46,425 assembled unigenes using the Nr, Nt, Swissprot, KEGG, COG, Interpro, GO and other protein databases (Table 1). A total of 40,261 proteins were obtained after removing repeated proteins. Of the predicted proteins, 68.40% had the highest homology to Citrus sinensis, 25.68% to Citrus clementina and 5.92% had closest genetic relationships to other species (Supporting Information Figure S1). Pummelo belongs to the Citrus genus, and the predicted proteins were highly homologous to its species.

Figure 1. Differentially expressed genes between ‘Guanxi’ pummelo and red-fleshed pummelo. The genes were divided into three classes. Red dots indicate genes that were upregulated, i.e. the gene’s expression in red-fleshed pummelo was higher than in ‘Guanxi’ pummelo. Blue dots indicate genes that were downregulated, i.e. the gene’s expression in ‘Guanxi’ pummelo was higher than in red-fleshed pummelo. Grey dots indicate genes that were not differentially expressed. The horizontal coordinates refer to the average FPKM value (A) and the vertical coordinates refer to the log2ratio (M).

Table 1. Annotation summary of predicted protein numbers from various databases.

Values	Total	Nr	Nt	Swissprot	KEGG	KOG	Interpro	GO	Intersection	Overall
Number	46,425	34,658	39,422	23,261	23,476	24,536	25,223	21,885	11,858	40,261
Percentage	100%	74.65%	84.92%	50.10%	50.57%	52.85%	54.33%	47.14%	25.54%	86.72%

DEG analysis

To identify DEGs between ‘Guanxi’ pummelo and its red-fleshed mutant, putative DEGs were identified based on the FPKM values calculated from the read counts mapped onto the reference transcriptome. A total of 5143 unigenes were differentially expressed between ‘Guanxi’ pummelo and red-fleshed pummelo according to the comparison of expression levels with false discovery rate ≤ 0.001 and |log2Ratio| ≥ 1 (Figure 1). The highest log₂Ratio value reached 10.84, whereas the lowest value was –12.62. Using ‘Guanxi’ pummelo as a reference, 3722 upregulated unigenes (with higher levels of expression in the mutant red-fleshed pummelo) and 1421 downregulated unigenes (with higher levels of expression in ‘Guanxi’ pummelo) were identified. Significantly more unigenes were upregulated in its mutant than in ‘Guanxi’ pummelo.

The GO analysis assigned 11,983 DEGs to 53 subcategories. Among them, 3582, 4722 and 3679 DEGs were involved in biological processes, cellular components and molecular functions, respectively. Within the biological process category, the majority were related to ‘metabolic processes’ (1056 DEGs) and ‘cellular processes’ (1037 DEGs). Within the cellular component category, most were enriched in the subcategories ‘membrane’ (1097 DEGs) and ‘membrane part’ (1038 DEGs). Within the molecular function category, the greater proportion was involved in ‘binding’ (1649 DEGs) and ‘catalytic activity’ (1632 DEGs; Figure 2). A total of 27,024 DEGs were annotated, and 4142 were mapped onto 130 pathways in the KEGG database. Most DEGs were involved in ‘metabolic pathways’ (2874), and ‘genetic information processing’ (692). The KEGG pathways ‘plant hormone signal transduction’, ‘ABC transporters’, ‘benzoxazinoid biosynthesis’, ‘diterpenoid biosynthesis’ and ‘carotenoid biosynthesis’ showed the greatest differences between ‘Guanxi’ pummelo and its red-fleshed mutant.

Figure 2. GO classifications of DEGs. Unigenes were assigned to three categories: cellular components, molecular functions and biological process.

The pathway ‘plant hormone signal transduction’ includes the signal transduction pathways of jasmonate, salicylate, auxin, cytokinin, gibberellin and abscisic acid. These plant hormones are believed to regulate fruit development and shape [38–42]. The pathway ‘ABC transporters’ is mainly involved in plant defence and the symbiosis between plants and microorganisms [43]. The pathway ‘benzoxazinoid biosynthesis’ plays an important role in host resistance against microbial pathogens and insects, and in allelopathic effects. The three pathways are not directly related to secondary metabolism. It is possible that they can indirectly regulate diterpenoid or carotenoid biosynthesis. Terpenoid backbone biosynthesis has two pathways, the mevalonate and the MEP/DOXP pathway. The mevalonate pathway showed significant differences between ‘Guanxi’ and red-fleshed pummelo, whereas no significant differences in the MEP/DOXP pathway. The mevalonate pathway produces two five-carbon building blocks, isopentenyl pyrophosphate and dimethylallyl pyrophosphate [44], which are the basic building blocks of carotenoid biosynthesis [45]. They undergo several reactions, which results in the production of the carotenoid precursor, geranylgeranyl pyrophosphate [45]. Changes in the mevalonate pathway could activate the carotenoid biosynthesis pathway.

Unigenes in the carotenoid biosynthesis pathway and their expression pattern

The FPKM values of genes were determined, and their relative expression levels were calculated in the carotenoid biosynthetic pathway. In total, 42 unigenes were involved in carotenoid biosynthesis. Although 12 genes related to carotenoid biosynthesis were used in the BLASTx analysis, only 1-deoxy-d-xylulose 5-phosphate synthase (DXS), 1-deoxy-d-xylulose 5-phosphate reductoisomerase (DXR), GGPS, PSY, PDS, ZDS, CRTISO, β-LCY, ε-LCY, β-CHX, ZEP and NCED have homologous genes in the assembled unigene database. The phylogenetic trees were constructed using the deduced amino acid sequences of 42 unigenes in pummelo and structural genes regulating the carotenoid biosynthesis in other species (Supporting Information Figures S2–S13). The same structural genes were assigned to the same branch (Figures S2–S13). The FPKM values of structural genes were determined to evaluate their relative expression levels. The carotenoid biosynthetic pathway was exhibited in pummelo, as described by Guo et al. [22] (Figure 3 and Supporting Information Table S2). The upstream structural genes had greater expression levels in red-fleshed pummelo than in ‘Guanxi’ pummelo, except GGPS and PSY, whereas the downstream structural gene β-CHX had lower expression levels in red-fleshed pummelo (Figure 3). Only a few genes in other pathways were identified in this experiment. Genes homologous to ε-CHX were not identified in the unigene database.

Figure 3. Differences in the expression levels of genes in the carotenoid biosynthetic pathway.

Note: Arrows show the metabolic stream; left or upward arrows represent the genes catalysing the process; the light color indicates that these genes are in the assembled unigenes; the numbers represent the increased expression in red-fleshed pummelo compared with ‘Guanxi’ pummelo.

The key bHLH transcription factor regulating carotenoid biosynthesis

Previous studies demonstrated that transcription factors MYB and bHLH regulated carotenoid biosynthesis in Citrus species [22,23,46,47]. In this research, MYB protein CsMYB [47] and bHLH protein CubHLH1 [23] were used as the reference sequences for blastx analysis to identify the candidate genes regulating carotenoid biosynthesis. Five MYB transcript Unigenes were selected to construct the phylogenetic tree (Figure 4). Phylogenetic analysis showed five MYB transcript Unigenes were not assigned to the same branch as CsMYB involved in carotenoid biosynthesis, suggesting that MYB could not be the key gene for the red flesh trait in the study (Supporting Information Figure S14). bHLH transcription factors negatively regulate carotenoid accumulation in Solanum lycopersicum [23,46]. CubHLH1 has a similar function to Arabidopsis AIFs and was directly involved in carotenoid metabolism in citrus fruit [23]. Two bHLH transcriptor Unigenes (CL4778.Contig1 and CL4778.Contig2) were selected with scores over 1166 to construct the phylogenetic tree containing CsbHLH and the bHLH transcriptors from Arabidopsis and Citrus. CL4778.Contig1 and CL4778.Contig2 were assigned to the same branch as CubHLH1, which implies that the CmbHLH (CL4778.Contig2) in the carotenoid biosynthesis pathway had been isolated from the RNA-Seq database. In addition, the amino acid sequences of CL4778.Contig1 and CL4778.Contig2 were aligned with CubHLH1 (Supporting Information Figure S15). The deduced amino acid sequences alignment demonstrated that the typical bHLH domains of CL4778.Contig1 and CL4778.Contig2 were coincident with CubHLH1 [18]. When compared with red flesh, CL1085.Contig2 had a higher transcript level in ‘Guanxi’ pummelo, as they were up-regulated 360 times. This suggested that the Unigene CL4778.Contig2 was the candidate gene for carotenoid biosynthesis, thus providing a theoretical basis for developing new cultivars with higher carotenoid levels.

Figure 4. Phylogenetic relationships of the deduced amino acid sequences of bHLH proteins in pummelo and bHLHs in other species.

Conclusions

In this study, transcriptome analysis was employed to discover the genetic differences between ‘Guanxi’ pummelo and its red-fleshed mutant. Here, a total of 46,425 unigenes, with an average length of 1208 bp, were assembled, in which 40,261 were annotated. A total of 5143 unigenes were differentially expressed between ‘Guanxi’ pummelo and red-fleshed pummelo, and 3722 upregulated unigenes (with higher levels of expression in the mutant red-fleshed pummelo) and 1421 downregulated unigenes (with higher levels of expression in ‘Guanxi’ pummelo) were identified. The mevalonate pathway involved in terpenoid backbone biosynthesis was significantly different, which was closely associated with carotenoid biosynthesis. The expression levels of the upstream genes had higher expression levels in red-fleshed pummelo than in ‘Guanxi’ pummelo, whereas the downstream gene β-CHX was downregulated in the carotenoid biosynthetic pathway. Two bHLH transcript factors, CL4778.Contig1 and CL4778.Contig2, were of the same branch as CubHLH1, regulating carotenoid biosynthesis, and Unigene CL4778.Contig2 had significant expression difference in the two kinds of pummelo. Therefore, we suggest that the bHLH transcript factor regulates the carotenoid accumulation in pummelo.

Disclosure statement

The authors have declared that no competing interests exist.

Additional information

Funding

This research was financially supported by the QingHai Science and Technology Department [grant no. 2020-ZJ-972Q], West Light Talent Program of the Chinese Academy of Sciences, the Open Project of State Key Laboratory of Plateau Ecology and Agriculture, Qinghai University [grant no. 2018-KF-06].