Transcriptome analysis provides new insight into the molecular mechanisms in regulating carotenoid accumulation in pericarp of fruits from two pummelo cultivars Citrus grandis (L.) Osbeck

Background Carotenoid improves fruit external quality and benets to human health. Although the mechanism underlining carotenoid metabolism is clear, the regulation of carotenoid accumulation remains poorly understand. Recently, ‘Huangbao’ pummelo was selected from the bud mutation of ‘Guanxi’ pummelo (yellow pericarp). The pericarp of mutant-type fruits is characterized by red colour development during fruit ripening stage. The aim of the study is to explore the molecular mechanisms in regulating pericarp colouration based on transcriptomic prole analysis of the two cultivars. Results Lycopene content in the pericarp of ripe fruits from ‘Huangbao’ (160.29 µg g-1 FW) was signicantly higher than that of from ‘Guanxi’ (0.91µg g-1 FW). Lycopene is a main contributor for pericarp colouration. 929 differentially expressed genes (DEGs) were indentied through transcriptional proling in the pericarp of ripe fruits from the two cultivars. Expression levels of the Genes encoding for fructose-bisphosphate aldolase, enolase and pyruvate kinase isozyme were signicantly higher in ‘Huangbao’ than that of in ‘Guanxi’, which possibly promote 3-phosphate-glyceraldehyde and pyruvic acid biosyntheses. This two productions act as precursors for the biosyntheses of Isopentenylpyrophosphate (IPP) and Dimethylallyl pyrophosphate (DMAPP) in Methylerythritol phosphate (MEP) pathway. Genes ecoding for phytoene synthase and prolycopene isomerase were unexpected down-regulated in ‘Huangbao’. Gibberellin-2-β-dioxygenase (GA2OX) gene was specially expressed in ‘Huangbao’, while the transcript amounts of genes related to auxin response and auxin transport as well as several negative regulators in ethylene signaling in ‘Huangbao’ were all down regulated. Additionally, Transcription factors (DELLA, NAC, MADS-box, WRKY) and post translational modication proteins (histone acetyltransferase and E2 ubiquitin-conjugating enzyme) involved in regulating ethylene and carotenoid metabolisms were also differentially expressed in the two pummelo cultivars. Conclusions The main differences in the carotenoid metabolism in the pericarp of ripen fruits from ‘Huangbao’ and ‘Guanxi’ were observed in the carotenoid precursor biosynthetic pathway. Differentially expressed genes in gibberellin, auxin and ethylene metabolisms and plant hormone signaling pathways, as well as several transcription factors and post translational modication protein genes involved in regulating ethylene and carotenoid metabolisms provide targets for further exploration in revealing mechanisms underlying fruit colour development. Abstract Background: Carotenoid improves fruit external quality and benets to human health. Although the mechanism underlining carotenoid metabolism is clear, the regulation of carotenoid accumulation remains poorly understand. Recently, ‘Huangbao’ pummelo was selected from the bud mutation of ‘Guanxi’ pummelo (yellow pericarp). The pericarp of mutant-type fruits is characterized by red colour development during fruit ripening stage. The aim of the study is to explore the molecular mechanisms in regulating pericarp colouration based on transcriptomic prole analysis of the two cultivars. Results: Lycopene content in the pericarp of ripe fruits from ‘Huangbao’ (160.29 µg g -1 FW) was signicantly higher than that of from ‘Guanxi’ ( 0.91µg g -1 FW). Lycopene is a main contributor for pericarp colouration. 929 differentially expressed genes (DEGs) were indentied through transcriptional proling in the pericarp of ripe fruits from the two cultivars. Expression levels of the Genes encoding for fructose-bisphosphate aldolase, enolase and pyruvate kinase isozyme were signicantly higher in ‘Huangbao’ than that of in ‘Guanxi’, which possibly promote 3-phosphate-glyceraldehyde and pyruvic acid biosyntheses. This two productions act as precursors for the biosyntheses of Isopentenylpyrophosphate (IPP) and Dimethylallyl pyrophosphate (DMAPP) in Methylerythritol phosphate (MEP) pathway. Genes ecoding for phytoene synthase and prolycopene isomerase were unexpected down-regulated in ‘Huangbao’. Gibberellin-2-β-dioxygenase (GA2OX) gene was specially expressed in ‘Huangbao’, while the transcript amounts of genes related to auxin response and auxin transport as well as several negative regulators in ethylene signaling in ‘Huangbao’ were all down regulated. Additionally, Transcription factors (DELLA, NAC, MADS-box, WRKY) and post translational modication proteins (histone acetyltransferase and E2 ubiquitin-conjugating enzyme) involved in regulating ethylene and carotenoid metabolisms were also differentially expressed in the two pummelo cultivars. Conclusions: The main differences in the carotenoid metabolism in the pericarp of ripen fruits from ‘Huangbao’ and ‘Guanxi’ were observed in the carotenoid precursor biosynthetic pathway. Differentially expressed genes in gibberellin, auxin and ethylene metabolisms and plant hormone signaling pathways, as well as several transcription factors and post translational modication protein genes involved in regulating ethylene and carotenoid metabolisms provide targets for further exploration in revealing mechanisms underlying fruit colour development. Indole-3-acetic acid; carotenoid EOLs: Ethylene-overproducer1 (ETO1)/ETO1-like;


Background
It was found that Carotenoids are involved in developmental process in organisms, and specially synthesized in plants and microorganisms. As one kind of light harvesting pigments, Carotenoids impact the photosynthetic activity in plants [1,2]. Carotenoids elevate oxidation resistance in plants due to their roles in scavenging active oxygen species, and preventing their biosynthesis [3,4]. Carotenoids affect phytohormone biosynthesis, for examples, Isopentenylpyrophosphate (IPP) is the common precursor for carotenoid and gibberenllin (GA) biosyntheses [5]; the intermediate product in the downstream metabolic pathway of carotenoids is the precursor for abscisic acid biosynthesis [6,7]. Carotenoid composition determines plant tissue colours, and more important, carotenoid components are bene cial for hunman health [8]. Due to the fact that carotenoids participating in metabolic processes in plants and humans, much extensive attention is concentrated on exploring the knowledge of carotenoid metabolism.
The whole process of carotenoid biosynthesis can be devided into two periods: the previous and following carotenoid precursor biosynthesis. IPP and Dimethylallyl pyrophosphate (DMAPP) are synthesized through mevalonic acid and methylerythritol phosphate (MEP) pathways, and then act as precursors for carotenoid biosynthesis. MVA and MEP pathways exist in the cytosol and the plastid of plants, respectively. Carotenoids are synthesized mainly through MEP pathway in plants [5]. IPP and DMAPP are catalyzed by Geranylgeranyl pyrophosphate synthase to form geranylgeranyl diphosphate (GGPP). Phytoene synthase ( PSY) catalyzes condensation of two molecules of GGPP, forming the phytoene. Phytoene is desaturated by phytoene desaturase (PDS) to yield ζ-carotene, and then the desaturation of ζ-carotene is catalyzed by ζ-carotene desaturase (ZDS) to produce tetra-cis-lycopene. Carotenoid isomerase (CRTISO) converts tetra-cis-lycopene to all-trans-lycopene. Lycopene β-cyclase ( β-LCY) catalyzes the cyclation of all-trans-lycopene to form β-carotene, however, when the cyclation of alltrans-lycopene is catalyzed by β-LCY and lycopene ε-cyclase (ε-LCY) , α-carotene is formed. β-ring carotene hydroxylase (BCH) and ε -ring carotene hydroxylase (E CH) catalyze the hydroxylation of αcarotene, yielding the lutein. β-cartoene is hydroxylated only by BCH, converting to zeaxanthin. Zeaxanthin epoxidase converts zeaxanthin to antheraxanthin and violaxanthin. Violaxanthin is converted to neoxanthin by neoxanthin synthase at the last step [9].
The modulation of carotenoid accumulation is not limited in the carotenoid metabolic pathway. The carotenoid distribution in cell organelles also impacts carotenoid levels in plants. Chromoplasts act as the sink for carotenoids, their structure changes parallel carotenoid accumulation. Most chromoplasts in mature green fruits of red-eshed papaya (Carica papaya L. cv. Pococí) were characterized by widely undifferentiated stroma, then the plastoglobules could be observed in chromoplasts during fruit development, at ripening stage, crystal-type lycopene formed in numerous chromoplasts, meanwhile, plastoglobules and tubular structures could still be indenti ed [10]. The nal structural characteristics of chromoplasts are correlated with carotenoid composition. For example, lycopene accumulation in pulp of fruits from red-eshed 'Cara Cara' (Citrus sinensis L. Osbeck) and Star Ruby (Citrus paradise Macf.) was linked to the formation of needle-like crystal in chromoplasts, β, β-xanthophylls accumulation in 'Navel orange' (Citrus sinensis L. Osbeck) fruit pulp led to the formation of plastoglobule in chromoplasts, however, the elevated colourless carotene contents in 'Pinalate' (Citrus sinensis L. Osbeck) fruit pulp were associated with the development of numerous vesicles in chromoplasts [11]. The component analysis of chromoplasts in tomato (Solanum lycopersicum) fruits showed that the increased plastid-localized lipid monogalactodiacylglycerol content in chromoplast membrane, was correlated with the promoted capacity for β-carotene and lycopene in chromoplasts [12]. These reports indicate that the concurrent changes in chromoplast structures with carotenoids may offer feedback information for carotenoid biosynthetic metabolism in plants. Some other metabolic pathways are capable of impacting carbon skeleton supply for carotenoid biosynthesis or controlling carotenoid degradation, consequently, in uence carotenoid accumulation. [7]. Environmental factors also play important roles in regulating carotenoid metabolism. Light quality and illumination time regulate carotenoid contents and transcript levels of the genes involved in carotenoid biosynthesis and degradation [13,14]. In an addition, soil [15] and mineral elements [16][17][18] also exhibite effects in modulating carotenoid contents in plants. The uctuations of carotentoid contents in plants cultivated in different geographical conditions may be caused by the in uences of external environmental factors at least in part [19]. Even the carotenoid biosynthesis in fruits during storage is also in uenced by environment. The previous reports have been showed that high storage temperature (12-20 °C) conditioning signi cantly promoted carotenoid biosynthsis in fruits when compared to low storage temperature (4 °C) [20,21]. Therefore, the approaches in regulation of carotenoid accumulation are not limited in the core carotenoid metabolic pathway, which can be extend to the carotenoid precursor metabolic pathway and the carotenoid distribution in plastids. The above states can partially explain why the transcript levels of genes involved in core carotenoid metabolism are not always coincident with carotenoid accumulation [22,23].
it is necessary to carry out experiments using plant materials with different colour development to explore the mechanisms in regulating carotenoid accumulation, due to the fact that both of genetic background and environmental conditions in uence carotenoid metabolism. 'Huangbao' pummelo was derived from the bud mutation of 'Guanxi' pummelo. The pericarp of the mutant-type fruits is characterized by the red colour development during fruit ripening stage, while the wild-type pericarp still keeps yellow even at the fully ripe stage. The aim of the present study is to explore molecular metabolisms in regulating carotenoid accumulation based on transcriptome analysis in the pericarp of ripe fruits from 'Huangbao' and 'Guanxi' pummelo.

RAN-sequencing and assembly of unigenes
Six cDNA libraries were created using 'Huangbao' and 'Guanxi' samples (i.e. each cultivar with three independent biological replicates). Raw reads were generated from the Illumina sequencing, and clean reads were obtained by ltering the adaptor sequences and the low quality reads. The average GC content of clean reads was nearly the same (45%), and the average Q30 base percentage was more than 92%.
Those high quality reads were de novo assembled using trinity method, generating 97377 unigenes with an mean length of 730 nt, and a N50 length of 1188 nt. The size distribution of unigenes was showed in Fig. 1. More than 78% of the clean reads per sample could mapped on to the assemlbled unigenes ( Table  1 in Supplementary Files).

Functional annotation of unigenes
A total of 55820 (57.32%) all unigenes and 826 (88.91%) differentially expressed unigenes (DEGs) were annotated based on sequence alignment in 8 public databases, including NCBI non-redundant protein sequences (NR), Protein family (Pfam), a manually annotated and reviewed protein sequence database (Swiss-Prot), gene ontology (GO), Clusters of Orthologous Groups (COG), euKaryotic Orthologous Groups (KOG), Evolutionary Genealogy of Genes: non-supervised orthologous groups (eggNOG) and Kyoto Encyclopedia of Genes and Genomes (KEGG) ( Table 2 in Supplementary Files).
All unigenes were annotated within GO database into three main categories, including biological process (BP) (20 terms), molecular founction (MF) (17 terms) and cellular component (CC) (19 terms) (Fig. 3). In the BP category, the main enriched GO terms included metabolic process, cellular process, response to stimulus, biological regulation, localization and signaling, which indicated that the metabolic activities still maintained high levels in the pericarp even at fruit ripening stage. In case of MF, most of the genes belonged to catalytic activity, binding, transporter activity, nucleic acid binding transcription factor activity, enzyme regulator activity, and protein binding transcription factor activity, which suggested that those genes possibly involved in regulation of metabolisms. As far as CC, the predominant subcategories were cell part, organelle, membrane, macromolecular complex, membrane part, membrane-enclosed lumen, cell junction and symplast.
Functional annotation of unigenes in COG, KOG and eggNOG databases indicated that large number of unigenes were calssi cated into several functional categories, including general function prediction only, posttranslational modi cation, carbohydrate transport and metabolism, signal transduction mechanisms, translation and transcription (Additional le 1). It was also found that large number of unigenes were encriched in carbon metabolism and ribosome pathways based on the KEGG pathway enrichment analysis (Additional le 2).

Global detection of DEGs
The comparison of transcriptomes of the two pummelo cultivars identi ed 929 DEGs showing a signi cant change in expression level (FDR<0.01 and fold change ≥ 2), including 326 up-regulated genes and 603 down-regulated genes (additional le 3). The fold changes of the majority of those DEGs were distributed in the range of 2-10. It was investigated that 130 and 150 DEGs specially expressed in the pericarp of ripe fruits from 'Huangbao' and 'Guanxi', respectively. 90% of those specially expressed DEGs exhibited a low transcript level (FPKM < 5). Cluster analysis of the total DEGs based on log2(FPKM) was provided in Fig. 4.
GO enrichment analysis of DEGs GO enrichment analysis was performed using a threshold value of p≤ 0.05. In the group of up-regulated genes, there were 42, 14 and 37 subcategories were signi cantly enriched in BP, CC and MF, respectively, while in the down-regulated gene group, there were 55, 13 and 58 subcategories were signi cantly enriched in BP, CC, and MF, respectively. TopGO analysis showed that cell surface receptor signaling pathway (GO: 0007166) is one of the most signi cantly enriched subcategories in the BP category. In the case of CC, plastid (GO: 0009536), chloroplast (GO: 0009507), plastid part (GO: 0044435) and chloroplast part (GO: 0044434) were the most signi cantly enriched subcategories. As far as the MF, catalytic activity (GO: 0003824) and DNA binding (GO: 0003677) were the highest represented subcategories. The top ten classes from each categories were displayed in Additional les 4-6.
In addition to the above results, other signi cantly enriched subcategories were also obtained. In BP category, almost all of the genes calssi ed into translation (GO: 0006412) and isoprene biosynthetic process (GO: 0043612) were specially identi ed in the pericarp of ripe fruits from 'Huangbao'. Other relevant over-represented subcategories included: light harvesting (GO: 0009765), metabolic process (GO: 0008152), oxidation-reduction process (GO: 0055114), auxin transport (GO: 0010540) and transmembrane receptor protein serine/threonine kinase signaling pathway (GO: 0007178). As far as the MF category, most of the genes categorized into structural constituent of ribosome (GO: 0003735) and ribulose-1,5-bisphosphate carboxylase/oxygenase activator activity (GO: 0046863) only expressed in the pericarp of ripe fruits from 'Huangbao'. Monooxygenase activity (GO: 0004497) and potassium ion binding (GO: 0030955), carotenoid isomerase activity (GO: 0046608), metal ion binding (GO: 0046872) and hydroxymethylglutaryl-CoA reductase activity (GO: 0042282) were also the signi cantly enriched groups.

DEGs involved in carotenoid and its precursor metabolisms
Phytoene synthase and prolycopene isomerase are two enzymes that catalyze chemical reactions previous to the lycopene biosynthetic step. In this study, expression levels of genes ecoding for the two enzymes were signi cantly lower in the pericarp of ripe fruits from 'Huangbao' compared with 'Guanxi'.
Three genes (c36907.graph_c0, c54686.graph_c0, c38780.graph_c0) involved in MVA pathway were down regualted in 'Huangbao' and annotated as Acetyl-CoA acetyltransferase, Hydroxymethylglutaryl-CoA synthase and 3-hydroxy-3-methylglutaryl-coenzyme A reductase, respectively. Transcript levels of the three genes indicated that MVA pathway activity in 'Huangbao' was lower than that in 'Guanxi'. The above genes were listed in Table 3 (in Supplementary Files.

DEGs involved in Other metabolic pathways
The DEGs analyses showed that all of the genes involved in biosynthesis and degradation of fatty acid were down regulated in 'Huangbao'. In contrast, most of the genes annotated as 40S and 60S ribosomal proteins were up regulated in 'Huangbao'. In addition, all of the genes encoding for elongation factor and cytochrome c oxidase subunit were also up regulated. Those genes exhibited distinct differences in transcript levels were listed in Additional le 7.
DEGs associated with signal transduction DEGs related to auxin and ethylene responsive proteins were down regulated in 'Huangbao'. Four genes were identi ed for auxin responsive proteins, one of them was annotated as putative indole-3-acetic acidamido synthetase (GH3). In plants, this protein plays an important role in balancing active and inactive auxin contents through by catalyzing the synthesis of indole-3-acetic acid-amido. Two genes (c67065.graph_c0b, c69145.graph_c0) were not identi ed in 'Huangbao', one of which (c67065.graph_c0b) was annotated as EIN3-binding F-box protein 1 (EBF1), a negative factor in ethylene signaling pathway, while the other one encoded ethylene insensitive 3 (EIN3), a positive factor in ethylene signal transduction ( Table 5 in Supplementary Files).
Many DEGs were found to be associated with other signaling pathways. Those genes were mainly annotated as calmodulin, tyrosine kinase, leucine-rich repeat receptor-like protein kinase, cysteine-rich receptor-like protein kinase, serine/threonine-protein kinase, histidine kinase, rapid alkalinization factor, protein phosphatase and non-speci c lipid-transfer protein. It was observed that a largest number of genes were identi ed for serine/threonine-protein kinase. Most of these genes involved in other signaling pathways were down regulated in 'Huangbao' (Table 6).

Differentially expressed transcription factors
Some transcription factors related to fruit ripening process, phytohormone signal transduction, pigment metabolism and terpene biosynthesis, such like ethylene response transcription factor, MYB, DELLA, MADS-box and WRKY, were down regulated in 'Huangbao'. However, most of the NAC domain-containing protein genes showed up regulation ( Table 7 in Supplementary Files).
DEGs participated in post-translational modi cation Transcript levels of genes encoding for histone acetyltransferase and monothiol glutaredoxin were up regulated in 'Huangbao', in contrast, 14-3-3-like protein gene showed down regualtion. It was found that two genes were annotated as ubiquitin-conjugating enzyme E2, the highly expressed one was up regulated in 'Huangbao'. The two genes encoding for glutathione S-transferase showed the same trend ( Table 8 in Supplementary Files).
DEGs related to substance transport Nine genes encoding for the members of ATP binding cassette (ABC) family, except one up regulated gene, the rest of them were down regulated in 'Huangbao'. Gene expression of ion transporter in 'Huangbao' showed down regulation. Water and sugar transporter genes were up or down regulated in 'Huangbao' (Tables 9-12 in Supplementary Files).

Mesurement of carotenoid content
Carotenoid levels in the pericarp of ripe fruits from 'Huangbao' and 'Guanxi' were identi ed using high performance liquid chromatography (HPLC). In 'Huangbao' pummelo, lycopene was the main contributor for carotenoid accumulation, lutein and zeaxanthin were not indenti ed. However, minor amounts of carotenoids were detected in 'Guanxi' pummelo (Table 13 in Supplementary Files). HPLC chromatograms and standard curves were displayed in Additional le 8.

Validation of DEGs by quantitative real time PCR (qRT-PCR)
Eight genes were randomly selected and used to verify the stabilization of the transcriptome analysis results. qRT-PCR assay revealed that the relative expressions of the tested genes were coincident with the transcriptome pro les (Fig. 5). Primer and gene sequneces were listed in Additional le 9.

Discussion
Carotenoid precursor metabolisms in uence carotenoid accumulation Glyceraldehyde-3-phosphate and pyruvate are the precurcors for IPP and DMAPP synthesized through MEP pathway [5]. Fructose-1,6-phosphate conversion to glyceraldehyde-3-phosphate is catalyzed by fructose-bisphosphate aldolase. Enolase catalyzes the formation of phosphoenolpyruvate, and then, pyruvate kinase converts phosphoenolpyruvate to pyruvate. The expression levels of those genes encoding for the three enzymes were signi cantly higher in 'Huangbao' than that in 'Guanxi'. Those results indicated that the biosyntheses of the original precursors for carotenoids may be promoted in red pericarp of ripe fruits from 'Huangbao'.
The transcript levels of the genes annotated as PSY and CRTISO were distinctly lower in 'Huangbao' than that of in 'Guanxi'. It was investigated that the transcript abundance of PSY was not always coincident with the lycopene accumulation in tomato cultivar, its expression level at fruit ripening stage decreased 4 folds when compared with the colour break stage [26]. It has been reported that amino acid sequences in uenced the activities of PSY and CRTISO [29][30][31], which suggests that translation regulation should not be ignore. It has been found that down regulation of β-LCY and capsanthin/capsorubin synthase gene involved in the downstream of carotenoid biosynthetic pathway contributed to lycopene accumulation in the pulp of 'Hong Anliu' sweet orange fruits [32]. Taken together, it seemed that expression levels of genes encoding for PSY and CRTISO had no direct impact on lycopene accumulation in citrus fruits.

Candidate genes involved in regulating carotenoid metabolism
Genes participated in phytohormone metabolisms and hormone signaling pathway Carotenoid accumulation in tomao fruits during ripening stage was accompanied by the decreases of GA contents [33]. The application of exogenous GA led to increases in chlorophyll contents and decreases in carotenoid levels in the pericarp of lemon fruits [34]. GA are not likely to facilitate carotenoid accumulation in fruits. In this study, GA2OX and GRP were specially expressed in 'Huangbao' and 'Guanxi', respectively. Active GA (GA1 and GA4) were decreased by GA2OX reaction [35]. Overexpression of GRP caused decrease in transcript abundance of a gene encoding for the rst-step enzyme in GA biosynthetic pathway, accordingly, the contents of GA4 and intermediate products of GA decreased [36].
The results of this study showed that the genes involved in regulating active GA contents were different between the two pummelo cultivars.
Exogenous cytokinin was capable of increasing both contents of chlorophyll and carotenoids [37][38][39]. Inhibition of expression of cytokinin dehydrogenase (CKX) gene was accompanied by the increase in endogenous cytokinin content [40]. A gene annotated as CKX was down regulated in 'Huangbao' pummelo. This result indicated that the rate of cytokinin degradation in the pericarp of 'Huangbao' ripe fruits might be lower when compared to 'Guanxi' and consequently improve carotenoid accumulation in part.
Ethylene palys a key role in regulating fruit ripening process. It has been observed that improvement of nutritional quality and colour development in fruits were in uenced by ethylene [41]. Less ethylene production led to less carotenoid content in tomato fruits, correspondently, the orange fruit could not turn into red even at ripening stage [42]. The application of 1-methylcyclopropene (an inhibitor for ethylene biosynthesis) delayed colour development in citrus fruits [43]. By contrast, exogenous ethylene application elevated carotenoid contents in loquat fruits and accelerated fruit colour change [44]. ACO catalyzes the nal step of ethylene biosynthetic pathway. ACO3 was down regulated by virus induced gene silencing in hot pepper (Capsicum frutescens) fruits, resulting in inhibition of fruit colouration, which reveals the relationship between ethylene metabolism and fruit colour development [45]. ACS catalyzes the formation of 1-aminocyclopropane-1-carboxylate (ACC) [46]. ACO catalyzes the conversion of ACC to ethylene, and acts as a positive regulator in ethylene biosynthesis [47]. In the present study, although the expression level of ACS was lower in 'Huangbao' than that of in 'Guanxi', ACO was highly transcripted in 'Huangbao', the FPKM value was up to 1458.83. Those results indicated that ethylene is synthesized in the pericarp of fruits from both of the two cultivars. Ten DEGs were identi ed in ethylene signaling pathway, one of them was EBF1, and the rest of them were enthylene responsive transcription factors, including members in ERF and RAP2 subfamilies and EIN3. They were all lowly expressed in 'Huangbao'when compared with 'Guanxi', in particular, EBF1 and EIN3 were only detected in 'Guanxi'. ERF and RAP2 subfamilies belong to the family of ethylene responsive element binding protein and contain many members. For example, a total of 102 ERFs and 18 AP2s were identi ed in 'Newhall' sweet orange [48]. In tomato, RNAi repression of SlERF6 resulted in increases in carotenoid and ethylene contents [49]. The RNAi approach was also used to reduce SlAP2a expression level, the results showed that ethylene was over producted but carotenoid levels were decreased in tomato fruits [50]. Those reports indicate that ERF and AP2 both are negative regulators in ethylene signaling pathway, but play opposite roles in modulating carotenoid accumulation in fruits. EIN3 is a critical transcription factor in ethylene signal transduction pathway [51]. In addition, it is also involved in regulation of carotenoid metabolism. In papaya fruits, elevated transcript abundance of CpEIN3a accompanied with increased carotenoid content and up regulated genes in carotenoid biosynthetic pathway (CpPDS2/4, CpZDS, CpLCY-e and CpCHY-b) [52]. SlEIN2 silenced tomato fruits could not turn into red at ripening stage due to the inhibition in lycopene biosynthesis [53]. EBF expression is in uenced by feedback loop of ethylene signaling, and then involved in ne-tuning of ethylene signal transduction. EIN3 transcription factor could interact with EBF2 promoter and directly elelvated its expression level [51]. EBF belongs to F-box protein family and acts as a component in SCF complex, this complex could degrade EIN3 protein [54]. Recently, it has been reported that SlEBF3 could interact with all known ethylene insensitive-like proteins (EIL) in tomato fruits and overexpression of SlEBF3 reduced the total SlEIL protein levels [55]. Consequently, EBF protein in uences the sensitiveness of plant to ethylene. Taken together, the results obtained in the present study suggested that the high expression levels of negative transcription factors in ethylene signal pathway were capable of inhibiting ethylene signaling, nally, at least in part, led to carotenoids were weakly synthesized in the pericarp of ripe fruits from 'Guanxi'.
Application of indole-3-acetic acid (IAA) decreased transcript levels of several genes (PSY1, PSY3, PDS, ZISO and CRTISO) involved in the up-stream carotenoid metabolic pathway, but induced expression of down-stream genes (β-LCY1 and CRTR-β1). Consequently, tomato fruit could not turn into red due to the signi cant decrease in lycopene content and remarkable increase in violaxanthin and neoxanthin contents [56]. GH3 playes an important role in auxin homeostasis by converting active auxin into inactive IAA-amino acid conjugate [57]. A sum of 24 members in GH3 family were indenti ed in tomato fruits, in SlGH3-2 silenced fruits, the peak of ethylene production was delayed, nally, lycopene content decreased remarkably [58]. It can be observed that GH3 in uences fruit carotenoid content through by modulating auxin and ethylene contents. Genes annotated as GH3 and auxin responsive protein were detected in 'Huangbao' and 'Guanxi', all of them were down regulated in 'Huangbao'. It has been reported that ABC transporter B family member 19 (ABCB19) is capable of transporting auxin in plant [59]. In this study, ABCB19 wasonly detected in 'Guanxi'. Taken together, although GH3 was highly expressed in the pericarp of 'Guanxi' ripe fruits, the carotenoid contents were signi cantly lower in 'Guanxi' than that in 'Huangbao'. However, low expression levels of genes encoding for Auxin responsive protein and ABCB19 were associated with high carotenoid contents in 'Huangbao'. Auxin responsive protein and ABCB19 seem to act as negative regulators of carotenoid accumulation in fruits.

Transcription factors
The two opposite roles of DELLA protein in regulating GA metabolic pathway were observed in different plant species. In pea (Pisum sativum L.), DELLA protein down regulated the expression of PsGA2ox1 (a GA deactivation gene) [60]. The loss-of-function of DELLA protein induced dwarf phenotype in tomato was demonstrated by CRISPR/Cas9 approach [61]. On the contrary, a functionless SLN1 (a DELLA protein) protein could activate GA signaling pathway in the barley aleurone [62]. In this study, GA2OX and DELLA protein genes were specially expressed in the pericarp of ripe fruits from 'Huangbao' and 'Guanxi', respectively, which indicated that the GA activity might be inhibited in 'Huangbao'. Due to the effect of GA inhibition on carotenoid accumulation, GA deactivation should promote pericarp colouration in 'Huangbao'. It is still unclear whether DELLA protein gene is involved in regulating carotenoid accumulation.
NAC family genes in association with fruit ripening. In banana (Musa paradisiaca), MaNAC1/2 could physically interacte with MaEIL5, a downstream component of ethylene signaling pathway [63]. It has been found that the SlNAC1 gene negatively regulated the carotenoid and ethylene biosyntheses in tomato fruits (Solanum lycopersium). Overexpression of this gene was accompanied with decrease in transcript levels of SlPSY1, SlACO1 and SlACS2 and increase in expression levels of SlLCYb, SlLCYe and SlCYCB. which consequently decreased the total carotenoid and ethylene contents in fruits [64]. The in uences of SlNAC1 in carotenoid and ethylene metabolisms could be further demonstrated by antisence suppression of it, as a consequence, lycopene and ethylene contents were much higher in the transgenetic fruits than that in the wild type fruits [65]. Another member of NAC family, SlNAC4 acted as a positive regulator in upstream of the ripening regulator RIN [66].When SlNAC4 was silenced by virusinduced gene silencing, fruit ripening was delayed and ethylene production and lycopene accumulation were reduced in tomato fruits [67]. The previous research reveals that NAC family genes play different roles in regulating ethylene and carotenoid metabolic pathways. Four genes encoding for NAC domaincontaining protein were indenti ed in this study, while three of them were up regulated in 'Huangbao' , the last one were down regulated. Further assays need to be carried out to de ne clear roles of those NAC family genes in modulating carotenoid accumulation in pummelo fruits.
Just like NAC family genes, the members of MADS-box gene family act as a positive regulator or a repressor in regulating carotenoid and ethylene metabolisms. Ripening inhibitor (RIN) is a MADS-box transcription factor. It is involved in regulating carotenoid and ethylene metabolisms in both climacteric and non-climacteric fruits [68,69]. MADS-box proteins encoded by SlMBP15, FUL1 and FUL2 can interact with RIN, when those three genes were silenced individually, the carotenoid and ethylene contents as well as the expression levels of ethylene biosynthetic genes decreased in tomato fruits [70,71]. CsMADS6 directly bound to the promoters of PSY, PDS, LCYB1 and CCD1 and trans-activated their expressions, nally, caused the increase in carotenoid content in sweet orange calli [72]. On the other hand, some other members of MADS-box gene family play a reverse role in modulating carotenoid and ethylene accumulation. For example, when SlMBP8 or SlMADS1 was silenced, genes linked to ethylene (ACO, ACS, ERF1, E4 and E8) and carotenoid pathway (PSY1, PDS and ZDS) were up regulated, which was accompanied by improvement in carotenoid and ethylene biosyntheses [73,74]. Three MADS-box transcription factors were identi ed in this study, and all of them were specially expressed in 'Guanxi'.
WRKY transcription factor family is involved in response to stresses in plants [75]. PqWRKY1 acted as a positive regulator in triterpene ginsenoside biosynthetic metabolism in Panax quinquefolius and improved insensitiveness of transgenic Arabidopsis plants to abiotic stesses [76]. It was found that OfWRKY3 could bind to the promoter region of carotenoid cleavage dioxygenase 4 (OfCCD4) to transcriptionally activate its expression in sweet osmanthus (Osmanthus fragrans) [77]. CCD4 catalyzes the cleavage of carotenoid substrates [78]. When InCCD4 was silenced by CRISPR/Cas9 system, the total carotenoid content in petals of Ipomoea nil was increased 20-fold in comparison with non-transgenic plants [79]. Those previous reports suggest that WRKY transcription factor seems to facilitate isoprene consumption in terpene biosynthesis and accelerate carotenoid degradation, which lead to weakened accumulation of carotenoids. The present study showed that WRKY was signi cantly down regulated in 'Huangbao' compared with 'Guanxi'.
MYB transcription factor family contains a large number of members. For instance, a total of 1986 MYB and MYB-related putative proteins were identi ed in Gossypium hirsutum genome [80]. Some of MYB members participate in regulating anthocyanin and avonoid metabolisms [81,82], while some others are involoved in modulating carotenoid metabolism. Reduced expression of Carotenoid Pigmentation 1 (an R2R3-MYB transcription factor) was associated with improved carotenoid biosynthesis in Mimulus lewisii [83]. Another R2R3-MYB transcriptional factor (CrMYB68) in a stay-green mutant of Citrus reticulata cv Suavissima directly and negatively regulated CrBCH2 and 9-cis-epoxycarotenoid dioxygenase 5, as a consequence, the transformation of α-carotene and β-cartoene and the ABA biosynthesis were delayed [84]. In this study, one MYB44 transcription factor and two MYB related proteins were distinctly down regulated in 'Huangbao' in comparison with 'Guanxi'. It has been demonstrated that MYB44 mainly participates in stress response [85], but there isn't any information con rms that it is also involved in regulating carotenoid metabolism. Therefore, it remains unclear whether the three MYB related genes identi ed in this study are associated with colour development in fruit pericarp.

Genes associated with post translational modi cation
SlHDT3 is a histone deacetylase and belongs to tomato HD2 family, when SlHDT3 was silenced by RNAi technology, the amount of ethylene production and carotenoid accumulation decreased [86]. SlHDA3, a member of tomato RPD3/HDA1 family of histone deacetylases, played a negative role in regulating carotenoid accumulation in tomato fruits [87]. It can be ovbserved that members from different histone deacetylase families may paly contradictory roles in modulating carotenoid metabolic pathway. In this study, expression level of the gene encoding for histone acetyltransferase was 36.95-fold higher in 'Huangbao' than that of in 'Guanxi'. The present study results likely support the idea that histone acetylation promotes carotenoid accumulation in fruits.
EOLs can bind to ACS protein and degrade it, however, this effect can be reversed by 14-3-3 protein which interacts with EOLs and destabilizes them. 14-3-3 protein can also directly bind to ACS protein and maintain stabilization of ACS [88]. Transcript levels of E2 ubiquitin-conjugating enzyme genes were lower in yellow areas compared with orange areas in tomato fruits (Solanum lycopersicum cv Ailsa Craig) [70]. E3 ubiquitin ligase and E2 ubiquitin-conjugating enzyme both are involved in regulating fruit ripening. In this study, in 'Huangbao' pummelo, one 14-3-3-like protein gene was down regulated, additionally, one of the two genes annotated as E2 ubiquitin-conjugating enzyme was up regulated and another one was down regulated. The up regulated E2 ubiquitin-conjugating enzyme gene showed high transcript level in both of the two pummelo cultivars when compared to the down regulated one. Collectively, the up regulated gene encoding for E2 ubiquitin-conjugating enzyme may act as a positive regulator in modulating carotenoid metabolism.
Other metabolisms may provide a metabolic background for carotenoid accumulation Expression levels of DEGs associated with fatty acid metabolism and lipid transport were signi cantly lower in 'Huangbao' than that of in 'Guanxi'. It has also been investigated that the some fatty acid biosynthetic genes were down regulated in fruits from orange-pericarp mutant in pummelo (Citrus grandis (L.) Osbeck) [89]. Transcript abundance of genes annotated as sugar, water and metal element transporters showed signi cant changes. Some signi cantly up regulated genes encoding for ribosomal proteins, mRNA elongation factors and cytochrome c oxidase were observed in 'Huangbao', which suggested that the total metabolic activity was higher in the pericarp of 'Huangbao' ripe fruits compared with 'Guanxi'. Those DEGs involved in other metabolic pathways might provide a metabolic background for lycopene accumulation in the pericarp of ripe fruits from 'Huangbao'.

Conclusions
Large amount of lycopene accumulated in the pericarp of ripe fruits from 'Huangbao' pummelo was associated with signi cantly higher expression levels of genes encoding for fructose-bisphosphate aldolase, enolase and pyruvate kinase isozyme, however, transcript levels of genes encoding for phytoene synthase and prolycopene isomerase, those are two key enzymes in lycopene biosynthesis pathway, were signi cantly lower in 'Huangbao' than that of in 'Guanxi'. Fructose-bisphosphate aldolase, enolase and pyruvate kinase isozyme catalyze the formation of 3-phosphate-glyceraldehyde and pyruvic acid. The two productions are precursors for IPP and DMAPP biosyntheses in MEP pathway, subsequently, IPP and DMAPP act as precursors for carotenoid biosynthesis. Those results indicate that the regulation in carotenoid metabolism is expanded to the carbon skeleton biosynthetic pathway in the pericarp of pummelo fruits. Further experiments need to be carried out to verify whether the activities of phytoene synthase and prolycopene isomerase are regulated by post-translational modi cation. Gibberellin-2-βdioxygenase gene was up regulated in 'Huangbao', however, genes related to auxin response and auxin transport, as well as several negative regulators in ethylene signaling pathway were down regulated in 'Huangbao'. It has been reported that phytohormone is involved in regulating carotenoid metabolism. Therefore, the DEGs related to phytohormone metabolism or its singal transduction may participate in modulating carotenoid accumulation in fruits. Additionally, transcription factors (DELLA, NAC, MADS-box, WRKY) and genes encoding for post translational modi cation proteins (histone acetyltransferase and E2 ubiquitin-conjugating enzyme) involved in regulating ethylene and carotenoid metabolisms were identi ed in this study. Those ndings provide a promising new research direction for better understanding regulation of carotenoid accumulation in fruits.

Methods
Fruit samples 'Huangbao' was selected from bud mutation of 'Guanxi' pummelo (Citrus grandis (L.) Osbeck) for its red colour pericarp. This phenotype is observed during fruit ripening stage, which is distinctly different with the yellow pericarp of wild type fruits (Fig. 6). 'Huangbao' and 'Guanxi' pummelo cultivars were both grafted on 'Sour pummelo' (Citrus grandis) and planted at 3.0 m x 3.5 m in a commercial orchard located in Zhangzhou City, China. Ripe fruit samples of each cultivar were harvested at the outside canopy of three 4-year-old trees with same growth vigor. Three fruits were sampled from each tree and each tree was taken as one biological replicate. Pericarp of three fruits harvested from each tree was cut into small pieces, mixed together, immediately frozen in liquid nitrogen and then stored at - Subsequently, the selected cDNAs were ampli ed by PCR using phusion high-delity DNA polymerase, universal PCR primers and index (X) primer. Finally, cDNA libraries were puri ed with AMPure XP system and quanti ed using the Agilent 2100 bioanalyzer. The cDNA libraries were sequenced on an Illumima Hiseq 4000 platform.
Sequence quality control and ungiene assembly Raw data (raw reads) of fastq format were rstly processed through in-house perl scripts. In this step, low quality reads and reads containing adapter or ploy-N were removed, and the remaining reads were taken as clean reads. Q30 and GC-content of the clean data were calculated. All the following analyses were based on clean reads with high quality. The clean reads were deposited in the National Center for Biotechnology Information (NCBI) short Read Archive database (https://www.ncbi.nlm.nih.gov/sra) with accession number of SRP131839. All clean reads were used for de novo assembly using Trinity software  Dalian, China). The assay mixture was in a nal volume of 25 µL, which consisted of 12.5 µL of SYBR Premix DimerEraser (2×), 1 µL of each primer (10 µM), 1 µL of cDNA template, and 9.5 µL of RNase-free water. Cycling protocol consisted of an initial step of 95 °C for 30 s, then followed by 40 cycles of 95 °C for 5 s, 55 °C for 30 s and 72 °C for 30 s, and completed with a melting curve analysis. The speci city of qRT-PCR primers was con rmed by melting curve. The primer sequences used in this study were designed using primer3 (http://bioinfo.ut.ee/primer3-0.4.0/primer3) and listed in Additional le 9. The fold change of gene expression was calculated by using

Availability of data and materials
The clean reads are available in the NCBI short Read Archive repository (https://www.ncbi.nlm.nih.gov/sra), available from accession SRP131839. All supporting data were included within the manuscript and its additional els. Xinkun Lu obtained the permission of fruit sample collection from the orchard owner and undertook the formal identi cation of the samples and provide details of any voucher specimens deposited.