Extensive Transcriptome Changes Underlying the Flower Color Intensity Variation in Paeonia ostii

Tree peonies are a group of traditional ornamental plants, especially in East Asia, with Paeonia ostii as one of the most important ancestral species. P. ostii has flowers with varying colors, ranging from nearly white, light pink to deep pink. However, few studies have been done to unravel the molecular mechanisms underlying the flower color intensity variation in plants. Based on comparative analyses of the pigment composition and transcriptomes of P. ostii with different flower color intensities, we found that the anthocyanin concentration was significantly correlated with the flower color intensity in P. ostii. Transcriptome analysis by RNA-Sequencing revealed 7187 genes that were differentially expressed between flowers with different color intensities. Functional enrichment analysis of differentially expressed genes revealed multiple pathways possibly responsible for color intensity variation in P. ostii, including flavonoid biosynthesis, fatty acid oxidation, carbohydrate metabolism, and hormone-mediated signaling. Particularly, while anthocyanin biosynthesis genes showing positive correlations between their expression and anthocyanin concentration in flowers, two transcription factors, PoMYB2 and PoSPL1, seem to negatively regulate anthocyanin accumulation by affecting the activation capacity of the MYB-bHLH-WDR complex, exhibiting an inverse relationship between their expression and anthocyanin accumulation. Our results showed that, although anthocyanin biosynthesis had a direct effect on the pigmentation of the P. ostii flower, other metabolic and hormone-mediated signaling pathways were also contributed to the flower color intensity variation in P. ostii, suggesting complex coordinated changes in the transcriptional network. Differential expression of genes encoding anthocyanin repressors seems to be the major factor responsible for the intensity variation in anthocyanin pigmentation in P. ostii.


INTRODUCTION
Flower color is one of the most attractive sceneries in nature. It confers flowers with diverse functions (Winkel-Shirley, 2001;Steyn et al., 2002;Nagata et al., 2003), and is of paramount importance to plant evolution (Davies et al., 2012;Schiestl and Johnson, 2013;Sobel and Streisfeld, 2013). Angiosperms exhibit an astonishing polymorphism in flower colors. The phenotypic polymorphism in flower pigmentation is typically manifested in two types: (i) variation in pigment intensity determined by the concentration of pigment, and (ii) variation in floral hue, which is generally determined by the distinction of pigment types or the absence/presence of co-pigments (Wessinger and Rausher, 2012;Sobel and Streisfeld, 2013). In contrast to a wealth of knowledge on the genes related to qualitative variations between different pigment hues and between absence and presence of pigments (Zufall and Rausher, 2003;Hopkins and Rausher, 2011;Wessinger and Rausher, 2012;Sobel and Streisfeld, 2013), very few studies have been performed to characterize the genetic and molecular mechanisms determining the quantitative variation in flower color intensity (Schwinn et al., 2006;Ohno et al., 2013;Wang et al., 2014).
Paeonia ostii is a perennial shrub in the genus Paeonia and reproduces sexually (Li et al., 2011). It has been reported that P. ostii is one of the most important ancestral species of the cultivated tree peony, which is an important ornamental crop in the word and is crowned the "king of flowers" in China Zhou et al., 2014b). Flowers of P. ostii exhibit color polymorphism within populations, ranging from nearly white, light pink to deep pink. The pigmentation characteristics of P. ostii flowers make them an excellent model for studying the molecular basis of the intensity variation in pigmentation.
Generally, the variation in anthocyanin concentration is responsible for the flower color intensity (Grotewold, 2006;Tanaka et al., 2008). The quantitative change in anthocyanin intensity is accompanied by the amount alteration of flux through the pathway (Sobel and Streisfeld, 2013). It has been proposed that either increasing the functional activity of pathway enzymes and activation regulators or removing the repressors of anthocyanin production in flowers can result in an increase in anthocyanin intensity (Sobel and Streisfeld, 2013). It is unclear, however, whether intensity variation in flower pigmentation in P. ostii is generated by alterations of expression of anthocyanin biosynthesis genes or anthocyanin repressor genes. In this study, we compared the pigment composition and transcriptomes of P. ostii flowers with different intensity of coloration. We aimed to explore the correlations between color intensity and anthocyanin concentration, and to identify transcriptional changes and candidate genes potentially responsible for the control of pigmentation intensity in P. ostii. The results would provide insights into the molecular basis underlying the intensity variation in flower pigmentation in P. ostii.

Plant Materials
P. ostii was grown in the peony planting base of Fenghuangshan, Tongling, Anhui, China (Figure 1A). At full-bloom stage, flower color was analyzed following the International Commission on Illumination (CIE) system. The L * (lightness), a * (redness and greenness), and b * (yellowness and blueness) were measured using a hand-held spectrophotometer (NF333, Nippon Denshoku Industries Co., Ltd., Tokyo, Japan). For each flower, three areas of the medial surface were measured. For each individual plant, the measurement was performed with three petals from three independent flowers. The L * is an indicator of flower color intensity, as lower lightness generally means deeper color. Four classes of color intensity were chosen for this study: L * ranges of >85, 72-75, 65-68, and 57-60 corresponded to color class I, II, III, and IV, respectively. For each color intensity class, four petal samples were collected from different plants and immediately frozen in liquid nitrogen and stored at −80 • C until required for anthocyanin analysis and RNA extractions. Each sample contained three to four full-bloom flowers.

Pigments Extraction and High-Performance Liquid Chromatography (HPLC) Analyses
Petal samples from 16 plants belonged to four color intensity classes (four plants each) were extracted separately in acidic methanol solution (acetic acid: methanol: water = 1: 4: 5 v/v; 5 ml/100 mg tissue) for 24 h at 4 • C and passed through a syringe filter (0.22 µm; Pall). Extracts were analyzed on a Nova Pak C18 column (Waters) and detected using a photodiode array detector (Waters) in the range 190-600 nm. Peonidin-3,5-di-Oglucoside (Pn3G5G) and cyanidin-3,5-di-O-glucoside (Cy3G5G) were identified using the retention time by comparing with the HPLC data obtained from the commercially available standards (Polyphenols Laboratory, Sandnes, Norway). All determinations were performed with three replicates. Peak area was recorded.
Total anthocyanins and anthoxanthins (flavone and flavonol) were determined semi-quantitatively using the peak area normalization method by employing simple linear regressions using standards Pn3G5G and Cy3G5G for anthocyanins and Rutin for anthoxanthins at 525 and 350 nm, respectively. The content of anthocyanin and anthoxanthin was calculated in mg per 100 mg of fresh petals (as a sum of quantity of Pn3G5G and Cy3G5G mg/100 mg, and as a quantity of Rutin mg/100 mg, respectively).

RNA Extraction and Transcriptome Sequencing
Total RNA was extracted from the petals using TRIzol Reagent (Life Technologies Corp., Carlsbad, CA, USA) following the manufacturer's standard protocol, then treated with RNasefree DNaseI (Ambion, USA). RNA purity was assessed using a Nanodrop 2000C spectrophotometry (Thermo scientific, USA). Another precipitation step with 0.1 volume of 3 M sodium acetate and 2.5 volumes 100% (vol/vol) ethanol, was conducted if the 260/230 absorbance ratio of the total RNA was less than 1.5. RNA quality was assessed on a 2100 Bioanalyzer (Agilent Technologies, Santa Clara, CA, USA); RNA integrity number (RIN) value was greater than 7.5. , each value is shown as average ± standard deviation; different superscript letters on the horizontal axis labels indicate statistically significant differences between means of different color classes, as judged by t-test (P < 0.05). In plot (D), the R-and P-values indicate the correlation between L* and anthocyanin concentration.
We performed transcriptome sequencing on the bulked RNA of four plants for each flower color intensity class. cDNA library construction and Illumina sequencing were carried out at Beijing Genomics Institute (BGI)-Shenzhen, Shenzhen, China (http://www.genomics.cn/index.php) following the Illimina manufacturer's instructions (Illumina, San Diego, CA, USA). Briefly, poly-A RNA was enriched from 10 µg of total RNA using Magnetic beads with oligo (dT) and broken into short fragments with fragmentation buffer. Using these short fragments as templates, first-strand cDNA was synthesized using random hexamer primer. Then, second-strand cDNA was synthesized using buffer, dNTPs, RNase H (Invitrogen) and DNA polymerase I (Invitrogen). Following size selection and PCR amplification, the cDNA library was sequenced in a HiSeq 2000 to generate paired-end reads. After removing Illumina adapters and reads with unknown nucleotides larger than 5%, and trimming low-quality bases, the remaining high quality reads (clean reads) with an average length of 90 bp were used in this study. The raw sequence data sets were deposited in the US National Center for Biotechnology Information (NCBI) Sequence Read Archive (SRA, http://www.ncbi.nlm.nih.gov/Traces/sra) (Wheeler et al., 2008) under accession number SRP058369, including SRR2027814, SRR2027815, SRR2027817, SRR2027818, SRR2027819, SRR2027820, SRR2027821, and SRR2027822.

De novo Assembly and Functional Annotation
Transcriptome de novo assembly was performed with a short reads assembling program-Trinity (Grabherr et al., 2011). We further used a rapid clustering tool-TGICL (Pertea et al., 2003) to assemble unigenes from all four libraries to obtain a single set of non-redundant unigenes. Unigenes were annotated by BLASTX searches against the NCBI non-redundant protein (Nr) database (http://www.ncbi.nlm.nih.gov), Swiss-Prot protein database (http://www.expasy.ch/sprot) and Arabidopsis protein database at the Arabidopsis Information Resource (TAIR, http:// www.arabidopsis.org) with an E-value threshold of 10 −5 . Only the top hit for each sequence was extracted. According to the homology annotation against Arabidopsis protein database and NCBI Nr database, gene ontology (GO, http://www. geneontology.org) annotation of unigenes was obtained using Blast2GO program (Conesa et al., 2005). Unigene sequences were also aligned to the Kyoto Encyclopedia of Genes and Genomes Pathway (KEGG; http://www.genome.jp/kegg) database (Ogata et al., 1999) by BLASTx program using an E-value cutoff of 10 −5 to predict the metabolic pathway annotation.

Digital Gene Expression Profiling and Screening of Differentially Expressed Genes
To obtain digital gene expression profiles for different flower color intensity, all of the clean reads from four color classes were separately mapped to the non-redundant reference transcriptome sequences (all-unigenes). SOAPaligner/soap2 (Li et al., 2009) was used for mapping of the reads because of its high alignment accuracy and rapid alignment speed (Hatem et al., 2013;Shang et al., 2014). The number of unambiguously mapped clean reads for each gene in each sample was separately counted. The gene expression level was calculated using the Fragments Per Kilobase of exon per Million fragments mapped (FPKM) method (Mortazavi et al., 2008). Differentially expressed genes were identified based on the method described by Audic and Claverie (1997). The false discovery rate (FDR) was adopted to correct P-values in multiple hypothesis tests. A genes was judged to be differentially expressed if it had an FDR = 0.001 and the absolute value of log 2 Ratio ≥ 2. Unigenes with fewer than 50 reads in all samples were excluded from analysis.

Go Enrichment Analysis
GO categories enrichment analysis was implemented using the Bioconductor topGO package (Gentleman et al., 2004;Alexa et al., 2006), with the default arguments and Fisher's exact test to evaluate statistical significance.

Real-Time Quantitative Reverse Transcription Polymerase Chain Reaction (qRT-PCR) Validation and Expression Analysis
qRT-PCR analysis was performed to validate the expression pattern of selected genes identified by the digital expression analysis. The RNA separately extracted from all of the 16 petal samples was used. First strand cDNA was synthesized from 1 µg of total RNA using PrimeScript RT (Perfect Real Time) kit (TAKARA, Japan). The correctness of the gene sequences in the reference transcriptome was validated by reverse transcription PCR and TA cloning using PMD19-T vector kit (TAKARA), followed by sequencing. The qRT-PCR reactions with genespecific primers (Supplementary Table 1) were performed using a Real-time PCR Detection Systems (RocheCycler 480, Roche, Germany) and SYBRGreen PCR Master Mix (Roche, Germany) according to the manufacturer's instructions. Three independent biological replicates were performed for each reaction. Expression levels of the selected unigenes were normalized to that of an internal reference gene, glyceraldehyde 3-phosphate dehydrogenase (GAPDH), which was revealed as a stably expressed gene by our Illumina sequencing data and previous studies . As the PCR efficiency for all the gene-specific primers ranged between 93 and 107% over 1000-fold of cDNA dilution (Supplementary Figure 1 and  Supplementary Table 1), relative expression levels were calculated using the 2 − Ct method (Livak and Schmittgen, 2001).

Phylogenetic Analysis
Multiple sequence alignments were produced by ClustalW using default settings, and the phylogenetic trees were created using the neighbor-joining method and bootstrap analysis (1000 replicates) in MEGA5 software (Tamura et al., 2011). The tree in Supplementary Figure

Flower Pigmentation
The typical view of the growth location of P. ostii is shown in Figure 1A. Based on the lightness (L * ) of petals, the flowers of   (Figures 1B,C). Detection of anthocynins present in various flowers showed that the anthocyanin concentration was inversely correlated with L * (Figures 1B,D). A relatively high level of Peonidin-3,5-di-O-glucoside (Pn3G5G) and a trace amount of cyanidin-3,5-di-O-glucoside (Cy3G5G) were detected in medium and deep pink flowers. A low level of Pn3G5G was detected in the light pink flower, while Cy3G5G was almost undetectable. Both Pn3G5G and Cy3G5G were undetectable in the white flower (Supplementary Figure 2). Pn3G5G and Cy3G5G are synthesized based on peonidin and cyanidin, respectively. Cyanidin and peonidin are two similar anthocyanidins derived from the same branch of the anthocyanin biosynthesis pathway, differing in that peonidin has a methyl substitution on the B-ring.

Transcriptome De novo Assembly and Functional Annotation
In order to identify genes associated with flower color variation in P. ostii, we separately sequenced the petal transcriptomes of flowers of different color classes.  Table 2). 29,621, 29,433, 30,213, and 30,293 unigenes could be assigned to at least one Gene Ontology (GO) term in class I, II, III, and IV, respectively. Base on the GO term annotation, the unigenes were categorized into 48 functional groups (Figure 2). The distribution patterns of unigenes from the four classes were similar under different GO categories (Figure 2). Unigenes assigned to categories of "cellular process, " "catalytic activity, " and "metabolic process" were dominant. More than 470 unigenes were assigned to "flavonoid metabolic process" FIGURE 3 | Changes in gene expression profile between P. ostii petals of different color classes. The number of up-and down-regulated genes between class II and class I, class III and class I, and class IV and class I are summarized. For example, "892" means the expression of 892 genes was significantly higher in color class II than in color class I; and "1274" means the expression of 1274 genes was significantly lower in color class II than in color class I. and "pigmentation" categories in each class of flowers. KEGG pathway analyses showed that 18,130, 18,083, 18,637, and 18,770 unigenes in class I, II, III, and IV respectively, could be assigned to 128 KEGG pathways. The most represented pathways were "metabolism" (containing around 21% unigenes in each class), "biosynthesis of secondary metabolites" (around 11% unigenes in each class) and "plant hormone signal transduction" (around 5% unigenes in each class) (Supplementary Table 4).

Differentially Expressed Unigenes Between Flowers with Varied Color Intensities
Gene expression profiles of petals were compared between flowers with different color intensities. 55.50, 57.28, 55.27, and 55.47% reads of class I, II, III and IV libraries were aligned uniquely to the reference transcriptome obtained via combined assembly. A total of 7187 unigenes were differentially expressed between flowers with different color intensities (Supplementary Table 5). Most of the differentially expressed unigenes were found between pink and white flowers. The highest number of differentially expressed unigene was observed between flowers of class I and IV (Figure 3). One hundred and twenty-eight unigenes were found to be commonly up-regulated in various pink flowers, while 424 unigenes being commonly downregulated. To validate our expression data obtained by RNA sequencing, the differential expression patterns of 19 unigenes were verified by qRT-PCR (Supplementary Figure 3). A high correlation was observed between expression levels obtained by RNA-Seq and qRT-PCR (R = 0.787, P = 1.90E-13) (Supplementary Figure 3).
Genes involved in anthocyanin biosynthesis and transport were identified from the P. ostii transcriptome (Supplementary Figure 4 and Supplementary Table 6), and their expression levels in different classes of flowers were measured by RNA-Seq and qRT-PCR. The results showed that genes encoding chalcone isomerase (CHI; PoCHI3), flavanone 3-hydroxylase (F3H; PoF3H13, PoF3H18, and PoF3H20), flavonoid 3 ′ -hydroxylase (F3 ′ H; PoF3 ′ H11, and PoF3 ′ H12), dihydroflavonol reductase (DFR; PoDFR6), and leucoanthocyanidin dioxygenase (LDOX; PoLDOX1) were expressed at a higher level in various pink flowers than in white ones ( Figure 5). Notably, the expression levels of PoDFR6 and PoLDOX1 were well correlated with flower color intensity and anthocyanin concentration (Figure 5). These two genes were involved in the later steps of anthocyanin biosynthesis. The expression of genes encoding chalcone synthase (CHS; PoCHS1, and PoCHS2) and UDP flavonoid glucosyl transferase (UFGT; PoUFGT1, PoUFGT2, PoUFGT3, PoUFGT4, and PoUFGT5) showed no significant difference between pink and white flowers (Figure 5, Supplementary Table  6). The anthocyanins are stored in the central vacuole in plants (Kitamura, 2006). The activity of vacuolar flavonoid transporters, such as multidrug and toxin extrusion (MATE) and ATPbinding cassette (ABC) transporters, can affect anthocyanins accumulation in vacuole, which in turn affects the intensity of pigmentation in plant tissues (Goodman et al., 2004;Gomez et al., 2009;Zhao et al., 2011). The increasing tendency for anthocyanin accumulation in pink flowers of P. ostii was also accompanied by increased expression of genes encoding MATE transporters (PoMATE1, PoMATE2, and PoMATE3) and ABC transporters (PoMRP1 and PoMRP2) (Figure 6), suggesting their potential roles in promoting anthocyanins accumulation in P. ostii.
The MBW ternary protein complex composed of R2R3-MYB and basic helix-loop-helix (bHLH) transcription factors as well as WD-repeat (WDR) proteins has been documented as a primary regulator in anthocyanin biosynthesis (Heim et al., 2003;Baudry et al., 2004;Zimmermann et al., 2004;Koes et al., 2005;Ramsay and Glover, 2005;Lepiniec et al., 2006;Petroni and Tonelli, 2011). However, genes encoding the active MBW complex components were not significantly differentially expressed among flowers with different color intensities. Instead, a gene, PoMYB2, was identified to exhibit an inverse relationship between its expression and anthocyanin accumulation ( Figure 7A). Phylogenetic analysis showed that PoMYB2 was homologous to Petunia PhMYB27 and strawberry FaMYB1 (Supplementary Figure 5), which have been suggested to act as transcriptional repressors of anthocyanin accumulation (Aharoni et al., 2001;Albert et al., 2011). A well-conserved bHLH interaction motif [D/E]Lx 2 [R/K]x 3 Lx 6 Lx 3 R (Zimmermann et al., 2004) was found in the R3 repeat domain in the predicted protein sequence of PoMYB2 (Figure 7B). A putative repression FIGURE 4 | Overrepresented GO terms amongst differentially expressed genes between P. ostii petals of different color classes. Only GO terms are shown where P < 0.01 in at least one of the six gene sets (up-or down-regulated in any of these three comparisons: class II vs. class I, class III vs. class I, or class IV vs. class I).
Frontiers in Plant Science | www.frontiersin.org  Figure 1. For qRT-PCR, each value was normalized relative to PoGAPDH expression and is shown as average ± standard deviation (Continued) FIGURE 5 | Continued from three biological replicate sampling. Different superscript letters on the horizontal axis labels indicate statistically significant differences between means of different color classes, as judged by t-test (P < 0.05). The R-and P-values given in the plot indicate the correlation between anthocyanin concentration and gene expression. CHS, chalcone synthase; CHI, chalcone isomerase; F3H, flavanone 3-hydroxylase; F3 ′ H, flavonoid 3 ′ -hydroxylase; DFR, dihydroflavonol 4-reductase; LDOX, leucoanthocyanidin dioxygenase.  Figure 1. For qRT-PCR, each value was normalized relative to PoGAPDH expression and is shown as average ± standard deviation from three biological replicate sampling. Different superscript letters on the horizontal axis labels indicate statistically significant differences between means of different color classes, as judged by t-test (P < 0.05). The R-and P-values given in the plot indicate the correlation between anthocyanin concentration and gene expression. MATE, multidrug and toxin extrusion transporter; MRP, ATP binding cassette transporter.
motif LNLDL conforming to the LxLxL type of ERF-associated amphiphilic repression (EAR) motif was present in the C terminal region of the predicted PoMYB2 protein sequence (Kagale and Rozwadowski, 2011) (Figure 7B). In addition to PoMYB2, PoSPL1, a gene homologous to the miR156-targeted Arabidopsis SPL protein AtSPL13 Wu et al., 2009), also exhibited a negative relationship between its expression and anthocyanin concentration (Figure 8A, Supplementary Figure 6). A conserved miR156 target site in the coding sequence of PoSPL1 was also identified (Figure 8B, Supplementary Data Sheet 1).

DISCUSSION
Anthocyanin concentration was closely correlated with the intensity of petal color in P. ostii. Anthocyanins were synthesized via the flavonoid pathway. The corresponding genes have been well characterized and been divided into early biosynthetic genes and late biosynthetic genes in dicotyledonous plants (Martin and Gerats, 1993;Mol et al., 1998;Petroni and Tonelli, 2011). The early biosynthetic genes, including CHS, CHI, F3H, and F3 ′ H, located upstream of the anthocyanin biosynthetic pathway and led to the production of flavonols and other flavonoid compounds, while the late biosynthetic genes, including DFR and LDOX, were downstream genes in the anthocyanin biosynthetic pathway and were specifically for anthocyanin biosynthesis (Lepiniec et al., 2006;Petroni and Tonelli, 2011). In P. ostii, the late anthocyanin biosynthetic genes were differentially expressed among flowers with different color intensities, displaying a strong correlation between their expression levels and the concentration of anthocyanin in petals, whereas none of the early biosynthetic genes exhibited expression patterns significantly positively correlated with the concentration of anthocyanin (Figure 5). The increased 2008; Stracke et al., 2009). Substrate competition seems to exist between FLS and DFR, controlling the metabolic flux through the branches of the flavonoid biosynthetic pathway underlying flower color intensity variation in P. ostii. A significant positive correlation was also observed between the expression levels of anthocyanin transporter genes and the anthocyanin levels in petals of different flowers (Figure 6). On the contrary, the putative anthocyanin-biosynthesis repressor genes were expressed at a high level in the nearly white petals and gradually reduced in deep color flowers, exhibiting an expression pattern inverse to anthocyanin biosynthetic and transporter genes (Figures 7, 8). The coordinated expression of genes involved in anthocyanin biosynthesis, anthocyanin transport and suppression of anthocyanin synthesis suggested a sophisticated regulatory system underlying anthocyanin pigmentation in P. ostii. It has been revealed that the combinatorial action of R2R3-MYB and bHLH transcription factors, along with WDR proteins, played an important role in regulating the transcription of anthocyanin genes (Heim et al., 2003;Baudry et al., 2004;Zimmermann et al., 2004;Koes et al., 2005;Ramsay and Glover, 2005;Lepiniec et al., 2006;Petroni and Tonelli, 2011;Thévenin et al., 2012;Xu et al., 2014). The action of the MBW complex is primarily specified by the activity of R2R3-MYB factors in the complex. Diverged members of the R2R3-MYB gene family had distinct roles in determining the action of the complex either to promote or inhibit the transcription of anthocyanin biosynthesis genes (Matsui et al., 2008). In P. ostii, genes homologous to the subunits of the active MBW complex were not differently expressed in flowers with different color intensity. Instead, a new R2R3-MYB gene encoding PoMYB2 highly similar to PhMYB27 (Figure 7B, Supplementary Figure 5) which can reduce anthocyanin production either by repressing the expression of key factors required for MBW complex thus preventing the MBW complex formation or by incorporating to MBW activation complexes to convert them into repressive complexes , was highly expressed in the nearly white flower and gradually reduced expression in flowers with increasing color intensity ( Figure 7A). In line with the inverse expression pattern of PoMYB2 with that of the late anthocyanin biosynthetic genes, the PoMYB2 protein contains an amino acid motif [D/E]Lx 2 [R/K]x 3 Lx 6 Lx 3 R, which is required for interaction with bHLH partners in MBW complex (Zimmermann et al., 2004), and a C terminal EAR repressor motif (LNLDL), which mediate transcriptional repression in plants (Kagale and Rozwadowski, 2011) (Figure 7B). A similar pattern of expression was also observed for the gene PoSPL1, which was homologous to AtSPL13 ( Figure 8A, Supplementary  Figure 6). Although the functional evidence on AtSPL13 is lacking to date, previous studies have shown that AtSPL13 is one of the miR156-targeted SPL genes, which also include AtSPL3, AtSPL9, and AtSPL10 Wu et al., 2009). Previous functional analyses revealed that the miR156targeted SPLs shared some common regulatory functions (Yu et al., 2010). For example, both AtSPL9 and AtSPL13 are involved in the regulation of trichome production (Yu et al., 2010;Preston and Hileman, 2013). A well-conserved miR156 target site was also present in the coding sequence of PoSPL1 ( Figure 8B). It has been confirmed that AtSPL9 can repress anthocyanin production through its "squelching activity" causing the depletion of transcription factors essential for activation (Gou et al., 2011). That is, AtSPL9 can serve as a competitive inhibitor that acts by binding R2R3-MYB proteins required for the formation of MBW complex for regulation of anthocyanin production (Gou et al., 2011). The expression pattern of PoSPL1 exhibited a negative correlation with the anthocyanin concentration (Figure 8A), implying a possible role of PoSPL1 in repressing anthocyanin synthesis in P. ostii. All these findings suggested the existence of a complex regulatory network controlling anthocyanin synthesis in P. ostii. Negative regulators identified in P. ostii seem to be the major factor responsible for the intensity variation of P. ostii in anthocyanin pigmentation. It is likely that decreased expression of PoMYB2 and PoSPL1 released the repression of anthocyanin structural genes, that in turn increased the flower color intensity in P. ostii.
The decreased expression of PoMYB2 and PoSPL1 in deeper color flowers was unlikely caused by loss-of-function mutations because we did not find critical differences in the coding sequences of PoMYB2 and PoSPL1 from different intensity variants. Additionally, these two genes belong to different gene families and function through distinct mechanisms. The probability is small for PoMYB2 and PoSPL1 to mutate simultaneously to result in a coordinated decrease in expression in P. ostii. PoMYB2 and PoSPL1 might be coordinately regulated by other elements. Gene expression profiling showed that many genes involved in MAPKKK cascade, respiratory burst, phytohormone biosynthesis, and various organic substance metabolic processes were differentially expressed among different intensity variants (Figure 4). The expression levels of two genes (Unigene12751_All, Unigene10302_All) encoding sucrose synthases were positively correlated with the flower color intensity (Supplementary Table 5). Sucrose can promote the anthocyanin production by suppressing the expression of the transcription factor MYBL2 which negatively regulates anthocyanin biosynthesis, while concurrently inducing the expression of the positive regulators, including PAP1, TT8, and GL3 (Solfanelli et al., 2006;Jeong et al., 2010;Das et al., 2012). Previous studies in Arabidopsis showed that the accumulation of anthocyanin in leaves is positively correlated with the increase in endogenous sucrose content (Jeong et al., 2010). Genes involved in ethylene mediated signaling pathway were down-regulated in deeper color flowers (Figure 4), that is of particular interest because it has been revealed that ethylene plays a negative regulatory role in anthocyanin biosynthesis and functions by regulating the activity of positive and negative regulators at the transcriptional level (Jeong et al., 2010;Das et al., 2011Das et al., , 2012Qi et al., 2011). Ethylene also has a suppression effect on the sucrose-induced anthocyanin pigmentation (Jeong et al., 2010). Arabidopsis mutants defective in ethylene signaling or wild-type plants treated with ethylene biosynthesis and ethylene-binding inhibitors displayed a promotion on anthocyanin accumulation (Jeong et al., 2010). The alteration of the internal environment resulted from changes in sucrose and ethylene synthesis and signaling probably had an effect on transcriptional variation of anthocyanin repressor genes in different color variants. Fatty acid beta-oxidation was also found among the enriched GO categories in the upregulated gene sets in deeper color classes (Figure 4). The up-regulation patterns of genes encoding acyl-CoA oxidase (PoACX, Unigene12172_All) and peroxisomal 3-ketoacyl-CoA thiolase (PoPKT, Unigene29334_All, Unigene33907_All) were coordinated with the increase in flower color intensity (Supplementary Table 5). This result probably has implications for the re-allocation of acyl-CoA between fatty acid and anthocyanin biosynthesis. As the final product of fatty acid betaoxidation, acyl-CoA can be used as a substrate in anthocyanin biosynthesis. The induction of enzymes involved in fatty acid beta-oxidation may result in split-flowing of acyl-CoA and eventually enhance the production of anthocyanin in deeper color petals.
Tree peony is an important ornamental crop in the word, especially in East Asia. Classical breeding techniques have enabled the creation of a wide range of flower colors for tree peony cultivars. It has been look forward, however, to produce more varieties with desirable and novel flower colors by genetic engineering. P. ostii is one of the most important ancestral species of the cultivars of tree peony. The results of this study provided deep insights into the molecular basis underlying the variation in flower color intensity in tree peony, providing a foundation for developing novel flower colors through manipulating the anthocyanin biosynthesis pathway.

CONCLUSIONS
Flowers of P. ostii exhibited variations in color intensity, which was significantly correlated with the anthocyanin concentration. Genome-wide transcription analysis showed that, although anthocyanin biosynthesis had a direct effect on the pigmentation of P. ostii flowers, other metabolic and hormone-mediated signaling pathways were also contributed to the color intensity variation in P. ostii. Differential expression of genes encoding anthocyanin repressors, which negatively regulate the expression of anthocyanin biosynthesis genes by affecting the activation capacity of the MYB-bHLH-WDR complex, seems to be the major factor responsible for the intensity variation in anthocyanin pigmentation in P. ostii. This study provided insights into the potential key components related to the regulation of flower color intensity in P. ostii. As P. ostii is one of the most important ancestral species of the cultivars of tree peony, these results might provide a foundation for developing novel flower colors in ornamental cultivars through manipulating the anthocyanin biosynthesis pathway.

AUTHOR CONTRIBUTIONS
LG designed the research, performed the experiments and the data analysis, and drafted the manuscript. HY contributed analysis tools and participated in the data analysis and manuscript preparation. HL helped in the HPLC experiments and data analysis; JY participated in the design of the study, helped in data analysis and manuscript preparation. YH conceived the idea, and participated in the design of the study and in interpreting results, and manuscript preparation. All authors carefully read and approved the final manuscript.