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

Characteristics of late maturity in fruit are good agronomic traits for extending the harvest period and marketing time. However, underlying molecular basis of the late-maturing mechanism in fruit is largely unknown. In this study, RNA sequencing (RNA-Seq) technology was used to identify differentially expressed genes (DEGs) related to late-maturing characteristics from a late-maturing mutant ‘Huawan Wuzishatangju’ (HWWZSTJ) (Citrus reticulata Blanco) and its original line ‘Wuzishatangju’ (WZSTJ). A total of approximately 17.0 Gb and 84.2 M paried-end reads were obtained. DEGs were significantly enriched in the pathway of photosynthesis, phenylpropanoid biosynthesis, carotenoid biosynthesis, chlorophyll and abscisic acid (ABA) metabolism. Thirteen candidate transcripts related to chlorophyll metabolism, carotenoid biosynthesis and ABA metabolism were analyzed using real-time quantitative PCR (qPCR) at all fruit maturing stages of HWWZSTJ and WZSTJ. Chlorophyllase (CLH) and divinyl reductase (DVR) from chlorophyll metabolism, phytoene synthase (PSY) and capsanthin/capsorubin synthase (CCS) from carotenoid biosynthesis, and abscisic acid 8′-hydroxylase (AB1) and 9-cis-epoxycarotenoid dioxygenase (NCED1) from ABA metabolism were cloned and analyzed. The expression pattern of NCED1 indicated its role in the late-maturing characteristics of HWWZSTJ. There were 270 consecutive bases missing in HWWZSTJ in comparison with full-length sequences of NCED1 cDNA from WZSTJ. Those results suggested that NCED1 might play an important role in the late maturity of HWWZSTJ. This study provides new information on complex process that results in the late maturity of Citrus fruit at the transcriptional level.


INTRODUCTION
Fruit maturity date is an important economic trait and selection of varieties with different harvest time would be advantageous to extend their storage period and market share.
Citrus, one of the most important fruit crops, is a large-scale commercial production in the tropical and subtropical regions of the world. The total harvested area of citrus exceeds 8.8 million ha, with an annual yield of over 130 million tons in 2015 (Food and Agricultural Organization of the United Nations, 2014). Currently, harvest time for most citrus is mainly from November to December resulting in huge market pressure. Therefore, breeding of early-and late-maturing citrus varieties is essential to extend marketing season, meet the needs of consumers and ensure an optimal adaptation to climatic and geographic conditions.
Plant hormones play important roles in the regulation of fruit development and ripening (Kumar, Khurana & Sharma, 2014). Ethylene is known to be the major hormonal regulator in climacteric fruit ripening. In addition to ethylene, abscisic acid (ABA), auxin, gibberellin (GA) and brassinosteroid are involved in regulating fruit ripening. ABA plays an important role as an inducer along with ethylene signaling for the onset of fruit degreening and carotenoid biosynthesis during development and ripening process in climacteric and non-climacteric fruits (Leng et al., 2009;Sun et al., 2010;Jia et al., 2011;Romero, Lafuente & Rodrigo, 2012;Soto et al., 2013;Wang et al., 2016). ABA treatment can rapidly induce flavonol and anthocyanin accumulation in berry skins of the Cabernet Sauvignon grape suggesting that ABA could stimulate berry ripening and ripening-related gene expression (Koyama, Sadamatsu & Goto-Yamamoto, 2010). ABA also participates in the regulation of fruit development and ripening of tomato (Zhang, Yuan & Leng, 2009;Sun et al., 2011), cucumber (Wang et al., 2013), strawberry (Jia et al., 2011), bilberry (Karppinen et al., 2013, citrus (Zhang et al., 2014) and grape (Nicolas et al., 2014). Recent studies showed that ABA is a positive regulator of ripening and exogenous ABA application could effectively regulate citrus fruit maturation (Wang et al., 2016). Those results suggest that ABA metabolism plays a crucial role in the regulation of fruit development and ripening. In addition, fruit deterioration and post-harvest processes might influence fruit quality and ripening process. However, there are few reports involved in those processes. α-mannosidase (α-Man) and β-D-N-acetylhexosaminidase (β-Hex) are the two ripening-specific N-glycan processing enzymes that have proved that their transcripts increased with non-climacteric fruit ripening and softening (Ghosh et al., 2011). Genetic results have proved that 9-cis-epoxycarotenoid dioxygenase (NCED) is the key enzyme in ABA metabolism in plants (Liotenberg, North & Marion-Poll, 1999;Luchi et al., 2001). NCED1 could initiate ABA biosynthesis at the beginning of fruit ripening in both peach and grape fruits (Zhang et al., 2009). Silence of FaNCED1 (encoding a key ABA synthesis enzyme) in strawberry fruit could cause the ABA levels to decrease significantly and uncolored fruits and this phenomenon could be rescued by application of exogenous ABA (Jia et al., 2011). Suppression of the expression of SLNCED1 could result in the delay of fruit softening and maturation in tomato (Sun et al., 2012). Overexpression of ABA-response element binding factors (SlAREB1) in tomato could regulate organic acid and sugar contents during tomato fruit development. Higher levels of organic acid, sugar contents and related-gene expression were detected in SlAREB1-overexpressing lines in fruit pericarp of mature tomato (Bastías et al., 2011). However, there is little information available about the role of NCED1 genes in citrus fruit maturation (Zhang et al., 2014).
Bud mutant selection is the most common method for creating novel cultivars in Citrus. The 'Huawan Wuzishatangju' (HWWZSTJ) mandarin is an excellent cultivar derived from a bud sport of a seedless cultivar 'Wuzishatangju' (WZSTJ). Fruits of HWWZSTJ are mature in late January to early February of the following year, which is approximately 30 d later than WZSTJ (Qin et al., 2013;Qin et al., 2015). Therefore, the late-maturing mutant and its original cultivar are excellent materials to identify and describe the molecular mechanism involved in citrus fruit maturation. In this study, the highly efficient RNA-Seq technology was used to identify differentially expressed genes (DEGs) between the late-maturing mutant HWWZSTJ and its original line WZSTJ mandarins. DEGs involved in carotenoid biosynthesis, chlorophyll degradation and ABA metabolism were characterized. The present work could help to reveal the molecular mechanism of late-maturing characteristics of citrus fruit at the transcriptional level.

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-yearold 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 TM 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 2 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).

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

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).

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

ABA metabolism Cs3g23530
Abscisic acid 8 -hydroxylase (AB)  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.

DISCUSSION
Chlorophyll degradation, carotenoid biosynthesis and ABA metabolism play important roles in regulating citrus fruit maturation through a series of related genes or special signal network (Zhang et al., 2014). In this study, RNA-Seq technology was used to screen DEGs between a late-maturing mandarin mutant HWWZSTJ and its wild type WZSTJ during fruit maturation. DEGs between any two of the three libraries were significantly enriched in biological processes such as photosynthesis, phenylpropanoid biosynthesis, carotenoid biosynthesis, chlorophyll metabolism, ABA metabolism, starch and sucrose metabolism ( Table 3). Thirteen maturing-related transcripts involved in carotenoid biosynthesis, chlorophyll degradation and ABA metabolism were selected for further analysis. CLH is the key enzyme catalyzing the first step in the chlorophyll degradation. It can catalyze the hydrolysis of ester bond to yield chlorophyllide and phytol in the chlorophyll breakdown pathway (Jacob-Wilk et al., 1999;Tsuchiya & Takamiya, 1999). Jacob-Wilk et al. (1999) isolated a CLH encoding an active chlorophyllase enzyme and verified the role of CLH in chlorophyll dephytylation by in vitro recombinant enzyme assays. Expression level of CLH in Valencia orange peel was low and constitutive and did not significantly increase during fruit development and ripening (Jacob-Wilk et al., 1999). In the present study, a CLH was obtained from the transcriptome dataset. No difference was detected in the amino acid sequences of CLH between HWWZSTJ and WZSTJ. Expression levels of CLH were increasing prior to citrus fruit maturing, decreasing 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 that of WZSTJ (Fig. 3). Similar results were also observed in peels of the late-maturing mutant from Fengjie72-1 navel orange (Liu et al., 2006) and Tardivo clementine mandarin (Distefano et al., 2009). Those results suggested that CLH may balance between chlorophyll synthesis and its breakdown (Jacob-Wilk et al., 1999).
Citrus is a complex source of carotenoids with the largest number of carotenoids (Kato et al., 2004). Carotenoid contents and compositions are main factors that affect peel color of most citrus fruits (Tadeo et al., 2008). PSY is a regulatory enzyme in carotenoid biosynthesis (Welsch et al., 2000). PSY is present at low expression level in unripe (green) melon fruit, reaches its highest levels when the fruit turns from green to orange and persists at lower levels during later ripening stages (Karvouni et al., 1995). Liu et al. (2006) studied the mechanism underlying the difference between Fengwan (a late-maturing mutant) navel orange and its original cultivar (Fengjie72-1). The highest expression levels of some carotenoid biosynthetic enzymes in the peels of the late-maturing mutant occurred 30 d later than that of the original cultivar (Liu et al., 2006). In this work, PSY showed a trend from rise to decline at all fruit maturation stages of the late-maturing mutant HWWZSTJ and its original line WZSTJ. The expression levels of PSY3, PSY5, PSY6 and PSY7 in HWWZSTJ were lower than that in WZSTJ. These results demonstrated that the mutation in HWWZSTJ influenced the transcriptional activation of PSY.
ABA can be considered as a ripening regulator during fruit maturation and ripening. NCED, a key enzyme involved in ABA biosynthesis, plays an important role in fruit ripening of avocado (Persea americana) (Chernys & Zeevaart, 2007), orange (Citrus sinensis) (Rodrigo, Alquezar & Zacarías, 2006), tomato (Solanum lycopersicum) (Nitsch et al., 2009;Zhang, Yuan & Leng, 2009), grape (Vitis vinifera) and peach (Prunus persica) (Zhang et al., 2009). The NCED1 were expressed only at the onset stage of ripening in peach and grape, when ABA content became high (Zhang et al., 2009). Zhang et al. (2014 studied the mechanism of a spontaneous late-maturing mutant of 'Jincheng' sweet orange and its wild type through the comparative analysis. The highest expression of CsNCED1 was at 215 DAA in WT. In our study, expression levels of NCED1 increased prior to fruit maturing and decreased significantly thereafter in both HWWZSTJ and WZSTJ. The highest expression level of NCED1 was detected on the 305 th DAF of WT (WZSTJ). Our results were consistent with previous findings that NCED1 plays the most important role in the ABA biosynthesis pathway during the fruit maturing process (Zhang et al., 2014). Deletion of nucleotides could cause a shift of the reading frame and truncated protein, which can result in natural mutants. Compared with the cDNA sequence of NCED1 from WZSTJ, there were 270 consecutive bases missing in HWWZSTJ (Fig. 6). Those results suggested that NCED1 might play an important role in late-maturing of HWWZSTJ. A high-efficient regeneration system for WZSTJ has been established (Wang et al., 2015) and further study on the role of NCED1 in citrus is being carried out through genetic engineering.

CONCLUSION
RNA-Seq technology was used to identify pigment-related genes from a late-maturing mandarin mutant HWWZSTJ and its original cultivar WZSTJ. Thirteen candidate transcripts related to chlorophyll metabolism, carotenoid biosynthesis and ABA metabolism were obtained. NCED1, a gene involved in ABA metabolism, is probably involved in the formation of late maturity of HWWZSTJ based on sequence and expression analyses. The present study opens up a new perspective to study the formation of late maturity in citrus fruit.