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.

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₂HPO₄) were mixed and the absorbance at 510 and 700 nm (A₅₁₀ and A₇₀₀) 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₅₁₀ (pH 1.0) − A₇₀₀ (pH 1.0)] − [A₅₁₀ (pH 4.5) − A₇₀₀ (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.

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™ 2500.

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.

To annotate unigene sequences of ‘Furongli’ plums, blastx search (E < 10⁻⁵) 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⁻⁵. The distribution of the GO functional classifications of the unigenes was plotted using WEGO software (Ye et al., 2006).

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 (Robinson et al., 2010) with one normalization factor. Differential expression analysis of two samples was performed using the DEGseq R package (Wang et al., 2010). Unigenes differentially expressed between two samples were screened using false discovery rate (FDR) < 0.01 and absolute log₂ fold changes ≥1 as the threshold.

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

Nineteen candidate differentially expressed genes involved in anthocyanin biosynthesis were selected for validation by real-time 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⁴ 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™ 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.

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

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.

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⁻⁵. 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 trichocarpa (0.56), Hordeum vulgare (0.48), Ricinus communis (0.46), Medicago truncatula (0.38), and Glycine max (0.37%).

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.

The completeness of the plum transcriptome library and the validity of the annotations were further evaluated by COG annotation. Out of the 25,730 annotated unigenes, 8585 (33.37%) were clustered into 24 COG categories (Figure 5). The cluster for “general functional prediction only” (2260, 26.32%) represented the largest group, followed by “translation, ribosomal structure and biogenesis” (1095, 12.75%), “replication, recombination, and repair” (1074, 12.51%), “transcription” (1038, 12.09%), “Signal transduction mechanisms” (887, 10.33%), “posttranslational modification, protein turnover, chaperones” (840, 9.78%), “carbohydrate transport and metabolism” (652, 7.59%), “amino acid transport and metabolism” (553, 6.44%), and “energy production and metabolism” (516, 6.01%). “cell motility” (9, 0.1%) and “nuclear structure” (1, 0.01%) represented the smallest groups. No unigene was assigned to “extracellular structures.”

To better understand biological pathways involved in ripening of ‘Furongli’ plum fruits, all unigenes were searched against the reference canonical pathways in the KEGG database. 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.

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 expressed unigenes, 682 (1.31), 862 (1.65), 1280 (2.46), and 1.067 (2.05%) unigenes were expressed exclusively at 105, 115, 125, and 135 DAF, respectively.

By comparing the libraries, 507 differentially expressed unigenes (304 upregulated and 203 downregulated) were identified between 105 and 115 DAF, with 1674 differentially expressed unigenes (778 upregulated and 886 downregulated) between 115 and 125 DAF, as well as 1004 differentially expressed unigenes (523 up-regulated and 481 downregulated) between 125 and 135 DAF (Figure 8). Between 105 and 125 DAF, there are 2609 differentially expressed unigenes (1289 up-regulated and 1320 downregulated) showed differential expression. Between 105 and 135 DAF, there were 1743 differentially expressed unigenes (876 upregulated and 867 downregulated) show different expression. Between 115 and 135 DAF, there are 940 differentially expressed unigenes (467 upregulated and 473 downregulated). When considered together, 3548 unigenes were differentially expressed during fruit ripening. This result suggests that the developmental period with the most dynamic changes in the transcriptome was between 115 and 125 DAF, with almost half of the differentially expressed genes (47.18%) showing significant changes during this period.

To identify unigenes involved in anthocyanin biosynthesis, KEGG functional enrichment was analyzed to characterize the functions of differentially expressed unigenes (Table S3). A total of 11 unigenes encoding enzymes, including PAL (c38398.graph_c0), C4H (c23939.graph_c0), 4CL (c30378.graph_c0), CHS (c37054.graph_c0), CHI (c28749.graph_c0), F3H (c23888.graph_c0), F3′H (c38186.graph_c0), DFR (c24831.graph_c0), ANS/LDOX (c29583.graph_c0), and UFGT (c19095.graph_c0), were assigned to the anthocyanin biosynthetic pathway based on KEGG (Table 3; Table S3). GO annotation was carried out to further identify anthocyanin-associated unigenes without annotation information in KEGG database (Table S4). The unigene c8193.graph_c0 (CHI) was assigned to “anthocyanin-containing compound biosynthetic process.” The unigene c25083.graph_c0 (flavonoid 3′-monooxygenase) was predicted to have flavonoid 3′, 5′-hydroxylase activity (GO:0033772). Furthermore, a GST (c29416.graph_c0), with best match to AtGSTF12 (AT5G17220), also showed differential expression. As shown in Table 3, the expression of these 13 unigenes was significantly upregulated during fruit maturation.

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 a significant correlation with anthocyanin biosynthetic genes (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).

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.