Transcriptomic and metabolomic analysis provides insights into anthocyanin and procyanidin accumulation in pear

Background Pear is one of the most important fruit crops worldwide. Anthocyanins and procyanidins (PAs) are important secondary metabolites that affect the appearance and nutritive quality of pear. However, few studies have focused on the molecular mechanism underlying anthocyanin and PA accumulation in pear. Results We conducted metabolome and transcriptome analyses to identify candidate genes involved in anthocyanin and PA accumulation in young fruits of the pear cultivar ‘Clapp Favorite’ (CF) and its red mutation cultivar ‘Red Clapp Favorite’ (RCF). Gene–metabolite correlation analyses revealed a ‘core set’ of 20 genes that were strongly correlated with 10 anthocyanin and seven PA metabolites. Of these, PcGSTF12 was confirmed to be involved in anthocyanin and PA accumulation by complementation of the tt19–7 Arabidopsis mutant. Interestingly, PcGSTF12 was found to be responsible for the accumulation of procyanidin A3, but not petunidin 3, 5-diglucoside, opposite to the function of AtGSTs in Arabidopsis. Transformation with PcGSTF12 greatly promoted or repressed genes involved in anthocyanin and PA biosynthesis, regulation, and transport. Electrophoretic mobility shift and luciferase reporter assays confirmed positive regulation of PcGSTF12 by PcMYB114. Conclusion These findings identify a core set of genes for anthocyanin and PA accumulation in pear. Of these, PcGSTF12, was confirmed to be involved in anthocyanin and PA accumulation. Our results also identified an important anthocyanin and PA regulation node comprising two core genes, PcGSTF12 and PcMYB114. These results provide novel insights into anthocyanin and PA accumulation in pear and represent a valuable data set to guide future functional studies and pear breeding.


Background
Pear is an important fruit for human consumption, and its total global production is ranked third after grape and apple [1]. Pear is cultivated commercially in 76 countries or regions worldwide [2], among which China is the world's leading pear producer. In 2017, China produced 16.4 million tons (Mt) of pear fruits, accounting for 68% of global pear production (24.2 Mt) (FAOSTAT, 2017).
The anthocyanin and PA biosynthesis pathways have been well characterized in plants. Anthocyanin and PA are initially biosynthesized from phenylalanine, and share most steps in the flavonoid biosynthetic pathway. They are biosynthesized in the cytosol by enzymes including phenylalanine ammonia lyase (PAL), chalcone isomerase (CHI), chalcone synthase (CHS), flavonoid 3hydroxylase (F3H), and dihydroflavonol reductase (DFR) [6]. Leucoanthocyanidins and anthocyanidins are two important branch points between the anthocyanin and PA biosynthesis pathways. Downstream of these branch points, anthocyanins are synthesized by anthocyanidin synthase (ANS) and UDP-glucose flavonoid 3-Oglucosyl transferase (UFGT), and PAs are synthesized by leucoanthocyanidin reductase (LAR) and anthocyanidin reductase (ANR) [7,8]. O-methyltransferase (OMT) and glycosyltransferase (GT) are responsible for the elaboration of diverse anthocyanins and PAs [9,10]. After anthocyanins and PAs are synthesized in the cytosol, they are transported to their final destination, the vacuole. Some glutathione S-transferases (GSTs) and multidrug and toxic compound extrusion proteins (MATEs) are believed to function as anthocyanin and PA carrier proteins to sequester them into vacuoles [11].
The molecular mechanism underlying anthocyanin and PA accumulation has been extensively studied in numerous plants. Many anthocyanin and PA structural genes and their upstream regulators have been identified and characterized. Of these, R2R3-MYB TFs play key roles in controlling anthocyanin and PA accumulation by acting together with bHLH and WD40 proteins to regulate structural genes. Genes encoding R2R3-MYB TFs that contribute to anthocyanin accumulation include MdMYB1 in apple, PyMYB10 and PyMYB114 in pear, and VvMYBA1 in grape [12][13][14][15]. Some R2R3-MYB TFs regulate both anthocyanins and PAs, including VvMYB5a and VvMYB5b in grape, MdMYB9 and MdMYB11 in apple, PbMYB10b and PbMYB9 in pear, and PpMYB18 in peach [4,[16][17][18][19]. Other R2R3-MYB TFs are related only to the regulation of PAs, including VvMYBPA1 in grape and PpMYBPA1 in peach [16,20]. In addition, other TFs such as AUX and ERF also regulate anthocyanin or PA biosynthesis by directly or indirectly interacting with R2R3-MYB TFs and structural genes [21,22].
Recent technical advancements in transcriptome and metabolome analyses have provided effective ways to identify new genes and metabolites, and to elucidate complex secondary metabolic bioprocesses in plants. In fig (Ficus carica L.), integrated transcriptome and metabolome analyses have revealed genes in flavonoid and anthocyanin pathways that show differential expression between purple-and green-skinned cultivars [23]. Another study using combined transcriptome and metabolome datasets successfully constructed expression-anthocyanin metabolite networks in potato [24]. Recently, seven flavonoid metabolites and six genes were identified as candidates associated with the pigmentation of the red-fleshed and green-fleshed cultivars of Actinidia arguta [25].
Red pears have high nutritional and economic value because they are rich in anthocyanins and PAs. Therefore, studies on the regulation of anthocyanins and PAs are of great interest for the improvement of anthocyanin and PA production in pears. Considering the large number of anthocyanins and PAs, the molecular mechanisms of their biosynthesis and embellishment in pear might be more complex than expected. Bud mutation is an important method of selecting new red pear varieties. Red mutations are ideal materials to study the molecular mechanism of anthocyanin and PA accumulation because of their highly similar genetic backgrounds [13]. In the present study, we carried out metabolome and transcriptome analyses to identify candidate genes involved in anthocyanin and PA accumulation in pear using the young fruits of 'CF' and its red mutant 'RCF'. Correlation analyses between differentially expressed genes (DEGs) and anthocyanins/PAs revealed 203 candidate genes for the accumulation of 10 anthocyanins and seven PAs in pear. Of these, 20 genes were strongly correlated with all 10 anthocyanins and seven PAs. Thus, they seemed to be the core candidate genes related to anthocyanin and PA accumulation in pear. The GST gene PcGSTF12 was correlated with most anthocyanin and PA metabolites. PcGSTF12 was confirmed to play an important role in anthocyanin and PA accumulation in pear by functional complementation analyses. In addition, PcGSTF12 was found to be directly and positively regulated by PcMYB114, a well-known TF regulating anthocyanin accumulation in pear. These results greatly extend our knowledge of the molecular mechanism of anthocyanin/ PA accumulation in pear.

Results
Anthocyanin and PA profiles of 'CF' and its red mutant 'RCF' Except for skin color, no significant morphological differences were observed between the fruits of 'CF' and its red mutant 'RCF'. The young fruits of both 'CF' and 'RCF' initially had a deep green appearance. The color difference between 'CF' and 'RCF' became visible from about 5 days after full bloom (DAFB). The 'RCF' fruits quickly turned dark red, and retained their strong color until maturity. In contrast, the fruits of 'CF' only developed a slight red blush on the sun-exposed surface (Fig. 1a).

Transcriptome analysis
The RNA-seq process yielded 95.6 G clean bases and 637 million clean reads. The mean number of clean reads per sample was 53 million. Of the clean reads, 93.54% were mapped in total, and 90.72% were mapped uniquely against the improved apple reference genome sequence. In total, 14,514 genes were expressed with FPKM ≥10 (Supplementary Table S2).
The 203 transcripts were annotated with descriptions from the SwissProt and NR databases. Six transcripts have been functionally characterized to play roles in anthocyanin accumulation in pear previously: PcMYB10, PcMYB114, PcCHS, PcCHI, PcF3H, and PcANS (Supplementary Table S7). The rest were newly identified as candidate genes involved in anthocyanin and PA accumulation in pear. The 203 transcripts were grouped into two clusters (I-II) (Supplementary Table S7). Genes in cluster I were strongly correlated with anthocyanins. Cluster I comprised 183 genes (90.1%). Of these, 147 genes were correlated with a single anthocyanin: 142 genes were correlated with kuromanin chloride, three genes were correlated with cyanidin, one gene was correlated with malvidin-3-galactoside chloride, and one gene was correlated with pelargonin. The remaining genes in cluster I were closely correlated with two or more anthocyanins: 31 genes were commonly correlated with kuromanin chloride and pelargonin, two genes were correlated with pelargonin and cyanidin, two genes were correlated with malvidin-3-galactoside chloride and oenin, and one gene was correlated with cyanidin 3rutinoside, oenin, and cyanidin. Cluster II contained 20 genes (9.9%) that were strongly correlated with both anthocyanins and PAs. Of these genes, two phenylpropanoid structural genes (encoding 4CL1 and 4CL2), six flavonoid structural genes (encoding CHS, 3 CHIs, F3H, and ANS), six TF genes (encoding bZIP1, MYB3, MYB86, MYB111, MYB114, and KNAT1), two phytohormone signal transduction genes (encoding IAA13 and ERF003), two DNA-directed RNA polymerase genes (encoding rpoB and Rpb1) and one GST transporter gene (encoding GSTF12) were positively correlated with anthocyanins and PAs. One WRKY TF gene, WRKY28, was negatively correlated with anthocyanins and PAs (Supplementary Table S7, S8). Each gene in cluster II was strongly correlated with many metabolites. We found that these 20 genes were strongly correlated with all 17 anthocyanin and PA metabolites (Fig. 3a). Thus, they were considered to represent the core genes for anthocyanin and PA accumulation in pear. Of these, PcRPB1 (PCP004386) was correlated with the fewest metabolites: one anthocyanin and three PAs; and PcGSTF12 (PCP025171) was correlated with the most metabolites: seven anthocyanins and seven PAs (Fig. 3b).

qPCR analysis of DEGs related to anthocyanin and PA accumulation
To validate the RNA-seq data, we conducted qPCR analyses of 10 of the anthocyanin and/or PA candidate genes: PcCHI, PcC1, PcMYB114, PcHB7, PcGAI1, PcCHS, PcGSTF12, PcANS, PcHB12, and PcMYB10 (for  Table S9). The transcript profiles of all selected genes were very similar to those detected from the RNA-seq data (Fig. 4). The results showed that PcGSTF12 was most up-regulated in comparison groups 'RCF1 vs. CF1' and 'RCF2 vs. CF2'. This result was highly consistent with the results of RNA-seq, and provided further evidence for the crucial role of PcGSTF12 in anthocyanin and PA accumulation in pear. Thus, we conducted further analyses to confirm the function of PcGSTF12.

PcGSTF12-mediated anthocyanin and PA accumulation in pear
Our combined metabolite and transcriptomic analyses revealed a core set of genes closely correlated with pear anthocyanins and PAs, which strongly suggested that Fig. 3 Connection network between core genes and anthocyanin and procyanidin (PA) metabolites. a Networks between 20 core genes and 10 anthocyanins and seven PAs. b Networks between PcGSTF12 (PCP025171) and seven anthocyanins and seven PAs. PcGSTF12 is shown in red font they play key roles in anthocyanin and PA accumulation in pear. To test this, we focused on the most upregulated gene among the core set of anthocyanin and PA candidate genes, PcGSTF12, for functional analysis.

Functional analysis of PcGSTF12
The phylogenetic analysis showed that PcGSTF12 is a homolog of FvRAP in strawberry, Riant2 in peach, and MdGST in apple, all of which are in the phi subfamily [26] (Fig. 5a). Members of the phi subfamily are anthocyanin transporters. To test the potential role of PcGSTF12 in anthocyanin accumulation, 35S:: PcGSTF12 was transformed into the Arabidopsis mutant tt19-7 (for primers, see Supplementary Table S9). The tt19-7 plants showed a green hypocotyl phenotype, while the tt19-7-OE transgenic plants showed the red hypocotyl phenotype, like that of the wild type (WT) (Fig. 5b). However, the brown color of seed coats was not rescued in the tt19-7-OE lines (Fig. 5b). This result was consistent with the fresh seed phenotype obtained by transferring 35S:: RAP-RFP into Arabidopsis tt19-7 [26].
Next, we analyzed RNA-seq data to identify which genes were affected by PcGSTF12 in the seedlings of tt19-7-OE vs. tt19-7. In total, we found 28 strongly affected genes, which encoded proteins involved in anthocyanin and PA biosynthesis, regulation, and transport (Supplementary Table S10). These results showed that PcGSTF12 might not only be an anthocyanin and PA transporter, but may also participate in many other steps of anthocyanin and PA accumulation.

Upstream regulation of PcGST12
Correlation analyses showed that the transcript level of PcGSTF12 was significantly correlated with that of PcMYB114 (Fig. 5d). Further, several MYB-binding sites were found in the PcGSTF12 promoter, indicating that PcGSTF12 might be directly bound by, and regulated by, MYB transcription factors (Supplementary Table S11). Several R2R3-MYB genes are known to bind to MBS sites [27], and we found a MBS site within an 801-bp region upstream of the start codon. Thus, this MBS site was used in an EMSA assay (for primers, see Supplementary Table S9). The biotinylated probe was able to bind PcMYB114 protein, and the addition of a high concentration of cold probe significantly reduced the binding affinity of the biotinylated probe. To test whether PcGSTF12 could be regulated by PcMYB114, we further carried out luciferase reporter assay (for primers, see Supplementary Table S9). The relative LUC Phenotypes of tt19-7-OE seedlings and fresh seeds. c Heat map of anthocyanins and PAs of 7-day-old seedlings of tt19-7 mutant, tt19-7-OE, and WT (wild-type). tt19-7-OE refers to 35S::PcGSTF12 transgenic lines. Scale bars: 5 mm (B-C). Differentially accumulated metabolites between WT and tt19-7-OE marked by red star. d Correlation analyses between transcript levels of PcGSTF12 and PcMYB114. e EMSA assays. Probe was biotin-labeled fragment containing MBS motif. Competitor probe was non-labeled probe. Mutant probe contained two nucleotide mutations. Competitors, and mutant probes at a 10×, 25×, and 50× molar excess were present (+) or absent (−) in each reaction. Black arrows indicate increasing multiples. f Effects of PcMYB114 on promoter activities as demonstrated by luciferase reporter assays. Empty-vector was used as a control. Values are means ± SD of three independent biological replicates. Statistical significance: **P < 0.01 activity of the PcGSTF12 promoter was about 6-fold that of the empty vector control. These results showed that PcMYB114 could directly bind to the MBS site in the PcGSTF12 promoter (Fig. 5e) and positively regulate its activity (Fig. 5f).
We also found many cis-acting elements involved in auxin-, ethylene-, and gibberellin-signaling in the PcGSTF12 promoter. This indicated that PcGSTF12 might be the common downstream target of R2R3-MYBs, and auxin, ethylene, and gibberellin signals that regulate the anthocyanin and PA pathways (Supplementary Table S11). Together, these results provide new clues about PcGSTF12-mediated anthocyanin and PA accumulation in pear.

Discussion
We present the first genome-wide examination of anthocyanins, PAs, and the gene expression profiles of pear using young fruits of cv. 'CF' and its red mutant 'RCF'. Through combined transcriptomic and metabolic analyses, we found a core set of 20 candidate genes related to anthocyanin and PA accumulation in pear. These findings increase our understanding of the molecular mechanism of anthocyanin and PA accumulation in pear, especially during the early stage of fruit development.
The core set of candidate genes for pear anthocyanin and PA accumulation includes six flavonoid structural genes: PcCHI (PCP027877, PCP021141, and PCP039844), PcCHS (PCP023048), PcF3H (PCP013732), and PcANS (PCP027029). CHI, CHS, and F3H are common to both the anthocyanin and PA biosynthetic pathways [6]. ANS can catalyze the conversion of (+)-catechin to cyanidin and a procyanidin [28]. Functional analyses have confirmed the roles of PcCHI, PcCHS, and PcF3H in anthocyanin accumulation [29]. The results of those studies and our study provide evidence that the patterns of anthocyanin and PA accumulation are conserved among different plants.
In plants, TFs play important roles in flavonoid regulation. The R2R3-MYBs make up one of the largest TF families [30], and most of them play important roles in flavonoid accumulation [31]. PyMYB10 was the first R2R3-MYB TF identified to be involved in anthocyanin accumulation in pear [13]. Functional analyses have confirmed the roles of PyMYB10.1 and PyMYB114 in the regulation of anthocyanin accumulation in pear [14,32], and the roles of PbMYB10b and PbMYB9 in both anthocyanin and PA accumulation in pear [4]. In this study, the core set of candidate genes for anthocyanin and PA accumulation included four R2R3-MYB genes: PcMYB3, PcMYB86, PcMYB111 and PcMYB114. Their transcript levels were strongly positively correlated with anthocyanins and PAs. Although MYB114 was already known to function in anthocyanin accumulation in pear, it was unknown which particular metabolites were affected. Our results provide evidence that PcMYB114 functions in the accumulation of six anthocyanin metabolites: pelargonidin 3-Glu, malvidin-3-galactoside chloride, cyanidin 3-O-malonylhexoside, oenin, delphinidin and peonidin O-hexoside (Supplementary Table S8). We also detected positive correlations between PcMYB114 and procyanidin A, procyanidin B1, and procyanidin B3 (Supplementary Table S8), indicating that PcMYB114 also functions in PA accumulation in pear. Except for PcMYB114, the other MYB genes PcMYB3, PcMYB86, and PcMYB111 are newly identified as candidates involved in anthocyanin and PA accumulation in pear. Interestingly, these R2R3-MYBs were correlated with different anthocyanins and PAs, suggesting that they have undergone sub-specialization to play different and specific roles in anthocyanin and PA accumulation in pear.
Other TFs are also involved in anthocyanin and PA accumulation in pear. The HD family of TFs is unique to plants, and its members are proposed to play key roles in developmental processes such as root development, plant cell differentiation, fruit ripening, and leaf and flower senescence [33][34][35]. Members of the HD-Zip I and HD-Zip IV TF subfamilies also play key roles in anthocyanin accumulation. ANTHOCYANINLESS2 (ANL2) was the first HD-Zip IV gene found to be involved in the tissue-specific accumulation of anthocyanins. In Arabidopsis, ANL2 affects anthocyanin accumulation in subepidermal tissue on the abaxial side of rosette leaves and in epidermal tissue on the abaxial side of the leaves [36]. Recently, two HD-Zip I genes, MdHB1 and RhHB1 have been shown to influence anthocyanin accumulation in apple and rose, respectively. Overexpression of MdHB1 led to reduced anthocyanin accumulation in apple flesh. MdHB1 was found to repress the transcription of MdDFR and MdUFGT indirectly by interacting with MdMYB10, MdbHLH3, and MdTTG1 [37]. Consistent with the results in apple, silencing of RhHB1 in rose led to higher anthocyanin levels in the petals [38]. In this study, the core set of genes involved in anthocyanin and PA accumulation in pear included an HD TF gene, PcKNAT1. KNAT1 is a member of the Class I KNOX HD gene family and is thought to play a role in meristem development and leaf morphogenesis [39]. We found that PcKNAT1 was strongly positively correlated with five anthocyanins and seven PAs, implying that its function differs from the known functions of HD-Zip I and HD-Zip IV TFs in anthocyanin accumulation. Our results suggest that the Class I KNOX HD gene family might play important roles in anthocyanin and PA accumulation; this expands our knowledge of the function of the Class I KNOX HD gene family.
In our study, a bZIP TF gene, PcbZIP1, was positively and closely correlated with six anthocyanins and five PAs. The bZIP TFs harbor a highly conserved bZIP domain [40]. They are diverse transcriptional regulators, playing crucial roles in plant development, physiological processes, and biotic/abiotic stress responses [41]. Recently, two bZIP TFs, MdHY5 and MdbZIP44, were shown to promote anthocyanin accumulation in apple [42,43]. Our results provide further evidence that bZIP TFs act as positive regulators in anthocyanin accumulation. Thus, as well as functioning in anthocyanin accumulation, bZIP TFs may also play important roles in PA accumulation.
In this study, a WRKY TF gene, WRKY28, was closely negatively correlated with one anthocyanin and four PAs. WRKY41-1 in Brassica napus and WRKY75 in Arabidopsis are known repressors of anthocyanin biosynthesis [44,45], while MdWRKY40 in apple is a positive regulator of wounding-induced anthocyanin biosynthesis [46]. Our results are consistent with findings in B. napus and Arabidopsis and opposite to those in apple. It is possible that WRKYs have evolved different functions in anthocyanin accumulation, like the R2R3-MYB TFs. For example, some R2R3-MYB TFs are positive regulators of anthocyanin accumulation [13,14], while others are negative regulators [19]. Further studies are required to elucidate the complex roles of WRKYs in PA accumulation in fruit trees and other plants.
Auxin is known to suppress anthocyanin accumulation, and has been shown to decrease the expression of anthocyanin regulatory and structural genes in apple and Arabidopsis [47,48]. A recent study showed that the auxin factor MdARF13 negatively regulates the anthocyanin pathway in apple through interacting with MdMYB10 and binding to the promoter of MdDFR. The auxin/IAA protein MdIAA121 was shown to inhibit the recruitment of MdARF13 to the MdDFR promoter and weaken the inhibitory effect of MdAFR13 on anthocyanin accumulation [21]. Contrary to auxin, ethylene enhances anthocyanin and PA accumulation in pear and apple. In pear, the ethylene-responsive factor PyERF3 enhances anthocyanin accumulation via interacting with PyMYB114 [14]. In apple, MdERF1B regulates anthocyanins and PAs through interacting with MdMYB1, MdMYB9, and MdMYB11 [22]. In this study, an ethylene-responsive gene, PcERF003, was positively correlated with six anthocyanins and seven PAs, and an auxin-responsive gene, PcIAA13, was positively correlated with two anthocyanins and four PAs. These results are consistent with the previous findings that auxin represses anthocyanin accumulation via IAA genes, and ethylene enhances anthocyanin accumulation via ERF genes. Phytohormones play important roles in young fruit development. Young pear fruits contain high level of auxin, and low level of ethylene [49]. Further functional analyses of PcERF003 and PcIAA13 may help to elucidate the links between the effects of auxin and ethylene on anthocyanin and/or PA accumulation in young developing pear fruit. Plant GSTs are encoded by a large gene family, and are soluble and abundant in the cytosol [26,50,51]. They can be divided into eight subgroups, among which the tau and phi classes play key roles in flavonoid transport [52]. Bronze 2 (bz2) in maize was the first tau class GST reported to be involved in anthocyanin accumulation; bz2 produces yellow skin kernels because of disabled anthocyanin transport into the vacuole [53]. Anthocyanin deposition is also affected by genes in the phi class, such as AtGSTF12 in Arabidopsis [54,55], FvRAP in strawberry [26], CsGSTF1 in tea [52], and VvGST4 in grapevine [56]. Interestingly, these phi-class GSTs show extensive functional diversification. For example, AtGSTF12 plays a key role in both anthocyanin and PA accumulation in Arabidopsis [54], while CsGSTF1 in tea functions only in anthocyanin accumulation [52]. Functional divergence of GSTs has arisen through nonsynonymous mutations, especially at key amino acid sites [51]. For example, a single amino acid mutation (Arg39 to Trp39) was found to be responsible for the high enzymatic activity of Populus euphratica PeGSTU30 [51]. Li et al. [57] showed that one Trp to Leu substitution at amino acid 205 in AtTT19 resulted in an anthocyanin-deficient phenotype in Arabidopsis. Recently, Luo et al. reported that a single nucleotide polymorphism (C to T) in the second exon of FvRAP dramatically decreased the anthocyanin level in the petiole and fruit of strawberry [26]. In our study, a phi-class GST gene, PcGSTF12, was among the core genes for anthocyanin and PA accumulation in pear. PcGSTF12 was strongly associated with most of the studied metabolites in pear: seven anthocyanins and seven PAs. Similar to AtGSTF12, PcGSTF12 was functionally characterized as an anthocyanin and PA carrier. Interestingly, we detected different affinities for anthocyanin and PA between PcGSTF12 and AtGSTs in Arabidopsis. PcGSTF12 plays an opposite role to that of AtGSTs in procyanidin A3 and petunidin 3, 5-diglucoside accumulation. These results suggest that phi-class GSTs have undergone extensive functional diversification during evolution. Furthermore, transformation with PcGSTF12 affected genes encoding proteins involved in anthocyanin and PA biosynthesis, regulation, and transport. This functional analysis of PcGSTF12 deepens our understanding of the roles of phi-class GSTs in anthocyanin and PA accumulation in pear and other plants.
We found that PcMYB114 can positively regulate PcGSTF12 activity by directly binding to the MBS motif in its promoter. A recent previous study has shown that overexpression of PyMYB114 in young pear fruit significantly enhance anthocyanin accumulation by upregulating anthocyanin structural genes PyDFR, PyANS, and PyUFGT [14]. However, whether MYB114 functions in PA accumulation or affects anthocyanin and PA transportation are largely unknown. This result provides new evidence that PcMYB114 functions both in anthocyanin and PA transportation in pear via regulating PcGSTF12 activity, which may further affect anthocyanin and PA accumulation. Furthermore, we detected many cis-acting elements involved in auxin, ethylene, and gibberellin signaling in the promoter of PcGSTF12. Thus, we propose that PcGSTF12 might be the common downstream target of auxin-, ethylene-, and gibberellin-mediated anthocyanin and PA accumulation pathways.

Conclusions
In this study, we identified 4065 DEGs and 19 differentially expressed metabolites (12 anthocyanins and seven PAs) between young fruits of the green pear 'CF' and its red mutation 'RCF'. Based on correlation analyses between DEGs and anthocyanins/PAs, we found 203 candidate genes for the accumulation of 10 anthocyanins and seven PAs. We further identified a 'core set' of 20 candidate genes for pear anthocyanin and PA accumulation. Of these, PcGSTF12 was functionally characterized as an important anthocyanin and PA carrier in pear. We also identified an important pear anthocyanin and PA regulation node consisting of two core genes, PcGSTF12 and PcMYB114. These results provide novel insights into pear anthocyanin and PA accumulation. The candidate genes for pear anthocyanin and PA accumulation presented here represent a valuable data set to guide future functional studies.

Plant materials
The cultivar 'RCF' is a typical red pear sport of cultivar 'CF' that was discovered in the USA. The fruit of 'RCF' is initially green, then changes quickly to red within 1 week after full bloom, and remains red until the fruit ripens. The coloration pattern of 'RCF' differs from that of most pear species, which color at the ripening stage. Thus, 'CF' and 'RCF' are ideal materials to study the molecular mechanism of anthocyanin and PA accumulation in young pear. The pear cultivars 'CF' and 'RCF' were cultivated in the experimental orchard of Yantai Academy of Agricultural Science, Shandong province, Yantai, China (37°5′N, 122°1′W). The fruits of 'CF' and 'RCF' were collected in 2017 from 6-year-old trees grafted onto Pyrus betulaefolia rootstocks. The different fruit skin color phenotypes of 'CF' and 'RCF' were visible at 5 DAFB. Thus, the regulation of anthocyanin and PA accumulation at this developmental stage of fruits is important for fruit coloration. We collected fruits of 'CF' and 'RCF' at two early developmental stages (2 DAFB and 5 DAFB) for further analyses. Briefly, fruits of 'CF' and 'RCF' with a similar green color were first sampled on 19 April 2017; these samples were collected at 2 DAFB, and were named 'CF1' and 'RCF1', respectively. Fruits with significant differences in color between 'CF' and 'RCF' were sampled on 22 April 2017; these samples collected at 5 DAFB were named 'CF2' and 'RCF2', respectively. In each experiment, skins from 100 fruits per replicate were collected. Three independent biological replicates were collected for analyses. The fruit skin samples were immediately frozen in liquid nitrogen and stored at − 80°C until further metabolite, RNAsequencing (RNA-Seq), and qPCR analyses.

Metabolite extraction and separation
Metabolite extraction and separation were carried out as described by Wang et al. [23]. Briefly, the freeze-dried fruit skin was crushed into powder and then extracted overnight at 4°C in 1.0 ml 70% aqueous methanol. Following centrifugation at 10, 000 g for 10 min, the extracts were filtered and analyzed by HPLC.

Anthocyanin and PA identification and quantification
Anthocyanin and PA metabolites were annotated by comparisons against public databases including KNAPSAcK, MassBank, MoToDB, METLIN and HMDB, and were quantified using MRM as described by Wang et al. [23].

Total RNA isolation and RNA-Seq analysis
Total RNA was isolated using Trizol reagent (Invitrogen, Carlsbad, CA, USA) and its integrity was evaluated using a 2100 Bioanalyzer (Agilent Technologies, Santa Clara, CA, USA). The mRNA was purified from high-quality total RNAs using oligo (dT) magnetic beads and then broken into short fragments with fragmentation buffer. The cDNA was synthesized using a cDNA Synthesis Kit (TaKaRa, Dalian, China) and linked to sequencing adapters at both ends. The cDNA libraries were sequenced using the Illumina sequencing system (HiSeq™ 2000, Illumina, San Diego, CA, USA). Clean reads were obtained using the NGS QC Toolkit [58]. Differential expression analysis was carried out using the DESeq R package (2012). The threshold for significant differential expression was P < 0.05, and |log 2 fold change| ≥ 1 was used to identify the differentially expressed genes (DEGs) between two different cDNA libraries. Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analyses of DEGs were performed using the R platform as described elsewhere [59].

Integrated metabolome and transcriptome analyses
Pearson's correlation coefficients were calculated between the metabolome and transcriptome data. The coefficients were calculated from log2 (fold change) of each metabolite and log2 (fold change) of each transcript by the EXCEL program. Correlations with a coefficient of R 2 > 0.8 were selected. Metabolome and transcriptome relationships were visualized using Cytoscape (version 2.8.2).

qRT-PCR validation
The total RNA samples used for RNA-Seq were also used for cDNA synthesis using the PrimeScript™ RT Reagent Kit (TaKaRa) according to the manufacturer's instructions. The qRT-PCR analyses were conducted as described previously [13] using the primers shown in Supplementary Table 1. PcActin was used as the reference gene. Three biological replicates were analyzed.

Phylogenetic analysis
Phylogenetic tree analysis of PcGSTF12 and its homologs was carried out by MEGA 7 with bootstrap values calculated from 1000 replicate analyses.

PcGSTF12 promoter analysis
The 2000-bp upstream sequence of PcGSTF12 was analyzed using tools at the PlantCARE (http://bioinformat ics.psb.ugent.be/webtools/plantcare/html/) and PLACE Signal Scan Search databases (https://www.dna.affrc.go. jp/PLACE/). The functional motifs are listed in Supplementary Table 10. A MYB-binding site (MBS) sequence CAACTG was found at − 801 bp in the promoter of PcGSTF12 by PlantCARE. This MBS element is known to be a binding sequence of R2R3-MYBs [27,60]. Therefore, we selected this MBS element sequence to further analyze the interaction between PcMYB114 and the promoter of PcGSTF12 using electrophoretic mobility shift assays (EMSA) as described below.

Electrophoretic mobility shift assay
The EMSA was performed using the LightShift Chemiluminescent EMSA Kit (Thermo Scientific, Waltham, MA, USA) as described by Jiang et al. [61]. The recombinant protein was purified using the Niagarose His-Tagged Protein Purification Kit (CWbiotech, Beijing, China).

Luciferase reporter assay
Luciferase reporter assay was carried out as described by Wang et al. [62]. The CDS of PcMYB114 was recombined into the pHBT-AvrRpm1 effector. The promoter of PcGSTF12 was inserted into the pFRK1-LUC-nos reporter. The activities of LUC and GUS were detected using a Multimode Plate Reader (Victor X4, PerkinElmer, http://www.perkinelmer.com/).