Comparative transcriptome analyses of a late-maturing mandarin mutant and its original cultivar reveals gene expression profiling associated with citrus fruit maturation

Plant materials

The late-maturing mutant ‘Huawan Wuzishatangju’ (HWWZSTJ) (Citrus reticulata Blanco) and its original cultivar ‘Wuzishatangju’ (WZSTJ) were planted in the same orchard in South China Agricultural University (23°09′38″N, 113°21′13″E). Ten six-year-old trees of each cultivar were used in this experiment. Peels (including albedo and flavedo fractions) from fifteen uniform-sized fresh fruits were collected on the 275^th (color-break stage, i.e., peels turns from green to orange) and 320^th (maturing stage) days after flowering (DAF) of HWWZSTJ and 275^th (maturing stage) DAF of WZSTJ (Fig. S1) in 2012 and pools were named T3, T1 and T2, respectively. Peels from fifteen uniform-sized fresh fruits of HWWZSTJ and WZSTJ were collected on the 255^th, 265^th, 275^th, 285^th, 295^th, 305^th, 315^th and 320^th DAF in 2012 and used for expression analyses of candidate transcripts associated with chlorophyll, carotenoid biosynthesis and ABA metabolism. All samples were immediately frozen in liquid nitrogen and stored at −80 °C until use.

RNA extraction, library construction and RNA-Seq

Total RNA was extracted from peels according to the protocol of the RNAout kit (Tiandz, Beijing, China) and genomic DNA was removed by DNase I (TaKaRa, Dalian, China). RNA quality was analyzed by 1.0% agarose gel and its concentration was quantified by a NanoDrop ND1000 spectrophotometer (NanoDrop Technologies, Wilmington, DE, USA). RNA integrity number (RIN) values (>7.0) were assessed using an Agilent 2100 Bioanalyzer (Agilent Technologies, Santa Clara, CA, USA).

Construction of RNA-Seq libraries was performed by the Biomarker Biotechnology Corporation (Beijing, China). mRNA was enriched and purified with oligo (dT)-rich magnetic beads and then broken into short fragments. The cleaved RNA fragments were reversely transcribed to the first-strand cDNA using random hexamer primers. The second-strand cDNA was synthesized using RNase H and DNA polymerase I. The cDNA fragments were purified, end blunted, ‘A’ tailed, and adaptor ligated. The distribution sizes of the cDNA in the three libraries were monitored using an Agilent 2100 bioanalyzer. Finally, the three libraries were sequenced using an Illumina HiSeq™ 2500 platform.

Transcriptome assembly and annotation

Sequences obtained in this study were annotated in reference to the genome sequence of Citrus sinensis (Xu et al., 2013; Wang et al., 2014) using a TopHat program (Trapnell, Pachter & Salzberg, 2009). Functional annotation of the unigenes was performed using BLASTx (Altschul et al., 1997) and classified by Swiss-Prot (SWISS-PROT downloaded from European Bioinformatics Institute by Jan., 2013), Clusters of Orthologous Groups of Proteins Database (COG) (Tatusov et al., 2000), Kyoto Encyclopedia of Genes and Genomes Database (KEGG, release 58) (Kanehisa et al., 2004), non-redundant (nr) (Deng et al., 2006) and Gene Ontology (GO) (Harris et al., 2004). The number of mapped and filtered reads for each unigene was calculated and normalized giving the corresponding Reads Per Kilobases per Million reads (RPKM) values. DEGs between the two samples were determined according to a false discovery rate (FDR) threshold of <0.01, an absolute log2 fold change value of ≥1 and a P-value <0.01.

Gene validation and expression analysis

Data from RNA-Seq were validated using qPCR. All pigment-related (chlorophyll metabolism, carotenoid biosynthesis and ABA metabolism) uni-transcripts were selected to elucidate their expression patterns at all peel coloration stages of HWWZSTJ and WZSTJ with specific primers (Table S1). The citrus actin gene (accession No. GU911361.1) was used as an internal standard for the normalization of gene expression. Expression levels of all pigment-related uni-transcripts were determined using qPCR in an Applied Biosystems 7500 real-time PCR system (Applied Biosystems, CA, USA). A total of 20.0 µl reaction volume contained 10.0 µl THUNDERBIRD SYBR qPCR Mix (TOYOBO Co., Ltd.), 50×ROX Reference dye, 2.0 µl Primer Mix (5.0 µM), 6.0 µl ddH₂O, and 2.0 µl cDNA (40 ng). The qPCR parameters were: 94 °C for 60 s then 40 cycles of 95 °C for 15 s, 55 °C for 15 s and 72 °C for 30 s. All experiments were performed three times with three biological replicates. Relative expression levels of selected transcripts were calculated by the 2^−ΔΔCT method (Livak & Schmittgen, 2011).

All pigment-related genes (chlorophyll metabolism, carotenoid biosynthesis, ABA metabolism) were cloned using specific primers (Table S2). The 20.0 µl of reaction volume contained 2.0 µl 10×PCR buffer, 2.0 µl dNTP (2.0 mM), 0.2 µl of each primer (10 µM), 2.0 µl DNA (100 ng), 0.2 µl LA Taq and 13.4 µl ddH₂O. PCR reaction procedure was 94 °C for 4 min then 35 cycles of 94 °C for 30 s, 55 °C for 30 s and 72 °C for 2 min, with a final 72 °C for 10 min. Nucleotide sequences of the pigment-related genes were analyzed using the National Center for Biotechnology Information (NCBI) Blast program (http://www.ncbi.nlm.nih.gov/BLAST). ORFs were made using the NCBI ORF Finder (http://www.ncbi.nlm.nih.gov/gorf/gorf.html). Alignments were done using ClustalX 1.83 and DNAMan software. Phylogenetic analysis of deduced amino acid sequences were performed using MEGA (version 5.0) and the Neighborjoining method with 1,000 bootstrap replicates.

RNA-Seq analyses

To obtain differentially expressed genes (DEGs) between HWWZSTJ and WZSTJ, three libraries (T1, T2 and T3) were designed for RNA-Seq. As shown in Table 1, 26,403,257, 29,163,126, and 28,606,868 raw reads were obtained respectively from the three libraries. After removing low-quality bases and reads, a total of approximately 17.0 Gb clean reads were obtained. The GC contents for T1, T2 and T3 were 44.27%, 44.62% and 44.20%, respectively (Table 1). The range of most transcripts length was 100–200 bp (Fig. S2). Q30 percentage (percentage of sequences with sequencing error rate lower than 0.01%) for each sample was over 90% (Table 1).

A total of 44,664,047, 49,507,338 and 48,492,905 reads were mapped which accounted for 84.58%, 84.88% and 84.76% of the total reads, respectively (Table 2). Number of unique mapped reads accounted for 97.14% (T1), 97.25% (T2) and 97.19% (T3) of the total mapped reads compared with 2.86% (T1), 2.75% (T2) and 2.81% (T3) for multiple mapped reads, respectively. Those results suggested that the throughput and sequencing quality was high enough for further analyses.

Analyses of differentially expressed genes (DEGs)

DEGs were screened by comparison between any two of the three libraries using p < 0.01, FDR < 0.01 and Fold Change ≥ 2 as thresholds. A total of 2,687, 3,002 and 1,834 DEGs were obtained between the T1 and T3, T2 and T1, T2 and T3 libraries, respectively (Fig. 1A). Among those DEGs, 1,162, 1,567 and 770 were up-regulated and 1,525, 1,435 and 1,064 were down-regulated (Fig. 1B). Transcriptional levels of DEGs in HWWZSTJ on 320^th DAF were lower than that on 275^th DAF in HWWZSTJ suggesting that transcriptional levels of DEGs decreased during fruit maturation in HWWZSTJ (Fig. 1B).

Functional annotation of transcripts

A total of 299 new transcripts were annotated using five public databases (Nr, Swiss-Prot, KEGG, COG and GO). A summary of the annotations was shown in Table S3. Maximum number of annotation of differentially expressed transcripts (2,954) was in the Nr databases by comparison between T1 and T3, T2 and T1, T2 and T3, followed by GO databases (2,648) (Table S4). The differentially expressed transcripts were classified into three categories in GO assignments: cellular component, molecular function and biological process. DEGs between T1 and T3, T2 and T1, T2 and T3 were all significantly enriched in pigmentation, signaling and growth biological processes (Fig. S3A). Based on COG classifications, differentially expressed transcripts were divided into 25 different functional groups (Fig. S3B). DEGs between any two of the three libraries (T1-VS-T3, T2-VS-T1, T2-VS-T3) were assigned to 91, 100 and 91 KEGG pathways, respectively (File S1), and phenylalanine metabolism, porphyrin and chlorophyll metabolism, and flavonoid biosynthesis were the three significantly enriched biological processes (Table 3).

Verification of the accuracy of the RNA-Seq data using qPCR

Twelve DEGs with significant differences from the three libraries were selected for verification of RNA-Seq data by qPCR. Linear regression analysis showed an overall correlation coefficient of 0.828, indicating a good correlation between qPCR results and the transcripts per kilobase million from the RNA-Seq data (Fig. S4).

DEGs related to carotenoid biosynthesis, chlorophyll and ABA metabolism

Analyses of the expression data obtained through RNA-Seq revealed that most DEGs were involved in carotenoid biosynthesis, chlorophyll and ABA metabolism. The main transcripts involved in the three pathways were shown in Table 4 and heatmaps were made based on transcripts per kilobase million from the RNA-Seq data (Fig. 2). Three transcripts (Cs8g15480, Pheophorbide a oxygenase; Cs5g16830, Chlorophyllase type 0 and Cs3g19690, Chlorophyll synthase) involved in chlorophyll degradation, six transcripts (Cs3g19770, Delta-aminolevulinic acid dehydratase; Cs9g13460, Magnesium-chelatase subunit H; Cs2g05100, Magnesium-chelatase subunit ChlI-1; Cs7g19710, Magnesium-protoporphyrin O-methyltransferase; Cs6g16200, Magnesium-protoporphyrin IX monomethyl ester (oxidative) cyclase 1 and Cs1g06850, Protochlorophyllide reductase A) involved in chlorophyll biosynthesis, five transcripts (Cs6g15910, Phytoene synthase; Orange1.1t02108, phytoene synthase 2; Cs6g13340, Prolycopene isomerase 1; Cs4g14850/Orange1. 1t00772, Capsanthin/capsorubin synthase and Orange_new Gene_1755, Lycopene beta-cyclase) involved in carotenoid biosynthesis and four transcripts (Cs2g03270/Cs5g14370, 9-cis-epoxycarotenoid dioxygenase 2; Cs6g01180, Xanthoxin dehydrogenase; Cs7g14820, Carotenoid cleavage dioxygenase 4a and Cs6g19380, ABA 8&apos; -hydroxylase) involved in ABA metabolism were obtained (Fig. 2 and Table 4).

According to the result of expression and annotation analyses, thirteen transcripts i.e., BCHP, CRD1, CHLM, CHLH1, HEMFI/HEMF2, FC1, DVR, CAO, CLH, CCS, PSY, AB and NCED1 associated with chlorophyll metabolism, carotenoid biosynthesis and ABA metabolism were obtained with a Fold Change ≥2 and FDR <0.01 as screening standard (Table 5).

Expression analyses of candidate transcripts

Expression patterns of candidate transcripts associated with chlorophyll metabolism were analyzed between WZSTJ and HWWZSTJ at all fruit maturation stages (Fig. 3). Compared with WZSTJ, lower expression levels of ALAD1 and CLH were detected in HWWZSTJ at all fruit maturation stages. Expression of ALAD1 and CLH were increasing before fruit maturation and decreased thereafter in both WZSTJ and HWWZSTJ. The highest expression level of CLH was detected on the 295^th DAF in HWWZSTJ, which was 20 d later than WZSTJ. Expression levels of CAO1 and PAO in HWWZSTJ was higher than that in WZSTJ. FC1 showed a decrease trend during fruit maturation of WZSTJ and HWWZSTJ. As for GluRS, HEMF1, HEMG and CHLM, they showed irregular expression patterns in WZSTJ and HWWZSTJ (Fig. 3).

Six carotenoid biosynthesis-related transcripts showed a trend from rise to decline at all fruit maturation stages of WZSTJ and HWWZSTJ (Fig. 4). The highest expression level of CCS was detected on the 295^th DAF in HWWZSTJ, which was 20 d later than that of WZSTJ. Expression levels of PDS1, PSY3, PSY5, PSY6 and PSY7 in WZSTJ were higher than that of HWWZSTJ. PDS1 showed an increasing trend during fruit maturation of WZSTJ and HWWZSTJ and reached its maximum expression on the 295^th DAF. PSY5 showed the highest expression levels on the 275^th DAF compared to the highest expression levels of PSY3, PSY6 and PSY7 on the 265^th DAF in both WZSTJ and HWWZSTJ. Expression levels of PSY5 were increasing before the 275^th DAF and decreased thereafter. PSY3, PSY6 and PSY7 were up-regulation before the 265^th DAF and decreased gradually thereafter (Fig. 4).

Expression patterns of two candidate transcripts i.e., AB1 and NCED1 related to ABA metabolism were analyzed at all fruit maturation stages of WZSTJ and HWWZSTJ (Fig. 5). AB1 showed a trend from rise to decline during fruit maturation stages of WZSTJ and HWWZSTJ. The highest expression level of AB1 was obtained on the 295^th DAF in HWWZSTJ, which was 20 d later than WZSTJ. Similar expression patterns of NCED1 were observed before the 295^th DAF in HWWZSTJ and WZSTJ (Fig. 5). The expression level of NCED1 in HWWZSTJ was lower than that of WZSTJ during 275^th DAF to 305^th DAF. The highest expression level of NCED1 was detected on the 305^th DAF of WZSTJ and significantly deceased thereafter (Fig. 5). However, the highest expression of NCED1 was at 295 DAF in HWWZSTJ. Results from expression analyses of candidate genes suggested that NCED1 might play a leading role in late-maturing characteristics of HWWZSTJ.

Cloning and phylogenetic analyses of candidate genes

Full-length cDNA sequences of CLH and DVR from chlorophyll metabolism, PSY3, PSY5, PSY6, PSY7 and CCS from carotenoid biosynthesis, AB1 and NCED1 from ABA metabolism were cloned from HWWZSTJ and WZSTJ mandarins. There was one difference in base pair of CLH, PSY3 and PSY5 cDNA sequences between HWWZSTJ and WZSTJ (Figs. S5–S7). However, the amino acid sequences of CLH, PSY3 and PSY5 from HWWZSTJ was 100% identical to that from WZSTJ. There were 4, 6, 4, 3 and 17 bp difference between the sequences of DVR, CCS, PSY6, PSY7 and AB1 derived from HWWZSTJ and WZSTJ and this resulted in 2, 3, 3, 1 and 8 differences in the amino acids that would have been incorporated during translation of these transcripts (Figs. S8–S12). Compared with WZSTJ,there were 270 consecutive bases missing in cDNA sequence of the NCED1 from HWWZSTJ (Fig. 6). Phylogenetic analysis showed that CLH, DVR, PSY and NCED1 belonged to the same cluster, and their homology in comparison with similar sequences derived from other species is depicted in Figs. S13–S16. Results from sequence analyses suggested that deletion of 270 nucleotides in NCED1 maybe result in late-maturing characteristics of HWWZSTJ.