Transcriptome changes during fruit development and ripening of sweet orange (Citrus sinensis)

Plant material and RNA preparation

WT and MT plants were both cultivated in the same orchard at the Institute of Citrus Research (Guilin, Guangxi Province, China), with the same climatic conditions. Fruit samples were harvested at 120, 150, 190, and 220 DAF (days after flowering) from three different trees in 2009. At each developmental stage, ten representative fruits were sampled from each tree. The pulp was separated from the peel, and the pulp was sliced. The sliced WT pulps samples were combined with one another (similarly for the MT samples), snap-frozen in liquid nitrogen and kept at -80°C until required [24, 25]. One aliquot was used to extract RNA isolation, as described previously [26]. The remainder of the powder was used for the determination of sugar and organic acid composition and concentration, and the content of H₂O₂.

RNA-seq and functional assignment

The WT and MT fruit pulp harvested at 120, 150, 190, and 220 DAF was subjected to RNA-seq using an Illumina Genome Analyzer at Beijing Genomics Institute (Shenzhen) in 2009. The abundance of each tag was normalized to one transcript per million (TPM) for between sample comparison purposes. The raw data was filtered to remove low quality sequences including ambiguous nucleotides, adaptor sequences, and below 3 TPM, as described previously [25]. The sequencing data can be accessed at the website: http://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?token=dxqjxoygumyauzm&acc=GSE22505. To link the expressed signatures to known genes from orange, the TIGR unigene dataset (http://compbio.dfci.harvard.edu/tgi/cgi-bin/tgi/gimain.pl?gudb=orange) was used as a reference database. The Z-score method using the p-value as a statistical significance index [27] was applied to identify differentially expressed genes. A cluster analysis was performed according to Eisen et al. [28]; the log₂ of TPM for each gene was used for the hierarchical clustering analysis. Gene Ontology (GO) categorization was carried out as described previously [25]. The ultra-geometric test was applied to perform GO enrichment analysis. In the significance analysis of the enrichment of a GO item, the p-value represents the probability of satisfying the hypothesis that the designated genes involved in the GO item has not been enriched (statistical significance at P = 0.05).

Real-time quantitative RT-PCR

The differential expression of a selection of the genes identified as being differentially expressed was validated by applying real-time quantitative RT-PCR (qRT-PCR). The sequences of the primer pairs (designed using Primer Express 3.0 (Applied Biosystems, Foster City, CA, USA)) are listed in additional file 1. All qRT-PCRs were performed using an ABI 7500 Real Time System (Applied Biosystems) using the actin gene as the reference [20]. Primers for both the target gene and the reference were diluted in SYBR GREEN PCR Master Mix (Applied Biosystems) and 20 μL of the reaction mix were added to each well. Reactions were performed via an initial incubation at 50°C for 2 min and at 95°C for 10 min, and then cycled at 95°C for 15 s, and 60°C for 60 s for 40 cycles. The resulting data were handled by the instrument on-board software Sequence Detector Version 1.3.1 (Applied Biosystems).

Analysis of sugar, organic acid and H₂O₂

Soluble sugar and organic acid composition and concentrations were determined by gas chromatography (GC) using 3 g of the powdered pulp as described previously [29] with minor modifications. The powder was suspended in chilled 80% methanol and then held in a 75°C water bath for 30 min. After a 2 h ultrasonic extraction and centrifugation at 4000g for 10 min, the supernatant was collected and 1 mL internal standard (2.5% w/v phenyl-β-D-glucopyranoside, 2.5% w/v methyl-α-D-glucopyranoside) was added. The solution was made up to 50 mL with 80% methanol, and a 2 mL aliquot was centrifuged at 12000g for 15 min. A 0.5 mL aliquot of this final supernatant was vacuum-dried and then re-dissolved in 800 μL 2% w/v hydroxylamine hydrochloride in pyridine at 75°C for 1 h. Then 400 μL hexamethyldisilazane and 200 μL trimethylchlorosilane were added and the sample was held at 75°C for 2 h. A 0.5 μL aliquot was used for GC analysis in an Agilent 6890N device (Santa Clara, CA, USA) equipped with a flame ionization detector. A capillary column (HP-5, 5% phenyl-methylpolysiloxane, 30 m × 25 μm i.d. × 0.1 μm) was employed, with nitrogen as the carrier gas at a flow rate of 45 mL/min, and flow-rates of hydrogen and air set to 40 mL/min and 450 mL/min, respectively. Sugars and organic acids were identified through a comparison of retention times using standard compounds from Sigma (St. Louis, MO, USA). The concentration of H₂O₂ was measured using a hydrogen peroxide detection kit supplied by Nanjing Jiancheng Institute of Biological Technology (Nanjing, China). A 0.8 g sample of powdered pulp was suspended in 7.2 ml saline (0.90% w/v of NaCl) and centrifuged for 10 min at 10, 000g. The intensity of yellow complex formed by the reaction of molybdate and H₂O₂, as measured spectrophotometrically at 405 nm, was used to assess the concentration of H₂O₂. Three replicates were conducted for each sample.

The fruit transcriptome sampled at four developmental stages

Transcriptome changes during fruit development and ripening

To map tags to known genes, a reference citrus unigene dataset containing 26,826 contigs and 73,607 singletons was used. The procedure identified between 68.1% and 76.2% of the tags (additional file 4), of which 20,155 to 36,173 (31.8% to 33.1%) produced unambiguous identifications (one tag mapping to one gene). The libraries were relatively uniform with respect to mapping efficiency. A total of 18,829 genes were detected in at least one of the four stages in the wild type sweet orange, of which 8,825 genes were expressed in all the four stages. In this study, we solely used the wild type sweet orange as a model to demonstrate the transcriptome changes during fruit development and ripening. Three genes were most highly expressed in wild type, two of which were encoded a stress-response protein (one 22 kDa polypeptide response to low temperature stress) and a heat shock protein, while the function of the third one (F28C11.8) is unknown. Changes in the transcriptome during fruit development and ripening were examined by cluster analysis of gene expression patterns, which arranged the 18,829 genes into 22 groups (Figure 1); the 10,005 (53.1%) genes expressed in three or less of the three stages fell into groups 1 to 11. The largest group (20) comprised the 3,075 (16.3%) genes whose expression increased continuously during fruit development and ripening; this group included the genes encoding sucrose phosphate synthase, cysteine proteinase and a sucrose transporter. The second largest group (10) contained the 2,970 (15.8%) genes which were not expressed at 120 DAF but maintained a stable expression level at other three developmental stages. The 2,618 (13.9%) genes in group 2 were not expressed at 120 and 190 DAF. The cluster analysis also revealed that the abundance of 89.7% of the transcripts detected in the WT pulp varied over the course of fruit development and ripening (Figure 1). Many of the transcripts were single stage-specific (additional file 5). A comparison of expression patterns between WT (Figure 1) and MT (additional file 6) revealed that 20 of the groups were common to both, while 97.5% of the genes expressed in MT showed a similar expression pattern to that in WT (the exceptions belonged to groups 13 and 19).

Differentially expressed genes during fruit development and ripening

Of the 18,829 genes detected in the WT, 9,377, 7,886, and 7,757 were differentially expressed between 120 DAF and 150 DAF, 150 DAF and 190 DAF, and 190 DAF and 220 DAF, respectively. Of these, 36.7% were assigned to one of 18 GO categories (Figure 2). The categories "metabolic process", "cellular process", "establishment of localization", "localization", "biological regulation", "pigmentation", and "response to stimulus" based on biological process captured most of these genes. Several of the differentially expressed genes were associated with cell wall metabolism and softening (Table 2). A gene encoding pectinesterase (TC14614) showed a decrease in expression during fruit development and ripening, consistent with the modification of pectin which occurs during ripening. A xyloglucan endotransglycosylase gene (TC9277) was up-regulated during fruit development and ripening, while another (TC25537) was only detectable at the early stage of fruit development. The abundance of three expansin (EY741042, TC20787, and EY738078) and three alpha-expansin (EY727139, EY743186, and EY699153) transcripts increased during fruit development and ripening. Some genes involved in sucrose metabolism, the TCA cycle and carotenoid biosynthesis were found to be differentially expressed (Table 2). The expression of a sucrose phosphate synthase gene family member (EY677217) increased during ripening, and its expression profile correlated well with the rise in total soluble solids. Related genes encoded CS (EY703799), aconitase-iron regulated protein 1 (CB293814), phytoene synthase (PSY) (EY722043) and phytoene desaturase (PDS) (TC3). Some of the differentially expressed genes were associated with stress response (Table 2): these included genes encoding a heat shock protein and various metallothioneins. Most of the metallothionein genes were up-regulated in the later stages of ripening. Three cysteine protease genes (EY718215, TC22740, and EY707198), were also highly expressed at the late stage of fruit development and ripening. Also, the expression of a valencene synthase gene (TC363) increased during fruit development and ripening.

Developmental difference between fruit transcriptome of MT and WT

The comparison between the transcriptomes of WT and MT identified 6,540 genes showing a two-fold or greater difference in expression level at 120 DAF, 3,529 at 150 DAF, 4,601 at 190 DAF, and 3,289 at 220 DAF. Application of the Z-score method [27] suggested that the expression of 634, 568, 540, and 616 of these was significantly different at p < 0.05 and |log₂Ratio| ≥ 1 in the four developmental stages, respectively (additional file 7). The top ten differentially expressed genes in each of the four developmental stages have been listed in additional file 8. Many encode stress-related products, such as cysteine protease Cp5 and metallothionein-like protein. At all the four developmental stages the number of up-regulated genes was less than that of down-regulated genes (Figure 3).

When the expression profiles of the genes differentially expressed between MT and WT were subjected to a cluster analysis (additional file 9), over one half (492/883) turned out to be up-regulated in MT at all the developmental stages except 150 DAF. Some examples of this large cluster included genes encoding Cu/Zn superoxide dismutase (TC12069), ascorbate peroxidase (EY685405, TC22775 and TC14669), MYB1 (EY649968) and plastidic glucose 6-phosphate/phosphate translocator (EY722703). Only five genes were detected as differentially expressed at all four stages; one encoded a cysteine protease (Cp5) (TC5370), which had been reported exhibiting the acid-activatable cysteine protease forms [30], three had no assigned functions, and one shared no homology with any entry in. The five genes could be classified into two groups (two genes in one group, three genes in the other), with their expression patterns of being opposite to each other but with both groups having a turning point at 150 DAF (additional file 10).

Verification of differentially expressed genes during fruit development and ripening

RNA sampled from fruit of MT and WT at the four selected stages of fruit development and ripening provided the template for qRT-PCR validation of the sequence-based transcription profiles for 22 of the differentially expressed genes (Figure 4). Linear regression [(RNA-seq value) = a(RT-PCR value) + b] analysis showed an overall correlation coefficient of 0.8379, indicating a good correlation between transcript abundance assayed by real-time PCR and the transcription profile revealed by RNA-seq data (Figure 4E).

Changes in fruit soluble sugars, organic acids, carotenoids and H₂O₂ content

Since the expression of a number of genes implicated in carbohydrate metabolism and mitochondria-related citric acid metabolism varied greatly during fruit development and ripening, an attempt was made to monitor the dynamics of pulp soluble sugar and organic acid content. The content of soluble sugars increased markedly during the late stages of fruit development and ripening in both WT and MT (Figure 5A). The concentrations of glucose and fructose remained rather constant up to 150 DAF, and thereafter tended to be higher in MT than in WT. The concentration of sucrose was higher in MT than in WT throughout fruit development and ripening. The citric acid content fluctuated, but was lower in MT than in WT at all developmental stages (Figure 5B). The concentration of quinic acid decreased substantially over the course of fruit development and ripening, but did not differ between WT and MT. Carotenoids and lycopene both accumulated over time in MT, but remained at a low level in WT (Figure 5C). Finally, H₂O₂ content fell as the fruit developed and ripened, but was higher in MT than in WT at 120 DAF (Figure 5D).

Transcriptome dynamics during fruit development and ripening

Fruit ripening is a highly coordinated, genetically programmed and irreversible process which involves a series of physiological, biochemical, and organoleptic changes allowing for the development of an edible ripe fruit [31]. Fruit ripening in citrus is accompanied by carbohydrate build-up, acid reduction, carotenoid accumulation and chlorophyll degradation [32]. The bud mutant MT produces fruit with high sucrose and lycopene, but low citric acid content [20]. Here, transcriptome changes over the course of fruit development and ripening in MT and WT were monitored and annotated. We must point out that we collected samples from three different trees and pooled for RNA-seq analysis for each developmental stage. We did not include biological replicates considering that we used a pair of genotypes (WT and MT) and the data from WT and MT could corroborate to each other to some extent. We sequenced each pool once technically since the next generation sequencing data are highly replicable with relatively little technical variation [33]. We also used real-time qRT-PCR to verify the transcription profile revealed by RNA-seq data (Figure 4), with an overall correlation coefficient of 0.8379 which indicated the RNA-seq data was reliable. The RNA-seq approach detected a similar number of genes (18,829 genes in WT and 19,440 genes in MT) in both genotypes, with no apparent difference between the global number of genes expressed in the two genotypes at any of the four stages of fruit development and ripening sampled (additional file 4). The WT and MT patterns of gene expression were also largely alike (20 of the 24 expression patterns were similar). Thus, we believed that overall picture of the transcriptome captured by RNA-seq is robust.

An overwhelming proportion of the genes identified (84.8% in MT and 89.7% in WT) varied in their level of expression over the course of fruit development and ripening, reflecting the occurrence of a massive genetic re-programming. A large number of these genes were expressed in a stage-specific manner, which implicates their involvement in physiological processes which take place only at a specific developmental stage(s). A major group of the differentially expressed genes was involved in cell wall modification, which is not surprising since the major textural changes associated with the softening of fruit are due to enzyme-mediated alterations in the structure and composition of the cell wall [34], especially the cell walls of juice sacs in citrus. Changes in the activity of several cell wall-related genes were known to result in the abnormal development of juice sac granulation [35, 36], while modifications in cell wall structure or in the components of the membranes of the segments and juice sacs during fruit development and ripening clearly influenced the formation of the fruit pulp melting characteristic [37]. One of the cell wall-related genes revealed here encoded a pectinesterase, an enzyme which modifies the assembly and disassembly of pectin, a common component of the primary cell wall. In tomato, the gene for pectinesterase was highly expressed prior to ripening, and was down-regulated by ethylene as ripening begins [38]. Here, the expression of the gene encoding the pectinesterase PPE8B precursor decreased as the WT fruit matured (Figure 4A), while in the blood orange, the expression of a pectinesterase gene has been shown to increase during fruit development and ripening [39]. It might be due to the different members identified in the two studies. The hemicellulose xyloglucan is a common component of the cell wall, and is hydrolysed and transglycosylated by xyloglucan endotransglycosylase in growing tissues and ripening fruits [40]. Here, a gene encoding this enzyme was up-regulated during fruit development and ripening (Figure 4A), indicating its probable role in cell wall degradation during ripening. Both the early enlargement of the citrus fruit driven by cell expansion and the later ripening process require the presence of expansins to loosen the cell walls [41], and several genes encoding expansins were detected in the present study and their expression was higher at the later stages of fruit development and ripening.

The accumulation of carbohydrates represents one of the most obvious changes which occur during citrus fruit development and ripening. The perceived changes in expression of genes involved in carbohydrate metabolism here were consistent with the findings of other transcriptomic analyses in citrus [12]. The type of sugar deposited to a high level in the cell vacuole in citrus is predominantly sucrose, unlike in grape, where it is glucose and fructose [6]. Reflecting this difference, the expression profile of the gene encoding sucrose synthase in citrus [12] was rather different from that in the grape [6]. Nevertheless, in both species, sugar is important for the regulation of colour development, and perhaps also for other ripening processes [42]. It was notable that the gene encoding a sucrose phosphate synthase was induced during fruit development and ripening, perhaps because the activity of this enzyme may be required for the re-synthesis of sucrose to allow its further transport to the vacuole [43]. Of some interest also was the behaviour of a gene encoding valencene synthase, an enzyme the activity of which was known to be induced by ethylene, and which was part of the natural ripening process in citrus [11]. Valencene is an important component of the aroma of the ripe sweet orange fruit. Both in the present experiments (Figure 4A) and elsewhere [44], valencene synthase transcript accumulated in the ripening fruit.

Molecular processes involved in the formation of the fruit traits of red flesh sweet orange

Suppressive subtraction hybridization, in combination with cDNA microarray analysis, has identified a set of 267 genes which were differentially expressed between MT and WT [24]. RNA-seq technology is more capable of identifying a nearly complete inventory of transcripts and by this method a total of 582 genes were found to be differentially expressed between WT and MT at the stage of 170 DAF [25]. In the current study, we tested on four developmental stages of WT and MT and extended our understanding of the global and dynamic changes during fruit development and ripening in MT and WT. Almost all of the members of the 267 gene set revealed by SSH strategy were also identified in the present study, although some discrepant expression patterns were apparent between this new data set and previously assembled microarray-based set. For example, the abundance of the transcript encoding cysteine protease, which appeared to differ between MT and WT in both studies, was documented by the microarray analysis as being lower at all developmental stages, whereas it appeared to be higher at 150 DAF in MT in the present study (additional file 10).

The major biological processes occurring in the mitochondria (TCA cycle, coupling electron transfer, and oxidative phosphorylation) were remarkably altered in MT. The intermediates of the TCA cycle can be channelled into the syntheses of fats, terpenoids, porphyrins, nucleotides, and amino acids. In MT, the level of citric acid, the major organic acid present in citrus fruit, was consistently around 25% that present in WT fruit (Figure 5). However, no major difference was detected in the expression of the genes encoding CS and aconitase-iron regulated protein (Figure 4A), two predominant enzymes involved in the TCA cycle. Five differentially expressed genes, all associated with mitochondria-related processes (coupling electron transfer and oxidative phosphorylation) were down-regulated in MT compared with WT (Figure 4B). These included the genes encoding a NADH-ubiquinone oxidoreductase 75kDa subunit (TC7962), cytochrome C oxidase (TC18256), and cytochrome oxidase subunit 3 (TC156), suggesting that MT mitochrondria were capable of less efficient electron transport than that WT ones. If, as a result, flux through the TCA cycle is decreased, the accumulation of citric acid is likely to be compromised. In addition, PDS and ZDS catalyzing desaturation of phytoene to lycopene involve net electron transfer [45]. In tomato, a NAD(P)H dehydrogenase complex which participating in electron transfer was involved in carotenoid biosynthetic pathway [46], suggesting the possible exist of cross-talk between electron transfer and carotenoid accumulation in sweet orange.

In the plant cell, the mitochondrial electron transport chain is a major site of reactive oxygen species (ROS) production [47]. Here, the concentration of one of the primary ROS molecules (H₂O₂) was higher in MT than that in WT at 120 DAF (Figure 5D). The delicate balance between antioxidant defence and ROS production can be disrupted by either compromised antioxidant defence or the inhibition of electron flow [48]. Here, the primary anti-oxidant enzymes (SOD, APX and GR) were more active in MT than in WT pulp, suggesting that the level of oxidative stress may be greater in MT than in WT [49]. The expression of a large number of stress-related genes was also substantially different in MT and WT (additional file 8). Functioning as an important ROS scavenging pathway, the ascorbate-glutathione cycle requires the cofactor NADPH, which is provided by the OPP pathway [50]. This pathway is also a major source of NADPH for many biosynthetic processes, including carotenoid biosynthesis [51], in non-photosynthetic organs such as the fruit. The gene encoding glucose 6-phosphate dehydrogenase, which is considered as the first and rate-limiting enzyme of the OPP pathway in all cells [52], was up-regulated in MT. A statistical analysis of the qRT-PCR result confirmed that the level of transcription of this gene was significantly higher in MT than in WT, especially at 120 DAF (Figure 4C). Carotenoid biosynthesis, another effective anti-oxidative process [53], was also higher in MT than in WT. Lycopene is the most potent antioxidant among the carotenoids [54]. The expression level of several carotenoid biosynthesis genes, encoding namely PSY, ZDS, lycopene β-cyclase (LCYb) and CCS, was greatly changed in MT. PSY is generally accepted to be a rate-limiting enzyme in carotenoid biosynthesis pathway. The CCS product is an enzyme which is mechanistically similar to LCYb, and the low transcript level of CCS may well be responsible for the accumulation of lycopene in red grapefruits (Citrus paradisi) [19]. The qRT-PCR analysis confirmed that both upstream genes (PSY and ZDS) were up-regulated and both downstream ones (LCYb and CCS) down-regulated in MT (Figure 4D), consistent with the mechanism regulating lycopene accumulation in tomato [55].