Transcriptome profiling reveals differential gene expression in proanthocyanidin biosynthesis associated with red/green skin color mutant of pear (Pyrus communis L.)

Anthocyanin concentration is the key determinant for red skin color in pear fruit. However, the molecular basis for development of red skin is complicated and has not been well-understood thus far. “Starkrimson” (Pyrus communis L.), an introduced red pear cultivated in the north of China and its green mutant provides a desirable red/green pair for identification of candidate genes involved in color variation. Here, we sequenced and annotated the transcriptome for the red/green color mutant at three stages of development using Illumina RNA-seq technology. The total number of mapped reads ranged from 26 to 46 million in six libraries. About 70.11–71.95% of clean reads could be mapped to the reference genome. Compared with green colored fruit, a total of 2230 differentially expressed genes (DEGs) were identified in red fruit. Gene Ontology (GO) terms were defined for 4886 differential transcripts involved in 15 Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways. Three DEGs were identified as candidate genes in the flavonoid pathway, LAR, ANR, and C3H. Tellingly, higher expression was found for genes encoding ANR and LAR in the green color mutant, promoting the proanthocyanidin (PA) pathway and leading to lower anthocyanin. MYB-binding cis-motifs were identified in the promoter region of LAR and ANR. Based on these findings, we speculate that the regulation of PA biosynthesis might be a key factor for this red/green color mutant. Besides the known MYB and MADS transcription families, two new families, AP2 and WRKY, were identified as having high correlation with anthocyanin biosynthesis in red skinned pear. In addition, qRT-PCR was used to confirm the transcriptome results for 17 DEGs, high correlation of gene expression, further proved that AP2 and WARK regulated the anthocyanin biosynthesis in red skinned “Starkrimson,” and ANR and LAR promote PA biosynthesis and contribute to the green skinned variant. This study can serve as a valuable new resource laying a solid foundation for functional gene identification in the anthocyanin pathway of red-skinned pear and provide a good reference for relevant research on molecular mechanisms of color variation in other pear species.


Introduction
Red coloration is an appealing feature in many flowers, fruits, and other plant tissues and is associated with anthocyanin accumulation. Anthocyanins also play an important role in plant disease resistance and protection against ultraviolet radiation (Bieza and Lois, 2001), and have been generally considered to have antioxidant capability (Veeriah et al., 2006), leading to the prevention of neuronal and cardiovascular illnesses, cancer, and diabetes in humans (Konczak and Zhang, 2004).
Pear fruit (Pyrus) is cultivated world-wide, with China ranking as the top producer of Asian pear. However, few red-skinned pears are cultivated and sold in China. Comparatively more redskinned cultivars are found in European pear (Pyrus communis), which provides a good resource for red skin color and the possibility for color improvement of Asian pear. It is wellknown that flavonoid biosynthesis, such as the anthocyanins in pear skins, is genetically catalyzed by structural genes encoding key enzymes (Holton and Cornish, 1995;Davies and Schwinn, 2003). In addition, these structural genes involved in anthocyanin biosynthesis have been reported to be controlled by a specific transcriptional complex, first reported in the model plant Arabidopsis thaliana (Zhang et al., 2003;Rowan et al., 2009). Subsequently, grapevine (Hichri et al., 2010), apple (Espley et al., 2007;Lin-Wang et al., 2011;Xie et al., 2012) Chinese bayberry (Liu et al., 2013c) have all been reported to have the transcription factors MYB, bHLH, and WD40 coordinated with each other to promote anthocyanin accumulation. Recently, it has also been reported that NAC transcription factors activate anthocyanin biosynthesis in blood-fleshed peach .
As for anthocyanin biosynthesis in other fruit species, the understanding of the molecular mechanism of anthocyanin biosynthesis in red-skinned pear has also advanced recently. The important structural genes phenylalanine ammonialyase (PAL), chalcone synthase (CHS), chalcone isomerase (CHI), flavanone 3-hydroxylase (F3H), dihydroflavonol-4-reductase (DFR), anthocyanidin synthase/leucoanthocyanidin dioxygenase (ANS/LDOX), UDP-glucose: flavonoid-3-O-glucosyltransferase (UFGT), and transcription factors MYB10, bHLH, and WD40 have been successfully isolated in pears (Fischer et al., 2007;Feng et al., 2010;Zhang et al., 2011;Yang et al., 2013). Previous studies reported that MYB10 plays an important role in anthocyanin biosynthesis of apples (Espley et al., 2007), strawberry (Medina-Puche et al., 2014), nectarine (Rahim et al., 2014), and other fruit trees, while in red-skinned pear, anthocyanin biosynthesis is regulated by a R2R3-MYB transcription factor PyMYB10 (Feng et al., 2010;Zhang et al., 2011;Yu et al., 2012). However, Pierantoni et al. (2010) reported that PcMYB10 is not directly responsible for red vs. yellow color in "Max Red Bartlett" and "Williams" pear varieties by map position, as the mutation underlying this difference maps to a different region of the pear genome. Our study on a pair of red/green pears also revealed that structural genes PAL, CHS, CHI, F3H, DFR, ANS/LDOX, and UFGT did not have mutations but were more highly expressed in the red pear. MYB10 was only significantly more expressed at the early stage, while the expression levels of bHLH and WD40 were higher at a later stage (Yang et al., 2013). These data indicated that the high expression of structural genes in the anthocyanin biosynthesis pathway led to the red skinned pear, and MYB10 does not appear to be the key transcription factor regulating the biosynthesis of anthocyanin and determining the red/green color mutant (Yang et al., 2013). In addition, we reported differential regulation mechanisms of anthocyanin biosynthesis and coloration pattern between occidental and oriental pears by the different co-expression of MYB10 and bHLH33 and expressions of WD40 (Yang et al., 2015). This all indicates that the molecular mechanisms of coloration can be different even between different varieties of the same fruit. So far, the red color trait in pear is comparatively more complex than expected, and the regulatory molecular mechanism has not been consistently concluded. Hence, the underlying genetics of anthocyanin biosynthesis and the metabolic pathway still need further exploration.
"Starkrimson" (Pyrus communis) is called "Early Red Doyenne Du Comice" in China, and was introduced from England to China by Changli Institute of Pomology (Hebei Academy of Agriculture and Forestry Sciences) in 2001. The fruit of "Starkrimson" has purplish red skin with the color covering the whole fruit, from the very beginning of fruit setting to maturation. The red coloration of fruit is not affected by growing conditions and show stable red color in different orchards distributed in different provinces. The green mutant of "Starkrimson" was obtained by mutagenesis with Co 60 -γ ray, and were stabilized for 5 years at very early stages of fruit development and continued into the harvesting stages (Wu et al., 2013b;Yang et al., 2013). So the color difference of the "Starkrimson" and variant strain is controlled by genetic factors, and not environmental conditions (Yang et al., 2013). The artificial induction of skin color variations offer a good opportunity to elucidate the role of genes controlling and regulating anthocyanin biosynthesis in pear fruit. In our previous study, we found the candidate gene PyMADS18, which might be involved in regulating anthocyanin biosynthesis, by screening the cDNA libraries of a pair of red and green pears (Wu et al., 2013b). We also cloned seven anthocyanin biosynthesis genes in red/green fruit skin mutant strains and detected no sequence differences between the color mutants, which indicated that the skin color change was not caused by a mutation of any of these genes. However, their expression levels differed, leading us to conclude that unknown genes might play an important role in formation of the red color, and new sequencing technologies provide an effective way of studying differentially expressed genes at the genome level.
High-throughput sequencing technologies can generate large amounts of sequence data cheaply and quickly, and have been applied widely to transcriptome analysis of plants and animals. The transcriptome is the complete set of expressed RNA transcripts during a developmental stage or in response to a particular physiological condition. It provides valuable information for identifying differentially expressed genes, and its impact on modern plant breeding gives it broad use and application, such as in M. sprengeri , pummelo (Citrus grandis) (Guo et al., 2015), grapevine (Vitis labrusca × V. vinifera) (Cheng et al., 2015), and peach (Zhou et al., 2014).
The aim of transcriptome analysis is to screen a series of candidate genes associated with the traits of interest, and to lay a good foundation for further identifying gene functions and carrying out genetic improvement of the traits. The first draft genome of the pear (Pyrus bretschneideri) was reported recently  and provides a good platform and reference for transcriptome analysis. Recently, the calyx abscission process of "Kuerlexiangli" pear (Pyrus sinkiangensis Yu) was reported (Qi et al., 2013) through RNA-seq analysis and identified candidate genes with highly dynamic changes in expression during the calyx abscission process. Subsequently, the analysis of surface brown spot formation in pear fruit was reported , exploring differentially expressed genes between russet and green pericarp offspring of a sand pear cv. "Qingxiang" × "Cuiguan" F1 population. Examination of a selected set of these categories revealed repressed expression of candidate genes for suberin, cutin, and wax biosynthesis in the russet pericarps. Also, the bud release following early defoliation of "Hosui" (Pyrus pyrifolia) was analyzed by RNA-seq  for gene expression in pear floral buds of completely defoliated plants after harvest. These transcriptome studies provide a comprehensive molecular biology insight via the differentially expressed genes into the target traits and screening for important candidate genes.
In this study, we sequenced the transcriptome of red-skinned "Starkrimson" and its green-skinned mutant using Illumina RNA-seq technology. A set of up-regulated and down-regulated genes between red-and green-skinned pear fruits were screened and help to reveal the molecular mechanism for the color mutation from red to green. Some important candidate genes were identified related to anthocyanin synthesis and red skinned fruit of "Starkrimson." The assembled annotated transcriptome sequences provide a valuable genomic resource to further understand the molecular basis of regulation of anthocyanin biosynthesis in pear. In addition, the results and strategy will also contribute to relevant research on molecular mechanisms of color variation in other fruit species.

Plant Materials and Samples Collection
The plant materials used in this study were the red-skinned pear "Starkrimson" and its green mutant, obtained from the orchard of Changli Institute of Pomology in Hebei Province during the 2012 growing season. The growth conditions of "Starkrimson" and its green mutant were described in our previous research (Wu et al., 2013b). The samples were collected at 40, 55, and 85 days after full bloom (DAFB). For these samples, 12 consistent fruits of each sample were randomly divided into three groups (for fruitlets at the early stage, 30 of each sample were selected). The fruit skin was peeled off and immediately frozen in liquid nitrogen and stored at −80 • C. Part of the samples were used to measure anthocyanin content, others were used to extract RNA and for transcriptome sequencing. For red-and green-skinned pears, a total of six independent libraries were sequenced, and samples were named as red 1 and green 1 at 40 DAFB; red 2 and green 2 at 55 DAFB; and red 3 and green 3 at 85 DAFB ( Table 1).

Measurement of Anthocyanin Content in Skins of Pear
Anthocyanin extraction of the samples was done using 1 g pear fruit skins in 5 mL 1% HCl-methanol (V/V) for 24 h at 4 • C with shading. After centrifugation for 20 min at 12 000 g, the upper aqueous phase was subjected to spectrophotometric quantification at 530, 620, and 650 nm using a UV-vis spectrophotometer (MAPADA UV-1800, China). The relative anthocyanin content was determined with the following formula: OD = (A 530 -A 620 )-0.1(A 650 -A 620 ) (Lee and Wicher, 1991). One unit of anthocyanin content was expressed as a change of 0.1 OD (unit × 10 3 g −1 FW). At least three independent samples from each group were used to obtain the mean anthocyanin content.

RNA Extraction and Transcriptome Sequencing Using Illumina Platform
The total RNA was extracted using Plant RNA Isolation Kit (Auto Lab), followed by RNA purification with RNeasy MiniElute Cleanup Kit (Qiagen) according to the manufacturer's instructions. RNA degradation and contamination was monitored on 1% agarose gels. RNA purity was checked using a NanoPhotometer R spectrophotometer (IMPLEN, CA, USA). RNA concentration was measured using Qubit R RNA Assay Kit in Qubit R 2.0 Fluorometer. The quality of total RNA was evaluated with an Agilent 2100 Bioanalyzer (Agilent Technologies, CA, USA). Only samples with RIN (RNA integrity number) ≥ 8 and 28S:18S RNA ≥ 1.5 were used for deep sequencing. The cDNAs were quantified with Qubit R RNA Assay Kit (Invitrogen, Foster City, CA), following the instructions of the manufacturer, with initial volume range of 0.1-4 µg. For mRNA library construction and deep sequencing, RNA samples were prepared using the TruSeq RNA Sample Preparation Kit according to the manufacturer's protocol, including polyA-mRNA purification and fragment, first strand cDNA synthesis, second strand cDNA synthesis, and repair, and adapter ligation, PCR enrichment, agarose gel purification, and library quality control. The library was sequenced using an Illumina HiSeq ™ 2000 (San Diego, CA, USA) by CapitalBio Corporation (Beijing, China).

Data Analysis
Raw reads obtained by the HT-2000 were filtered to exclude low complexity reads and reads containing adaptor sequences. The resulting clean reads were assembled with Trinity (Grabherr et al., 2011), and The Gene Index Clustering Tool (TGICL) (Pertea et al., 2003) was used to optimize the original Trinity assembly result by removing sequences that could not be extended on either end. The remaining high quality sequences (clean reads) were mapped to the assembled pear genome data  using Bowtie (Kanehisa et al., 2006). The assembled transcripts were also annotated using Blast2GO (Conesa et al., 2005) with GO (Ashburner et al., 2000) and KEGG (Kanehisa et al., 2008). The calculation of transcript expression was with the RPKM method (Reads Per kb per Million reads) (Mortazavi et al., 2008). Differentially expressed genes (DEG) analysis was performed using the method described by Audic and Claverie (1997) False discovery rate (FDR) (Benjamini and Yekutieli, "Raw reads" means the number of paired-end reads, "Spliced reads" means one read spliced-mapped to two exon reads and the rate (%) means the percentage of the spliced reads/clean reads. Red 1, red-skinned fruit at 40 DAFB; Red 2, red-skinned fruit at 55 DAFB; Red 3, red-skinned fruit at 85 DAFB; Green 1, green-skinned fruit at 40 DAFB; Green 2, green-skinned fruit at 55 DAFB; Green 3, green-skinned fruit at 85 DAFB.
2001) was used to determine the p-value thresholds in multiple testing. Finally, genes with a P ≤ 0.01 and Fold Change ≥ 2 were marked significantly different between the two libraries. In addition, the transcriptomic data supporting the results of this article are available at NCBI under BioProject with accession number PRJNA290937 with SRA Study accession number SUB1036370 (https://submit.ncbi.nlm.nih.gov/subs/biosample/).

Real-time PCR Analysis
To validate the expression patterns revealed by DEG results, 17 identified genes were analyzed using quantitative real-time PCR. RT-qPCR amplification and analysis were performed with the LightCycler 480 SYBR GREEN Master (Roche, USA), according to the manufacturer's instructions. The primers used for amplifying each gene are presented in Table 2. The raw data were analyzed with LightCycler 480 Software release version 1.5.0 (Roche), and the gene expression levels were determined with the 2 − T algorithm by normalizing to the Pyrus EFα1 (EFα1, accession number AY338250) and Pyrus TUB-b2 (TUB-b2, accession number AB239681) (Wu et al., 2012). qRT-PCR data are technical replicates with error bars, representing means ± SE (n = 3). Statistical and correlation analysis was performed with SPSS for Windows NT (release 8.0.0).

Changes of Anthocyanin Content in Skins during Fruit Development
The anthocyanin content of "Starkrimson" and the green variant were measured in our previous study (Yang et al., 2013). The anthocyanin content of the red-skinned pear "Starkrimson" was significantly higher than its green mutant, from ten-fold higher at 40 DAFB to seven-fold higher at 85 DAFB, with the highest difference between the red and the green color mutant at the early fruit development stage, that is, 40 DAFB. Anthocyanin content changes in fruit development in other red pear species show similar patterns. Wang et al. (2013) reported that the concentration of anthocyanin in "Max Red Bartlett" Pear fruits increased from 9 DAFB, peaking at 45 to 55 DAFB, and then gradually decreased through 85 DAFB to the mature fruit stage. Yang et al. (2015) reported that the anthocyanin content of the occidental pear was highest at the early stage of fruit development, and then showed a tendency to drop during fruit development and maturation. Oriental pears have lower anthocyanin content than occidental pears in general, and the anthocyanin content first increased and reached maximum values at 85 DAFB and then decreased in the later development stages in different cultivars of pears. These results indicate that anthocyanin accumulation in pear fruit skin mainly appears in early developmental stages, with fruit color gradually fading in later stages.

DEG Library Sequencing and Mapping Sequences to the Reference Transcriptome Database
High-throughput sequencing of RNA was performed to obtain a global view of gene expression difference related to the anthocyanin biosynthesis between "Starkrimson" and its green skin mutant. The three red-and three green-skinned pear samples were used for DEG library sequencing analysis. The number of raw reads for each library ranged from 26 to 46 million. The total length of clean reads ranged from 46 to 84 million and the rate of clean reads/raw reads ranged from 87.30 to 91.15% ( Table 1). The clean reads of the six libraries were mapped to the assembled pear genome reference sequence of "Dangshansuli" , with a mapping rate from 70.11 to 71.95% (Table 1). Mapping region classifications included exon, intron, intergenic, and spliced. In the six libraries, the ratios of exons were the highest, from 42.95 to 56.98%, and the ratios of introns were the lowest, from 2.19 to 5.42% ( Table 2). In addition, the Q30 percentages (a Q-score of 30 corresponds to an error rate of 1 per 1000) of all six libraries were above 95%, which indicated that the sequencing and RNA qualities were high, and the data obtained was reliable enough for further profile studies on gene expression.

Identification of Genes Showing Differential Expression between Red and Green Fruit
To compare differential expression between red and green coloration in pear fruit, we used RPKM to investigate transcript enrichment, which can eliminate the influence of differences in gene length and sequencing level. Furthermore, we employed "Total mapped position" means the numbers and the rate (%) of the mapped to the positions of clean reads mapped to the reference genome; "Exon" means the numbers and the rate (%) of clean reads mapped to exon of the reference genome; "Intron" means the numbers and the rate (%) of clean reads mapped to the intron of the reference genome; "Intergenic" means the numbers and the rate (%) of clean reads mapped to the intergenic of the reference genome.
FIGURE 1 | The numbers of DEGs between red/green skin color mutant pairs. Up-regulated (red), down-regulated (blue), and Total DEGs (green) were quantified. The results of four comparisons are shown.
IDEG6 to identify mRNAs showing statistically significant differences based on their relative abundance. We compared the samples from different colored fruits at the same developmental stage, so that three pairs of comparisons were implemented (Figure 1). Among these comparisons, we found that the most differentially expressed transcripts were between red 1 vs. green 1, with 1032 up-regulated unigenes and 1226 down-regulated unigenes. Simultaneously, anthocyanin content reached peak values at the early stages of fruit development. Fewer unigenes were observed between red 2 and green 2, in which only 835 unigenes were identified. Of these genes, 528 were up-regulated and 307 were down-regulated. In general, compared with green color fruit, a total of 2230 unigenes were significantly differentially expressed. Among all unigenes, 1680 up-regulated unigenes and 550 down-regulated unigenes were identified in red skinned fruits. The plots of unigenes between red and green revealed unigenes with both fold change and significance (Figure 2).

GO Classification and the Enrichment Analysis of Differentially Expressed Genes (DEGs)
Gene Ontology (GO) is an international standardized gene function classification system that describes properties of genes and their products in any organism. In order to create a profile of gene expression in red-skinned pear, we used WEGO (Web Gene Ontology Annotation Plot) for gene annotation analysis. To screen the candidate genes of DEGs from the transcripts of functional annotation, we analyzed the 4886 differential transcripts at three developmental stages from a total of 174,810 transcripts in the red/green color mutant pair. The DEGs were categorized into 28 functional groups on the basis of their biological processes (Figure 3). The major subcategories were as follows: Seven subcategories for cellular location; eight subcategories for molecular function, and 13 subcategories for biological process. The DEGs in "binding, " "catalytic activity, " "metabolic process, " and "cellular process, " "cell, " "cell part, " and "pigmentation" played important roles during the pigment metabolic process. The results provided a comprehensive view for screening candidate genes related to anthocyanin biosynthesis and metabolic process of red/green color pear fruit. Taking into account that anthocyanin content of red color fruit at different developmental stages is higher than in the green mutant, we searched for DEGs at the FIGURE 2 | Comparison of transcript expression between red and green fruit. The abundance of each gene was normalized as reads per kb per Million reads (RPKM). The differentially expressed genes are shown in red and green, while black indicates genes that were not differentially expressed (not DEGs) between red and green fruit.

Kyoto Encyclopedia of Genes and Genomes (KEGG) Enrichment Analysis of DEGs
KEGG analysis provides information and further understanding on how "Starkrimson" pear regulates its biological functions and synthesizes secondary metabolites, including anthocyanin, at the molecular level. Usually, unigenes in the same pathway cooperate with each other. In this study, there were 42,684 unigenes that mapped to 248 KEGG pathways using Blast X homology search. Furthermore, there were 783 DEGs mapped in 101 KEGG pathways, including metabolism (595 DEGs), genetic information processing (68 DEGs), environmental information processing (58 DEGs), cellular processes (36 DEGs), and human diseases (26 DEGs) ( Table 4). Among them, 12 pathways were significantly enriched with more than 15 DEGs (Q ≤ 0.05) ( Table 5). In the metabolism  category, the most abundant DEGs were found in the amino acid metabolism sub-category (PATH:ko00270, PATH:ko00330, PATH:ko00250; PATH:ko00071), followed by carbohydrate metabolism (PATH:ko00010), lipid metabolism (PATH:ko00561; BR:ko01004), and metabolism of other amino acids (PATH:ko00982; PATH:ko00980) ( Table 4). Five hundred ninety five DEGs were related to catalysis of metabolism processes or generation of energy for primary and secondary metabolite production. We also found 3 DEGs in the secondary metabolite pathway of flavonoid biosynthesis (PATH:ko00941): p-coumarate 3-hydroxylase (C3H), Anthocyanidin Reductase (ANR), and leucoanthocyanidin reductase (LAR). In addition, the genetic information processing category included transcription (21 DEGs) (BR:ko03000) ( Table 6), folding, sorting, and degradation (31 DEGs), replication, and repair (16 DEGs). The DEGs in this category mainly function in ensuring correct transcription and translation processes, and the transcription factors indirectly involved in fruit development transcription regulation and chromosome (BR:ko03036) determine metabolism, cellular processes, and environmental information processing.

Genes Involved in Flavonoid Biosynthesis, Metabolism, and Transportation
Flavonoid biosynthesis is a dynamic and complex processes catalyzed by a series of enzymes. Flavonoids are a diverse group of plant secondary metabolites with various biological functions that play important roles during plant development. Proanthocyanidins (PAs, also known as condensed tannins) are components of metabolites synthesized through the general flavonoid biosynthesis pathway (Figure 4). In a previous study, it was reported that leucoanthocyanidin reductase (LAR), anthocyanidin synthase (ANS; also called leucoanthocyanidin dioxygenase, LDOX), and anthocyanidin reductase (ANR; in Arabidopsis, the product of the BANYULS gene) were the three principal enzymes for flavan-3-ols biosynthesis. The synthesis of PAs and anthocyanins share common steps leading to flavan-3,4-diols (such as leucoanthocyanidin), which can be converted to catechin (2,3-trans-flavan-3-ol) by LAR  or to anthocyanidin by ANS . Anthocyanidin then either serves as the substrate for the synthesis of epicatechin (2,3-cis-flavan-3-ol) by ANR (Xie et al., 2003) or can otherwise be converted to anthocyanin by glycosylation (Schijlen et al., 2004). Recently, it was proven that transgenic tobacco overexpressing TcLAR had decreased amounts of anthocyanidins and increased PAs. Overexpressing TcLAR in an Arabidopsis ldox mutant also resulted in elevated synthesis of not only catechin but also epicatechin (Liu et al., 2013a,b). In strawberry, it was demonstrated that redirection of the anthocyanin pathway to flavan-3-ols was observed by the down-regulation of an anthocyanin glucosyltransferase in ripening strawberry fruit (Griesser et al., 2008), while the downregulation of anthocyanin reductase (ANR) induced a redirection of the proanthocyanidin pathway, leading to premature, and ectopic anthocyanin biosynthesis via enzymatic glycosylation as the alternative pathway (Fischer et al., 2014). In our study, ANR showed relatively high expression levels in green color mutant of "Starkrimson" pear, leading to unstable anthocyanidin through ANS catalysis. This could not be converted to stable colored anthocyanin by 3,5-glycoside flavonoid transferase (UFGT), and thus anthocyanin did not accumulate in the pear fruit. Therefore, high ANR expression might be a main reason leading to the skin color change from red to green in this mutant. However, further study is needed to verify gene function and reveal the regulatory mechanism underlying the phenomenon. In addition, we detected significantly expressed p-coumarate 3-hydroxylase (C3H), which is in the flavonoid biosynthesis pathway. In a previous study, C3H was found to be a key gene playing an important role in lignin biosynthesis metabolism. Dardick et al. (2010) reported that p-coumarate 3-hydroxylase is involved in stone formation in peach fruit by microarray studies using a developmental series from young fruits. Recently, Xue et al. (2014) reported synthesis and codon-optimiztion of an Arabidopsis thaliana ref8 gene encoding a C3H for enhanced expression in Synechocystis. This heterologous pathway enabled Synechocystis to produce caffeic acid. As this is the first report of a C3H being related to flavonoid metabolism in pear, the gene function still needs further study.
GST is a family of multifunctional enzymes catalyzing detoxification reactions, and are also involved in endogenous FIGURE 4 | The flavonoid biosynthesis pathway leading to PA production in pear fruit. PAL, phenylalanine ammonialyase; C4H, cinnamate 4-hydroxylase; 4CL, 4-coumarate coenzyme A ligase; CHS, chalcone synthase, CHI, chalcone isomerase; F3H, flavanone 3-hydroxylase; FLS, flavonol synthase; DFR, dihydroflavonol-4-reductase; LAR, leucoanthocyanidin reductase; ANS/LDOX, anthocyanidin synthase/leucoanthocyanidin dioxygenase; ANR, anthocyanidin reductase; OMT, O-methyltransferases; UFGT, UDP-glucose: flavonoid-3-O-glucosyltransferase; RT, rhamnosyltransferase; GST, glutathione S-transferase. The blue arrow indicates a process detected between red and green skinned pear, the black arrow was clarified in other plant. In addition, the rhombus shape indicates anthocyanin accumulation, the rectangle shape proanthocyanidin biosynthesis, and the ellipse shape flavonoid biosynthesis. metabolism, including functioning as GSH-dependent isomerases, non-catalytically acting as flavonoid-binding proteins, stress signaling proteins, and regulators of apoptosis (Jain et al., 2010). In previous studies, there have been many related reports showing GST encoded enzymes involved with anthocyanin transportation from cytoplasm to vacuole. Anthocyanin synthesis is in the cytoplasm, and is then transported to the vacuole and stored, displaying different colors with different ion concentrations and pH conditions in the vacuole (Tanaka et al., 2008). Until now, it has been confirmed that some members of GST are involved in anthocyanin transportation in Arabidopsis (Kitamura et al., 2004), petunia (Alfenito et al., 1998), carnation (Larsen et al., 2003), corn (Marrs et al., 1995), and grapes (Conn et al., 2008). In this study, we detected GST (Pbr012649.1) to be up-regulated in red-skinned "Starkrimson." However, little is known about GST genes related to anthocyanin transportation in pear, the full transport mechanisms still need further investigation.

Genes Encoding Transcription Factors
The regulation of gene expression at the transcription level has a profound role in the control of many biological processes. Transcription factors are the key switches for secondary metabolite gene regulation. Between the red and green color mutant, 21 DEGs were identified as transcription factors from the 595 DEGs annotated into metabolite sub-category of the KEGG pathway in this study ( Table 6). These were annotated as MYB, AP2, WRKY, and MADS transcription factors. Among the group of transcription factors, we identified eight genes belonging to the MYB family of transcription factors. Previous studies have mostly focused on the superfamily of MYB TFs, which has been shown to be involved in the control of many biological processes, including anthocyanin biosynthesis (Takos et al., 2006), for example, a R2R3 MYB transcription factor associated with regulation of the anthocyanin biosynthetic pathway in Rosaceae (Lin- . MYB90/PAP2, MYB6, MYB10, MYB1, MYBA, and MYB19 have been reported to be involved in the regulation of anthocyanin synthesis (Borevitz et al., 2000;Takos et al., 2006;Ban et al., 2007;Gonzalez et al., 2008;Yamagishi et al., 2010;Telias et al., 2011;Zhang et al., 2012). In apple (Malus × domestica), expression of MYB10 showed a strong correlation with anthocyanin content during fruit development of redfleshed apple "Red Field" (Espley et al., 2007). Two more apple TFs, MYB1 and MYBA, were also reported to regulate genes in the anthocyanin pathway in red-skinned fruit (Ban et al., 2007). Both MYB1 and MYBA share identical sequences, while MYB10 and MYB1 genes are located at very similar positions on linkage group 9 of the apple genetic map (Chagné et al., 2008). In strawberry (Fragaria × ananassa), FaMYB1 plays a key role in down-regulating the biosynthesis of anthocyanins and flavonols (Aharoni et al., 2001). DkMYB4 is involved in proanthocyanidin biosynthesis in persimmon fruit (Akagi et al., 2009)

, while
AtMYBL2 encodes a protein with a single MYB domain that acts as a negative regulator of anthocyanin biosynthesis in Arabidopsis (Matsui et al., 2008). Recently, Guo et al. (2015) identified a number of MYB transcription factors that may be involved in the carotenoid regulation in orange-pericarp mutants and wild type in pummelo (Citrus grandis), via transcriptomic analysis. These results are examples of MYB genes' frequent involvement in regulating anthocyanin biosynthesis. In this study, we detected eight up-regulated MYB TFs in red skinned "Starkrimson, " but did not include MYB10, corresponding to the previous deduction that PcMYB10 is not directly responsible for red vs. yellow color in "Max Red Bartlett" and "Williams" pear (Pierantoni et al., 2010), and other key genes control the color mutant of "Starkrimson" pear and the green variant (Yang et al., 2015). In addition, anthocyanin accumulation is related to many biological processes, and as MYB TFs is a large regulatory gene family, some MYB genes might be involved in anthocyanin accumulation through regulation of related biological process. So, we should further verify these new MYB candidate genes to find one key gene controlling red coloration of pear. In addition, the promoter region (3000 bp upstream of ATG) of LAR and ANR genes were analyzed, and were found to have MYB-binding cis-motifs in the promoter region of both genes (Additional File 1). So we speculated that MYB TFs could bind to the promoter of LAR or ANR and regulate their expression, affecting fruit coloration in pears.
In previous studies on anthocyanin biosynthesis in the flavonoid pathway, genes have been found to be coordinately modulated by a conserved MYB-bHLH-WD40 (MBW) regulatory complex. In Arabidopsis, bHLH EGL3 and GL3 function in both the biosynthesis of anthocyanins and the formation of trichomes and root hairs (Ramsay and Glover, 2005). Xie et al. (2012) demonstrated that MdbHLH3 regulates LT-induced anthocyanin accumulation and fruit coloration in apple. MdbHLH300 is down-regulated in fruits grown in a hot climate compared with a temperate climate (Lin- . In this study, TFs bHLH, Pbr037113.1, and Pbr041849.1 were up-regulated and down-regulated, respectively, indicating their different regulation roles in anthocyanin biosynthesis. We also detected that a MADS-Box protein encoding the gene MADS15 (Pbr002427.2) was up-regulated in the early stage of fruit development. Wu et al. (2013b) screened a new transcription factor PyMADS18, which is likely to be involved in anthocyanin accumulation and regulation of anthocyanin synthesis in early pear fruit development. In this study, the gene (Pbr013927.1) was found to belong to AP2/ERF domain-containing transcription factors, and was up-regulated in red-skinned fruit at all three measured developmental stages. Four genes (Pbr030451.1, Pbr016185.1, Pbr022708.1, and Pbr037846.1) were up-regulated at both early and later stages of fruit development. In addition, we also identified four significantly differentially expressed WRKY transcription factors, WRKY9/33 (Pbr015939.1), WRKY17 (Pbr009294.1), WRKY51 (Pbr026903.1), and WRKY40 (Pbr026903.1), however, the function of these genes in the coloration of pear fruit still needs further study.

Real-time qPCR Validation of Differentially Expressed Transcripts from Transcriptome Analysis
To confirm the accuracy and reproducibility of the transcriptome analysis results, 17 genes having different expression patterns were used for real-time qPCR verification and their correlation evaluated. The results showed that although the exact fold changes for the selected genes at three stages varied between digital gene expression and qRT-PCR analysis, the trend of gene expression changes detected by the two different approaches were largely consistent (Figure 5). Pearson's correlation coefficients showed that the digital transcript abundance measurements and qRT-PCR data were highly correlated (Table 7), with R 2 -values ranging from 0.639 (Pbr027485.1) to 0.998 (Pbr004885.1), which was in agreement with previous reports (Pan et al., 2012).

R-ATGTTACTGTGAGGGACGTCTGG
In particular, the expression levels of screened candidate genes were in accordance with the transcriptome results. We verified candidate TFs such as AP2 (Pbr016185.1) and WARK (Pbr004885.1), which had similar expression trends with transcriptome analysis results, with higher expression in red-skinned "Starkrimson" than its green variant (Figure 5). This indicated that these genes were related to the anthocyanin biosynthesis and regulated the formation of red skinned pear. Meanwhile, LAR (Pbr027485.1) and ANR (Pbr000218.1) were more highly expressed in the green variant than red "Starkrimson" (Figure 5), which promote the process of PA pathway and contributed to the formation of green skinned pear. In addition, the trends of the genes Pbr008092.1 and Pbr021636.1 were consistent with anthocyanin biosynthesis and accumulation in fruit developmental stages, while the genes of Pbr019417.1 and Pbr041849.1 showed opposing changes. qRT-PCR further demonstrated that genes related to metabolism of cofactors and vitamins (Pbr019417.1), amino acid metabolism (Pbr036011.1, Pbr041919.1, Pbr012649.1), carbohydrate metabolism (Pbr008092.1), genetic information processing (Pbr041849.1), and other regulated gene (Pbr036167.1) showed significant difference between red-skinned pear "Starkrimson" and its green mutant, indicating that anthocyanin accumulation in pear fruit is related to many biological processes.

Conclusions
The present results demonstrate the usefulness of the digital transcript abundance measurement approach for identifying critical DEGs in the red-skinned pear "Starkrimson" and its green mutant. A list of candidate genes for functional studies involving anthocyanin/ proanthocyanidin biosynthesis and regulations was generated, among which LAR and ANR appeared to play a key role in promoting the proanthocyanin biosynthesis pathway. Their expression led to low anthocyanin accumulation in fruit skin, showing that LAR and ANR are associated with the red/green skin color mutant of pear. Furthermore, several transcription factors, MYB, AP2, WRKY, and MADS, were identified as differentially expressed in red/green color mutant, indicating that multiple transcription factors are involved in the regulation of anthocyanin/proanthocyanidin biosynthesis in pear fruit and lead to different fruit colors. The qRT-PCR results also indicated that the screening candidate genes AP2 and WARK were involved in red skinned pear formation, and LAR and ANR were related to the development of green skin. Further studies should concentrate on functional characterization of these genes. This study could lead to better understanding of the global network of differentially expressed genes for red or green coloration in pear fruits.