Identification of Candidate Anthocyanin-Related Genes by Transcriptomic Analysis of ‘Furongli’ Plum (Prunus salicina Lindl.) during Fruit Ripening Using RNA-Seq

Anthocyanins are important pigments and are responsible for red coloration in plums. However, little is known about the molecular mechanisms underlying anthocyanin accumulation in plum fruits. In this study, the RNA-seq technique was used to analyze the transcriptomic changes during fruit ripening in the red-fleshed plum (Prunus salicina Lindl.) cultivar ‘Furongli’. Over 161 million high-quality reads were assembled into 52,093 unigenes and 49.4% of these were annotated using public databases. Of these, 25,681 unigenes had significant hits to the sequences in the NCBI Nr database, 17,203 unigenes showed significant similarity to known proteins in the Swiss-Prot database and 5816 and 8585 unigenes had significant similarity to existing sequences in the Kyoto Encyclopedia of Genes and Genomes and the Cluster of Orthologous Groups databases, respectively. A total of 3548 unigenes were differentially expressed during fruit ripening and 119 of these were annotated as involved in “biosynthesis of other secondary metabolites.” Biological pathway analysis and gene ontology term enrichment analysis revealed that 13 differentially expressed genes are involved in anthocyanin biosynthesis. Furthermore, transcription factors such as MYB and bHLH, which may control anthocyanin biosynthesis, were identified through coexpression analysis of transcription factors, and structural genes. Real-time qPCR analysis of candidate genes showed good correlation with the transcriptome data. These results contribute to our understanding of the molecular mechanisms underlying anthocyanin biosynthesis in plum flesh. The transcriptomic data generated in this study provide a basis for further studies of fruit ripening in plum.


INTRODUCTION
The plum is one of the traditional fruit trees in China and is widely distributed in the world (Carrasco et al., 2012). Its fruit is highly appreciated by consumers. Plums are rich in bioactive substances such as vitamin C, carotenoids, polyphenols, and anthocyanins (Valero et al., 2013) and have great health benefits (Hooshmand and Arjmandi, 2009;Lee et al., 2009;Shukitt-Hale et al., 2009;Johnson et al., 2011). Color is one of the most important determinants of fruit quality. As a result of the influence of culture, Chinese people prefer red color. Anthocyanin accumulation is responsible for red coloration in plums (Cheng et al., 2015). Anthocyanins are widespread secondary metabolites that play an important role in the pigmentation of fruits. They are powerful antioxidants and naturally-occurring dietary anthocyanins are beneficial to human health (Nemie-Feyissa et al., 2015). Santhakumar et al. (2015) demonstrated that consumption of anthocyanin-rich plum juice attenuates thrombogenesis by reducing platelet activation/hypercoagulability and oxidative stress. It is therefore, desirable to increase the anthocyanin content in plums, especially in the flesh, through improved cultivation methods, postharvest handling, or breeding.
In addition, other regulatory factors also affect anthocyanin biosynthesis via interaction with MBW complexex or by modulating the transcription of structural genes directly.
MYBs that act as repressors of the anthocyanin pathway have been identified in several plants (Matsui et al., 2008;Salvatierra et al., 2013;Xu et al., 2013;Huang et al., 2014;Jun et al., 2015;Yoshida et al., 2015). Shin et al. (2007) found that PIF3 and HY5 regulated anthocyanin synthesis by binding directly to promoters of anthocyanin biosynthetic genes in Arabidopsis thaliana. The Arabidopsis COP1/SPA complex represses anthocyanin accumulation by degrading PAP1 and PAP2 proteins (Maier et al., 2013;Maier and Hoecker, 2014). Li et al. (2012b) also demonstrated that MdCOP1 negatively regulates accumulation of anthocyanin in apple peel by modulating the degradation of the MdMYB1 protein. ANAC078 has been shown to promote accumulation of anthocyanins by inducing the expression of flavonoid biosynthesis genes under high-light (Morishita et al., 2009). Recently, Zhou et al. (2015) found that NAC transcription factor BL, which controls the red-fleshed trait, interacts with PpNAC1 to form a heterodimer, and activate the transcription of PpMYB10.1 to induce anthocyanin accumulation. This process is repressed by a SQUAMOSA promoter-binding protein-like (SPL) transcription factor PpSPL1. SPLs have been shown to inhibit the expression of anthocyanin biosynthetic genes and negatively regulate anthocyanin accumulation through destabilization of MBW (Gou et al., 2011). MADS-box genes are reported to be involved in regulation of anthocyanin accumulation (Jaakola et al., 2010;Wu et al., 2013). Jasmonates induce anthocyanin biosynthesis through degradation of jasmonate-ZIM-domain (JAZ) proteins and the subsequent release of MBW (Qi et al., 2011). DELLA proteins promote anthocyanin biosynthesis by sequestering MYBL2 and JAZ proteins, which repress the activity of MBW, to release bHLH/MYB subunits in Arabidopsis (Xie et al., 2016). Furthermore, epigenetic mechanisms also play important roles in anthocyanin biosynthesis Zabala and Vodkin, 2014).
The 'Furongli' plum (Prunus salicina Lindl.), a red-skinned and red-fleshed cultivar, is native to Fujian, where it has been cultivated for more than 700 years. Fruit of 'Furongli' is popular for its attractive color, delicious taste, and healthpromoting nutrients (Fang et al., 2014). The fruit can be eaten fresh and is used for the production of candied fruits. In recent years, RNA-seq-based transcriptome analysis has been extensively used for identification of functional genes in fruit trees. Rodamilans et al. (2014) analyzed the transcriptome changes in response to Plum pox virus infection in P. cerasifera. Jo et al. (2015) reported the leaf transcriptome assembly of two different P. salicina cultivars. More recently, González et al. (2016a) reported the fruit skin transcriptomes of two Japanese plum cultivars with different skin color and developed candidate EST-SSR markers for marker-assisted selection of fruit skin color in the Japanese plum (González et al., 2016b). In the present study, we analyzed the transcriptomic changes during the ripening of 'Furongli' plum fruits to identify candidate genes involved in the biosynthesis of anthocyanins. Based on the RNA-seq datasets generated, we identify several potential structural genes, and transcription factor genes related to anthocyanin biosynthesis. The transcriptomic data generated in this study provide a basis for further studies of fruit ripening in plum, and identify candidate genes involved in sugar accumulation, organic acid degradation, fruit softening, and pigmentation.

Plant Materials
All samples were collected from 5-year-old field grown 'Furongli' plum (Prunus salicina Lindl.) trees in an orchard in Yongtai County, Fujian Province, China. Fruit samples were harvested at 105, 115, 125, and 135 days after flowering (DAF) from three trees in 2014 ( Figure 1A). Twenty representative fruits were sampled from each tree at each developmental stage and sliced. The sliced samples were combined and immediately frozen in liquid nitrogen and kept at −80 • C until use.

Determination of Anthocyanin Content
Anthocyanin content was quantified as described by Niu et al. (2010). Briefly, approximately 3 g of sample was ground to a fine powder in liquid nitrogen and extracted with 20 mL extraction solution (0.05% HCl in methanol) at 4 • C for 24 h. After centrifugation at 8000 × g for 20 min, the supernatant was transferred into a clean tube. One milliliter of supernatant and 4 mL of either bufferA (0.4MKCl, adjusted to pH 1.0 with HCl) or buffer B (1.2 N citric acid, adjusted to pH 4.5 with Na 2 HPO 4 ) were mixed and the absorbance at 510 and 700 nm (A 510 and A 700 ) measured for A and B buffers, respectively. The anthocyanin content was calculated according to Romero et al. (2008) using the following formula: TA = A × MW × 5 × 100 × V/e, where TA stands for total anthocyanin content (mg/100 g, as cyanidin-3-O-glucose equivalent), V for final volume (mL), and A = [A 510 (pH 1.0) − A 700 (pH 1.0)] − [A 510 (pH 4.5) − A 700 (pH 4.5)]. A molar absorptivity (e) of 26,900 and molecular weight (MW) of 449.2 were used according to Wrolstad et al. (1982). Three measurements were taken for each biological replicate.

RNA preparation, Library Construction, and RNA-Seq
Total RNA was extracted from each sample using a EZNA Plant RNA Kit (Omega Bio-tek). The concentration of RNA was quantified with a Qubit 2.0 Fluorometer (Invitrogen, Life Technologies, CA, USA), and RNA integrity was evaluated using an Agilent 2100 Bioanalyzer (Agilent Technologies, Santa Clara, CA, USA). Equal amounts of RNA from three samples at the same stage were mixed together. The RNA isolation, library construction, and RNA-seq were performed by staff at Beijing BioMarker Technologies (Beijing, China). cDNA libraries were constructed as described by Han et al. (2013). Poly-A mRNA was enriched using poly-T oligoattached magnetic beads, then broken into small pieces, and used as template for cDNA synthesis. First strand cDNA was synthesized using reverse transcriptase and random primers, followed by second strand cDNA synthesis using DNA Polymerase I and RNase H. The cDNA libraries were sequenced on an Illumina HiSeq TM 2500.

De novo Transcriptome Assembly
Before assembly, the raw reads in FASTQ format were filtered using in-house Perl scripts to discarding the reads containing the sequencing adapters and low-quality reads with ambiguous "N" bases or in which more than 10% of bases had a Q ≤ 20. The left files (read1 files) from all libraries/samples were pooled into one large left.fq file, and right files (read2 files) into one large right.fq file. De novo assembly of the clean reads was performed using Trinity (Grabherr et al., 2011) with min_kmer_cov set to 2 by default and all other parameters set to their default values.

Functional Annotation and Classification
To annotate unigene sequences of 'Furongli' plums, blastx search (E < 10 −5 ) was used to search against the NCBI nonredundant protein (nr), UniProt/Swiss-Prot, Gene Ontology (GO), Cluster of Orthologous Groups of proteins (COG), and Kyoto Encyclopedia of Genes and Genomes (KEGG) and Swiss-Prot databases and retrieve protein functional annotations based on sequence similarity. GO terms were assigned to the unigenes using Blast2GO (Conesa et al., 2005) with E ≤ 10 −5 . The distribution of the GO functional classifications of the unigenes was plotted using WEGO software (Ye et al., 2006).

Unigene Differential Expression Analysis
Reads from the four cDNA libraries were mapped to the assembled unigenes using Bowtie (Langmead et al., 2009). Unigene expression levels were quantified using fragments per kilobase of transcript per million mapped reads (FPKM ). FPKM values were calculated using RSEM (Li and Dewey, 2011). For each sequenced library, the read counts were adjusted using edgeR  with one normalization factor. Differential expression analysis of two samples was performed using the DEGseq R package . Unigenes differentially expressed between two samples were screened using false discovery rate (FDR) < 0.01 and absolute log 2 fold changes ≥1 as the threshold.

Prediction of Transcription Factors
To identify transcription factors expressed during ripening of 'Furongli' plum fruit, predicted peptide sequences of all unigenes were searched against transcription factors in PlantTFDB 3.0 using the Transcription Factor Prediction module (http://planttfdb.cbi.pku.edu.cn/prediction.php) with default parameters.

Correlation Analysis of Structural Genes and Transcription Factors
Correlation analysis of anthocyanin structural genes and transcription factors was performed as described by Ye et al. (2015). To exclude false positives, we only selected unigenes with a FPKM value ≥10 in at least one of the four stages during fruit ripening. Transcription factors with correlation coefficient values of ≥0.95 by t-test were considered to have an expression that was significantly correlated with the expression of genes in anthocyanin biosynthetic pathways. The formula used was t = (r √ (n-2))/ √ (1-r 2 ), at P < 0.05 and n = 4. A value of |t|>t 0.05, 2 = 4.303, implying r > 0.95, means significant correlation. Coexpression analysis was carried out using the "CORREL" function in Microsoft Excel 2003 and confirmed using in-house Perl scripts and IBM SPSS statistics software.

Real-Time Quantitative RT-PCR Analysis
Nineteen candidate differentially expressed genes involved in anthocyanin biosynthesis were selected for validation by realtime quantitative RT-PCR (qRT-PCR). Total RNA from fruit samples was extracted using a modified CTAB method as described by Xuan et al. (2015). The primer sequences used for qRT-PCR are listed in Table S1. The cDNA was transcribed from 500 ng of total RNA using the PrimeScript RT reagent kit with gDNA Eraser (Takara, Dalian, China). Quantitative RT-PCR was performed using the Eppendorf RealPlex 4 real-time PCR system (Hamburg, Germany) in a total volume of 20 µL in each well containing 10 µL of 2 × SYBR Premix Ex Taq TM II (Tli RNaseH Plus, TaKaRa), 1 µL of cDNA (in 1:10 dilution), and 0.4 µL 10 µM primers. Quantitative PCR conditions were 5 min at 95 • C, followed by 40 cycles of 5 s at 95 • C, 15 s at annealing temperatures listed in Table S1, and 30 s at 72 • C, followed by 60-95 • C melting curve detection. The actin gene was used as the reference. The expression levels were calculated as described by Fang et al. (2013). Three biological and three technical replicates were performed.

Anthocyanin Accumulation during Ripening of 'Furongli' Plum Fruits
As indicated in Figure 1A, the color of 'Furongli' plums changed from green to red during ripening and the flesh became pigmented before the skin. Anthocyanin accumulation is responsible for the red color of plums and the major anthocyanins in 'Furongli' plums are cyanidin 3-rutinoside and cyanidin 3-glucoside (Usenik et al., 2009;Zhang, 2009). The anthocyanin content of 'Furongli' plums increased from 0.13 mg/100 g FW to 7.0 mg/100 g FW as ripening proceeded ( Figure 1B).

RNA-Seq and De novo Transcriptome Assembly
Four cDNA libraries were constructed from the total RNA of 'Furongli' plums at 105, 115, 125, and 135 DAF. These libraries were subjected to RNA-seq using an Illumina HiSeq2500, generating 41,555,484, 43,021,234, 40,991,934, and 37,023,409 raw 100-bp paired-end raw reads, respectively ( Table 1). All of the raw reads are available in the NCBI SRA database (accession number SRP076001). After removing low-quality reads and trimming adapter sequences, 41,067,827, 42,872,666, 40,886,898 and 36,927,707 clean reads were obtained for the 105, 115, 125, and 135 DAF libraries, respectively.
Using Trinity, the clean reads from the four libraries were assembled into 100,711 transcripts with an average length of 1320 bp, and 52,093 unigenes with an average length of 872 bp were obtained ( Table 1). Unigenes shorter than 500 bp accounted for 55.43% of the total unigenes and 26.77% of the unigenes (13,944) were longer than 1 kb (Figure 2). These results suggest that the quality of the unigene data was high enough for the following analyses.

Functional Annotation and Classification
To annotate the transcriptome of 'Furongli' plums, 52,093 unigenes were searched against five public databases (nr, UniProt/Swiss-Prot, GO, COG, and KEGG) with a cutoff E-value of 10 −5 . The functional annotation results are listed in Table 2.
Only 49.4% of the unigenes (25,730) were identified. The remaining unigenes (50.6%) could not be annotated with known genes ( Table 2), most likely because oflimitations in the genomic information and the presence of short sequences. Only 947 (3.6%) of the unannotated unigenes were longer than 1000 bp, while 47.1% were shorter than 300 bp. The species distribution of the best match results in nr is indicated in Figure 3. The 'Furongli' plum unigenes showed the closest matches with Prunus persica (77.80%), followed by Fragaria vesca (5.32), Vitis vinifera (1.51), Zea mays (1.06), Theobroma cacao (0.71), Populus  The 'Furongli' plum unigenes were searched against the GO database to classify standardized gene functions. At least one GO term was assigned to 18,623 of the unigenes ( Table 2). Unigenes were assigned to three main GO categories, including biological process category, cellular component category, and molecular function category, and 58 subcategories shown in Figure 4. The terms "cell, " "cell part, " and "organelle" were dominant in the cellular component category, the term "binding" and "catalytic activity" was dominant in the molecular function category, and the terms "metabolic process" and "cellular process" were dominant in the biological process category.
Only 5816 (11.2%) unigenes had significant matches in the KEGG database and these were assigned to 115 KEGG pathways (Table S2). Among them, 2009 unigenes were assigned to metabolic pathways. As demonstrated in Figure 6A, 552 unigenes were assigned to carbohydrate metabolism, followed by Energy metabolism (548 unigenes), amino acid metabolism (426 unigenes), and lipid metabolism (195 unigenes). We concentrated on the "biosynthesis of other secondary metabolites" category in relation to fruit pigmentation. In this category, 119 unigenes were classified into seven subcategories ( Figure 6B). Among these, the cluster for "phenylpropanoid biosynthesis" was the most represented, followed by "flavonoid biosynthesis." "Flavone and flavonol biosynthesis" and "caffeine metabolism" appeared to be the smallest groups.

Changes in Gene Expression during Fruit Ripening
To study unigene expression during different developmental stages, the reads from the four libraries were mapped to the assembled transcriptome. As indicated in Figure 7, most of the unigenes (32,926, 63.21%) were expressed in all four stages of development. A total of 3891 unigenes (7.47%) were only detected at a single stage of development. Of these stage-specific
Transcription factors play important roles in the regulation of anthocyanin biosynthesis. In total, 791 unigenes (Table  S5) were predicted to encode transcription factors from 55 different families (Table S6) and 147 of them were differentially expressed (Table S7). To identify transcription factors that were coexpressed with the candidate enzymatic genes involved in anthocyanin biosynthesis, a transcription abundance correlation analysis was carried out between the differentially expressed transcription factors and structural genes from the anthocyanin biosynthetic pathway. This identified 37 transcription factors whose expression levels were highly correlated with those of the candidate structural genes ( Table 4; Table S8). Of these, 22 showed a significant correlation with five or more structural genes from the anthocyanin biosynthetic pathway. The identified transcription factors included homologs of Arabidopsis transcription factors that are implicated in regulating anthocyanin biosynthesis, such as MYB, bHLH, and NAC ( Table 3).
A total of 37 MYBs were differentially expressed during ripening of 'Furongli' plums and three of them (c39005.graph_c0, c29499.graph_c0 and c32850.graph_c0) were associated with the anthocyanin biosynthetic pathway ( Table 4; Table S8). The unigene c39005.graph_c0 (homologous to AtMYB113) was upregulated, while c29499.graph_c0 (homologous to AtMYB73) and c32850.graph_c0 (homologous to AtMYB102) were downregulated. The expression level of the homolog of AtMYBD (c28480.graph_c0) also increased. However, it was correlated with none of the structural genes. Seven of the differentially expressed transcription factors annotated as bHLH showed   Table 4; Table S8). Only two of the bHLH genes (c7988.graph_c0 and c18575.graph_c1) showed a positive correlation with the expression of structural genes, while most of them were negatively correlated with that of structural genes involved in the anthocyanin biosynthetic pathway. In addition, a plum bHLH (c36695.graph_c0, log2 fold change <1.0), which is the best BLAST match to Arabidopsis AtbHLH42, was significantly correlated with seven anthocyanin biosynthetic structural genes (Table S8). WD40 encoding unigenes (c28377.graph_c0 and c10590.graph_c1), which are the homolog of AtTTG1, did not show differential expression during the ripening process. Apart from the MBW components, other differentially expressed transcription factors, such as NAC, were also found to be potentially related to the anthocyanin pathway. Four plum NAC genes (c19087.graph_c0, c19209.graph_c0, c27539.graph_c0, and c37766.graph_c1) were significantly correlated with structural genes ( Table 4; Table S8). The unigenes c27539.graph_c0 (a  homolog of AtNAC100) and c19209.graph_c0 were upregulated during plum fruit ripening. Peach homolog of AtNAC100 have been reported to be involved in the regulation of anthocyanin accumulation in fruit flesh. GO annotation results indicated that c19209.graph_c0 is involved in "biological process: positive regulation of flavonoid biosynthetic process" (GO:0009963).
We further analyzed the expression profiles of 19 candidate unigenes (13 structural genes and six transcription factors) involved in anthocyanin biosynthesis using qRT-PCR. The results indicated that there is a good correlation between RNA-seq data and qPCR data for most of the genes (Figure 9).

DISCUSSION
In this study, we constructed a transcriptome of 'Furongli' plums during fruit maturation. In total, 52,093 unigenes were assembled with a mean length of 872 bp, which is comparable to 944 bp for sweet cherry (P. avium L.) (Wei et al., 2015). A large quantity of genomic data is available for many rosaceous plants, but only 49.4% of the plum unigenes were annotated to public databases (nr, Swiss-Prot, GO, COG, and KEGG). This means that more than half of the unigenes have no significant known homologs. The low rate of annotated unigenes could be a result of limitations in the genomic information available for P. salicina Lindl. as is the case in other non-model plant species (Yates et al., 2014;Wei et al., 2015;Wu et al., 2015). As expected, 77.80% of the unigenes annotated using nr show significant similarity to P. persica transcripts. The unannotated unigenes could be plum specific genes with novel functions. The transcriptome of 'Furongli' plums will serve as an important datasets for studying plum ripening processes such as sugar accumulation, organic acid degradation, fruit softening, and pigmentation.
Anthocyanin-rich plums are of great interest for their implications in human health (Santhakumar et al., 2015). The main objective of this study was to identify genes involved in anthocyanin biosynthesis in plums. RNA-seq-based comparative transcriptome analysis has been shown to be an efficient strategy for the investigation of genes involved in anthocyanin biosynthesis in several plants, such as kiwifruit (Li et al., 2015a), sweet cherry (Wei et al., 2015), zoysiagrass (Ahn et al., 2015), anthurium (Li et al., 2015b), and potato . In the later ripening stages, 'Furongli' plums accumulate anthocyanins rapidly (Figure 1). The expression of anthocyanin biosynthetic genes has been shown to be correlated with fruit anthocyanin content in Rosaceae such as apple (Feng et al., 2013(Feng et al., , 2014Vimolmangkang et al., 2014), pear (Li et al., 2012a;Yang et al., 2015), sweet cherry (Wei et al., 2015), strawberry (Xu et al., 2014), and plum (Cheng et al., 2015). In the present study, changes in gene expression between different stages of ripening were analyzed to identify differentially expressed genes implicated in anthocyanin biosynthesis, including PAL, C4H, 4CL, CHS, CHI, F3H, F3 ′ H, DFR, ANS/LDOX, UFGT, and GST, were significantly upregulated in the late stages of fruit maturation (Table 3; Figure 9).
Anthocyanin biosynthesis is regulated by several well-studied transcription factors such as MYB, bHLH, and WD40 (Gonzalez et al., 2008). MYB transcription factors have been reported to play a pivotal role in anthocyanin biosynthesis regulation in several fruit trees (Chagné et al., 2013;Ravaglia et al., 2013;Umemura et al., 2013;Lai et al., 2014;Shen et al., 2014;Tuan et al., 2015;Zhai et al., 2015;Jin et al., 2016). Lin-  demonstrated that R2R3 MYBs are highly conserved in rosaceous plants and MYBs from European plum and cherry plum are able to induce the anthocyanin accumulation in tobacco. Gu et al. (2015) proposed that constitutive activation of PcMYB10.6 is responsible for red pigmentation in purpleleaf plum. Cheng et al. (2015) indicated that PsMYB10 was involved in ethylene-regulated anthocyanin biosynthesis in plums. These studies suggested that MYBs play a role in anthocyanin biosynthesis in the plum. In this study, a plum MYB (c39005.graph_c0) was positively correlated with structural genes. Its Arabidopsis homolog (AtMYB113) is an activator of the anthocyanin pathway. Recently, AtMYBD was shown to enhance anthocyanin biosynthesis by repressing the negative regulator MYBL2 (Nguyen et al., 2015). The expression of an AtMYBD homolog (c28480.graph_c0) was upregulated in late ripening stages, but it does not show significant correlation with any structural genes. In addition, the anthocyanin biosynthetic pathway is also negatively controlled by MYB repressors in many plants, including Arabidopsis (Dubos et al., 2008;Matsui et al., 2008;Zhu et al., 2009), poplar (Yoshida et al., 2015), Medicago truncatula (Jun et al., 2015), grapevine (Cavallini et al., 2015;Pérez-Díaz et al., 2015), strawberry (Salvatierra et al., 2013), and apple (Lin- . We found that two MYBs (c29499.graph_c0 and c32850.graph_c0) were negatively correlated with anthocyanin biosynthetic genes. Another important component of the MBW complex, bHLH proteins have been shown to be essential for or to enhance MYB-induced anthocyanin accumulation in transient expression FIGURE 9 | Expression analysis of 19 differentially expressed genes related to anthocyanin biosynthesis in 'Furongli' plums during fruit ripening. Actin was used as the internal control. The error bars represent the standard error of three biological replicates. The numbers above the graphics correspond to values obtained with the Pearson correlation. Pearson correlation between the RNA-seq data and qRT-PCR data was calculated using the log2 value of FPKM and the relative expression level.
assays (Liu et al., 2013b;Rahim et al., 2014;Feng et al., 2015;Starkevič et al., 2015;Wei et al., 2015;Lai et al., 2016). Our results indicated that c36695.graph_c0, which is highly homologous to AtTT8 (AT4G09820), accumulates to higher levels in late stages of ripening and shows significant correlation with anthocyanin biosynthetic genes. Conversely, the expression of c33382.graph_c0 was repressed as ripening proceeded. The unigene c33382.graph_c0 is a homolog of AtbHLH14 (AT4G00870), which belongs to bHLH subgroup IIId. Song et al. (2013) demonstrated that Arabidopsis lines overexpressing bHLH17 showed jasmonate-induced anthocyanin accumulation. NAC proteins have also been reported to be involved in anthocyanin synthesis in Arabidopsis (Morishita et al., 2009). Recently, a peach homolog of AtNAC100 was shown to be responsible for regulation of anthocyanin accumulation in flesh . Our results indicated that a plum homolog of AtNAC100 (c27539.graph_c0) is upregulated and positively correlated with anthocyanin biosynthetic genes. It should be noted that coexpression analysis usually requires a large sample size and the small sample size in our study will reduce the reliability of our results. However, relatively small numbers of samples have recently been used to analyze the correlation of structural genes and regulators involved in biological processes, such as flavonoid biosynthesis (Zhai et al., 2015) and fruit ripening (Wu et al., 2016). The exact roles of these candidate transcription factor should be investigated in further studies.
In the current study, we used RNA-seq to analyze changes in the transcriptome during ripening of 'Furongli' plums. We generated 52,093 unigenes and over 50% of them were not annotated to public databases. Unigenes differentially expressed during fruit ripening were identified. Candidate genes encoding anthocyanin biosynthetic enzymes and transcription factors involved in anthocyanin biosynthesis were identified using functional annotation and coexpression analysis of differentially expressed genes. The expression patterns of some candidate genes encode anthocyanin biosynthetic enzymes and transcription factors were further validated by qRT-PCR. This provides an important datasets for studying fruit ripening processes, especially anthocyanin biosynthesis, in plums. Further studies are needed to determine whether the identified candidate genes are related to anthocyanin biosynthesis in plum.

AUTHOR CONTRIBUTIONS
This study was conceived by ZF and XY. The plant material preparation were carried out by ZF and SP. ZF, XY, CJ, and DZ analyzed the RNA-seq data. ZF, CJ, DZ, and SP performed the laboratory experiments and analyses. ZF drafted the manuscript. XY revised the manuscript. All authors read and approved the final manuscript.