Transcriptomic Analysis of Leaf in Tree Peony Reveals Differentially Expressed Pigments Genes

Tree peony (Paeonia suffruticosa Andrews) is an important traditional flower in China. Besides its beautiful flower, the leaf of tree peony has also good ornamental value owing to its leaf color change in spring. So far, the molecular mechanism of leaf color change in tree peony is unclear. In this study, the pigment level and transcriptome of three different color stages of tree peony leaf were analyzed. The purplish red leaf was rich in anthocyanin, while yellowish green leaf was rich in chlorophyll and carotenoid. Transcriptome analysis revealed that 4302 differentially expressed genes (DEGs) were upregulated, and 4225 were downregulated in the purplish red leaf vs. yellowish green leaf. Among these DEGs, eight genes were predicted to participate in anthocyanin biosynthesis, eight genes were predicted involved in porphyrin and chlorophyll metabolism, and 10 genes were predicted to participate in carotenoid metabolism. In addition, 27 MYBs, 20 bHLHs, 36 WD40 genes were also identified from DEGs. Anthocyanidin synthase (ANS) is the key gene that controls the anthocyanin level in tree peony leaf. Protochlorophyllide oxido-reductase (POR) is the key gene which regulated the chlorophyll content in tree peony leaf.


Introduction
Tree peony (Paeonia suffruticosa Andr.) is a woody shrub of the genus Paeonia and family Paeoniaceae, and is a traditional ornamental flower because of its large, showy, and colorful flowers [1]. Besides its beautiful flower, the leaf of tree peony has also good ornamental value owing to its leaf color change in spring.
Coloration of plant organs such as flowers, fruit and leaves is governed by pigments in plants [2]. There are four main classes of pigments for coloration in plants: flavonoids/anthocyanins, chlorophylls, carotenoids, and betalains [3]. Flavonoids, particularly anthocyanins, acting as major pigments that responsible for the orange to blue colors found in many plant flowers, leaves, fruits, seeds, and other tissues [2][3][4]. Chlorophylls including chlorophyll a and chlorophyll b, the photosynthetic pigments that play essential roles in plant photosynthesis, are main pigments in plant green leaf [5]. Carotenoid belong to the isoprenoid pigments, which are usually responsible for yellow-orange colors [6,7]. Carotenoids are also photosynthetic pigments that usually exist in plant leaf with chlorophylls. Betalains are the yellow and violet pigments belong to the order Caryophyllales, and they only exist in very few plants such as Beta vulgaris and Suaeda salsa [8,9]. In plants, anthocyanins and betalains are usually mutual exclusive, which are never found together in the same plant [10].
Prior to advancements in this field, researchers mainly focused on floral flavonoid/anthocyanin pigment components detection [11][12][13], or flavonoid/anthocyanin genes of petal analysis in tree peony [14][15][16][17]. As far as we know, there has been no report on leaf pigments components detection and related gene analysis in tree peony. In the absence of a complete genome sequence in tree peony, transcriptomic analysis is an effective method for gaining insights into differentially expressed genes (DEGs) between different color leaves. Transcriptome sequencing of the different leaf color stages of tree peony will provide useful insights into the molecular mechanisms of tree peony leaf color changes in spring. Therefore, the leaf pigment levels of different stages in tree peony was analyzed. The different leaf color stage transcriptomes were compared, and the DEGs were filtered. Our main objective was to annotate and analyze the DEGs involved in flavonoid/anthocyanin biosynthesis, chlorophyll and carotenoid biosynthesis to identify the candidate genes responsible for tree peony purplish red leaf in spring.

Pigments Level in the Leaf of Three Developmental Stages in Tree Peony
In spring, as the tree peony leaf developed two major physiological changes, color fading and yellowish greening were visually apparent ( Figure 1a). The tree peony leaf color was purplish red at S1, then changed to half purplish red and half yellowish green at S2, finally yellowish green at S3. As shown in Figure 1b, the anthocyanin contents were decreased from S1 to S3. While the total flavonoid, carotenoid, and chlorophyll contents were increased from S1 to S3. These indicated that the change in leaf color phenotypes of tree peony 'Man Yuan Chun Guang' was mainly caused by the decrease of anthocyanin content and increase of chlorophyll and carotenoid. peony [14][15][16][17]. As far as we know, there has been no report on leaf pigments components detection and related gene analysis in tree peony. In the absence of a complete genome sequence in tree peony, transcriptomic analysis is an effective method for gaining insights into differentially expressed genes (DEGs) between different color leaves. Transcriptome sequencing of the different leaf color stages of tree peony will provide useful insights into the molecular mechanisms of tree peony leaf color changes in spring. Therefore, the leaf pigment levels of different stages in tree peony was analyzed. The different leaf color stage transcriptomes were compared, and the DEGs were filtered. Our main objective was to annotate and analyze the DEGs involved in flavonoid/anthocyanin biosynthesis, chlorophyll and carotenoid biosynthesis to identify the candidate genes responsible for tree peony purplish red leaf in spring.

Pigments Level in the Leaf of Three Developmental Stages in Tree Peony
In spring, as the tree peony leaf developed two major physiological changes, color fading and yellowish greening were visually apparent ( Figure 1a). The tree peony leaf color was purplish red at S1, then changed to half purplish red and half yellowish green at S2, finally yellowish green at S3. As shown in Figure 1b, the anthocyanin contents were decreased from S1 to S3. While the total flavonoid, carotenoid, and chlorophyll contents were increased from S1 to S3. These indicated that the change in leaf color phenotypes of tree peony 'Man Yuan Chun Guang' was mainly caused by the decrease of anthocyanin content and increase of chlorophyll and carotenoid.

Library Construction and Transcriptome Sequencing
To understand the molecular basis of leaf color change in tree peony, the leaves of S1, S2, and S3 were used to build three libraries for high-throughput sequencing. 64,530,684; 67,008,409; and 69,283,538 raw reads were obtained from S1, S2, and S3 sequencing libraries, respectively. After removal of adaptor sequences, ambiguous reads, and low-quality reads, 62,835,627; 65,363,790; and 67,681,919 high quality clean reads comprising 9.43, 9.80,10.15 Gb nucleotides with Q20 > 98% were generated from S1, S2, and S3 transcriptome sequencing. In total, 53,323 unigenes with an N50 of 1215 were generated byTrinity 2.1.1 (Broad Institute, Cambridge, MA, USA). The min length, max length and average length of these unigenes were 201, 13121, and 801 nt, respectively. The size distributions of these unigenes are shown in Supplementary Materials Figure S1.

Identification and Analysis of DEGs between Three Leaf Color Stages
The thresholds false discovery rate (FDR) < 0.05 and absolute value of log2 ratio > 1 were used to judge the significance of the DEGs at three leaf color stages (S1 vs. S2, S2 vs. S3, S1 vs. S3). There were 797 DEGs between S1 and S2 libraries, in which 575 unigenes were significantly upregulated and 222 unigenes were downregulated. A total of 7587 DEGs were found between S2 and S3 libraries, with 3627 downregulated and 3960 upregulated. Finally, the greatest number of differentially expressed genes (8527 DEGs) occurred between S1 and S3 libraries, with 4225 downregulated and 4302 upregulated ( Figure 3). Therefore, DEGs between S1 and S3 libraries were used for further analysis.
Unigene0036137 and Unigene0052871 were annotated as UDP-glycosyltransferase and glucoside glucosyltransferase, these two genes were related with anthocyanidin glycosylation. In anthocyanin biosynthesis pathway, the spatial and temporal expression of structural genes such as F3H, DFR, and ANS is usually controlled by transcription factors from MYB, bHLH, and WD40 families [18][19][20][21]. In our study 27 MYBs, 20 bHLHs, and 36 WD40 were identified from DEGs ( Table 2). Among these transcription factors, 15 MYBs, 15 bHLHs, and 29 WD40 were upregulated in purplish red leaf. It indicated that these upregulation transcription factors were positively related with anthocyanin level in tree peony leaf.

Real-Time Quantitative PCR Analysis of the DEGs Involved in Pigment Biosynthesis
To test the reliability of the RNA-Seq data, 20 DEGs participated in tree peony leaf pigments biosynthesis including 6 DEGs involved in anthocyanin biosynthesis, 4 DEGs involved in chlorophyll biosynthesis, and 10 DEGs involved in carotenoid biosynthesis were selected for real-time quantitative PCR (q-PCR) analysis. Unigene0020903 (ubiquitin) was used as internal control. The expression patterns revealed by q-PCR analysis were similar to those revealed by RNA-Seq for the same genes ( Figure 4).

Discussion
Non-green leaves occur in many plants which are very popular in modern landscape garden with its unique ornamental values. Prior to the previous reports, there are three main pigments responsible for non-green leaves. Anthocyanidin accumulation was a major reason for the formation of red color leaves [25][26][27]. Carotenoid was usually responsible for the formation of yellow color [6,28]. Chlorophyll is mainly responsible for the formation of green color [28,29]. High concentrations of chlorophylls can mask the red or yellow color that is provided by anthocyanin or carotenoid [29,30]. In the present study, leaf color altered from purplish red to yellowish green complied with leaf development of tree peony. The anthocyanidins level of the purplish red leaf was significantly higher than that of yellowish green leaf, while chlorophylls and carotenoid level of purplish red leaf was significantly lower than that of the yellowish green leaf. These results indicated that purplish red leaf in tree peony could be attributed to anthocyanidin accumulation. The decrease of anthocyanidin, the increase of chlorophyll and carotenoid were mainly responsible for the alteration from purplish red to yellowish green of tree peony leaf in spring.
According to RNA-Seq data, four anthocyanin structural genes including F3H, F3'H, DFR and ANS were upregulated in the purplish red leaf and downregulated in the yellowish green leaf. While another two anthocyanin structural genes (CHI and HCT) were downregulated in the purplish red leaf and upregulated in the yellowish green leaf. Which were similar to those revealed by q-PCR analysis for the same genes. CHI is the enzyme that catalyze 2 ,4 ,4 ,6 -tetrahydroxy chalcone form to naringenin (a flavanone). HCT can catalyze p-Coumaroyl-CoA produce Caffeoyl-CoA, then synthesize eriodictyol (a flavanone). F3H belongs to the OGD family, which catalyzes naringenin yield dihydrokaempferol. Dihydrokaempferol is catalyzed by F3'H to produce dihydroquercetin [3]. Dihydrokaempferol and dihydroquercetin are catalyzed by DFR to form leucoanthocyanidins, and then generate colored anthocyanidins catalyzed by ANS [3]. High expression of CHI and HCT can produced more flavones, while anthocyanin production can be restrained by low expression of F3H, F3 H, DFR, and ANS. Owing to higher expression of CHI, HCT, and lower expression of F3H, F3 H, DFR, ANS in yellowish green leaf caused the higher accumulation of non-anthocyanin flavonoids (flavones or flavonols), which was verified by leaf pigment level detection (Figure 2). The pigment level analysis showed that anthocyanin content in leaf was consistent with the trend of leaf color variance. Therefore, co-expression of F3H, F3 H, DFR, and ANS was responsible for red pigmentation in the leaf of Paeonia suffruticosa 'Man Yuan Chun Guang'.
CHI is an upstream gene in anthocyanin biosynthesis that controls conversion of chalcones to (2S)-naringenin [3]. In this study, we found that CHI showed very higher expression level in yellowish green leaf stage compared with purplish red leaf stage. In P. suffruticosa 'Luoyang Hong', the expression level of F3H, F3'H, DFR, and ANS was higher in red petal, while the expression level of CHI in red petal was the lowest [14]. The similar expression pattern of CHI was also found in P. lactiflora [31]. While in the leaf of Ginkgo biloba, an organ lacking anthocyanin, the highly expression level of CHI was detected [32]. Li et al. (2006) reported that transgenic tobacco plants over-expressing sense Asteraceae CHI appeared normal in flower color along with an elevated accumulation of flavonoids, but not anthocyanin [33]. In addition, over-expression of CHI from tree peony can increased the level of total flavones and flavonols but reduced the intensity of flower pigmentation in transgenic tobacco [16]. The expression of CHI not only involved in anthocyanin biosynthesis, but also other secondary metabolites biosynthesis, such as flavones and flavonols [3]. All of these indicted that CHI is the key gene affecting the accumulation of flavonoids in plant, but not anthocyanin.
ANS is a downstream gene in anthocyanin biosynthesis, which encoded anthocyanidin synthase that catalyzes the oxidation of the colorless leucoanthocyanidin to the corresponding colored anthocyanidin [3]. In the present study, the expression level of ANS was 6.43-fold higher in purplish red leaf vs. yellowish green leaf, which was significantly higher than other three candidate structural genes of anthocyanin biosynthesis (Figure 4a), suggesting that ANS could be the key gene of anthocyanin accumulation in tree peony leaf. In herbaceous peony, ANS had the highest levels in red flower cultivar and the lowest levels in white flower cultivar, which resulted in the shift from white to pink and red in flowers [4]. In five cultivars of sweet potato, ANS expression was strongly associated with anthocyanin accumulation [34]. The similar result was also found in Chrysanthemum grandiflorum [35] and Malus sp. [36]. In addition, the silencing of ANS in red colored apple cultivar caused the shift of red leaf to green leaf, suggesting that anthocyanin biosynthesis was almost completely blocked [37]. All these implied that ANS was the key gene affecting the accumulation of anthocyanin in plant.
In our present study, eight DEGs (Unigene0045551, Unigene0031269, Unigene0022913, Unigene0032398, Unigene0036679, Unigene0039869, Unigene0026448, and Unigene0047139) potentially involved in porphyrin and chlorophyll biosynthesis were found and analyzed (Table 3). Among these genes, Unigene0031269, Unigene0022913, Unigene0032398, Unigene0036679, Unigene0026448, and Unigene0047139 showed higher transcript abundance in green leaf which has higher chlorophyll level compared with purplish red leaf (Table 3). It seemed that these six genes play positive role in chlorophyll accumulation in tree peony leaf. However, Unigene0026448 was annotated as HEMH, which encoded protoheme ferro-lyase that catalyzed protoporphyrin IX into protoheme, which will cause the decrease of protoporphyrin IX, an important intermediate product of chlorophyll biosynthesis [38]. Unigene0047139 was annotated as PAO, which will cause the opening of the chlorine macrocycle present in chlorophyll, and finally chlorophyll is broken down to colorless linear tetrapyrroles [29,39]. Therefore, Unigene0026448 and Unigene0047139 were not likely the key genes that affect chlorophyll level in tree peony leaf. Unigene0031269 (CRD1) encoded Mg-protoporphyrinogen IX monomethylester cyclase that catalyze Mg-protoporphyrin IX monomethylester formed into divinyl protochlorophyllide. Unigene0022913 (DVR) encoded divinyl reductase which catalyze divinyl protochlorophyllide produced protochlorophyllide. Unigene0032398 encoded protochlorophyllide oxidoreductase which can catalyze protochlorophyllide into chlorophyllide a or catalyze divinyl protochlorophyllide into divinyl chlorophyllide a depending on the availability of substrates [38]. Unigene0036679 (CAO) encoded chlorophyllide a oxygenase which catalyze chlorophyllide a formed to chlorophyllide b. q-PCR data shows that Unigene0031269, Unigene0022913, and Unigene0032398 were expressed higher level in the green leaf compared with the purplish red leaf, especially Unigene0032398 (POR) was expressed very highly level (44.82-fold higher) in the green leaf (Figure 4b). These indicated that POR was the key gene affecting chlorophyll level in tree peony leaf.
In the present study, 10 DEGs were predicted to participate in carotenoid biosynthesis metabolism in the tree peony leaf (Table 4). According to RNA-Seq data (Table 4) and q-PCR (Figure 4c), these 10 DEGs were all expressed higher in S3 stage which has higher carotenoid content than in S1. It seemed that these genes all play positive role in carotenoid accumulation in tree peony leaf. However, Unigene0051809 was annotated as carotenoid cleavage dioxygenase gene that can coded 9-cisepoxycarotenoid dioxygenase (NCED), which catalyze 9-cis-violaxanthin and 9-cis-neoxanthin cleaved to xanthoxin, the precursor of abscisic acid, (ABA) [3], thus caused the decrease of carotenoid accumulation. Therefore, Unigene0051809 (NCED) was not likely the key gene that affect carotenoid level in tree peony leaf.

Plant Material
Paeonia suffruticosa Andr. 'Man Yuan Chun Guang' plants were grown under field conditions in Northwest A&F University, Yangling Shaanxi, China. Leaf samples were collected at three different leaf color stages (S1, S2, S3) in the morning during the March and April in 2016 (Figure 1a). In order to avoid the effect of other organisms (arthropods, bacteria, fungi, etc.) on gene expression, only healthy, undiseased and non pest injured leaves were collected. All samples were collected from the same plant and washed three times with distilled water, and immediately frozen in liquid nitrogen, then stored at −80 • C for RNA extraction, and flavonoid/anthocyanin, chlorophyll, and carotenoid analysis.

Anthocyanin Level Measurement
Freeze-dried leaves were finely ground and 0.5 g was extracted with 15 mL acidic methanol (0.1% hydrochloric acid) at 4 • C in darkness for 24 h, then suspended by ultrasound for 1 h. After centrifugation at 5000 rpm for 1 min, the supernatant was filtered using a 0.22 µm membrane filter. The absorbance was read at 530 nm using a UV-Vis spectrophotometer (UV-3802, Unico, Suite E Dayton, OH, USA). A standard calibration curve was plotted using cyanidin-3-O-glucoside (Cy3G). The values were normalized to the dry weight of each sample to indicate anthocyanin level.

Flavonoid Level Measurement
Freeze-dried leaves were finely ground and 0.5 g was extracted with 15 mL methanol at 4 • C in darkness for 24 h, then suspended by ultrasound for 1 h. After centrifugation at 5000 rpm for 1 min, the supernatant was filtered using a 0.22 µm membrane filter. 2 mL supernatant, 2 mL 1.5% AlCl 3 , 3 mL 1 M Acetic acid sodium acetate (pH 5.0) were mixed and keep the volume to 10 mL. After incubation for 10 min, the absorbance was read at 415 nm using a UV-Vis spectrophotometer. A standard calibration curve was plotted using rutin. The values were normalized to the dry weight of each sample to indicate total flavonoid level.

Carotenoid Level Measurement
Freeze-dried leaves were finely ground and 0.5 g was extracted with 15 mL petroleum ether at 4 • C in darkness for 24 h, then suspended by ultrasound for 1 h. After centrifugation at 5000 rpm for 1 min, the supernatant was filtered using a 0.22 µm membrane filter. The absorbance was read at 440 nm using a UV-Vis spectrophotometer. A standard calibration curve was plotted using β-caroteneas. The values were normalized to the dry weight of each sample to indicate carotenoid level.

Chlorophyll Level Measurement
Freeze-dried leaves were finely ground and 0.5 g was extracted with 15 mL 80 % acetone at 4 • C in darkness for 24 h, then suspended by ultrasound for 1 h. After centrifugation at 5000 rpm for 1 min, the supernatant was filtered using a 0.22 µm membrane filter. Chlorophyll level was determined according to the protocol of Arnon.
Three biological replicates were performed for the above pigments measurement.

RNA Extraction, Library Construction, and Transcriptome Sequencing
Total RNA was extracted using the modified CTAB method. RNase-free DNase I (Tiangen; Beijing, China) was applied to remove contaminating genomic DNA. The RNA purity was determined with a NanoDrop 2000 spectrophotometer (NanoDrop Technologies; Wilmington, DE, USA), and 2% agarose gels were run to verify RNA integrity. The construction of the libraries and the RNA-Seq were performed by the Gene Denovo Corporation (Guangzhou, China). mRNA was enriched with oligo (dT)-rich magnetic beads. Then the enriched mRNA was fragmented into short fragments using fragmentation buffer and reverse transcripted into cDNA with random primers. Second-strand cDNA were synthesized by DNA polymerase I, RNase H, dNTP, and buffer. Then the cDNA fragments were purified with QiaQuick PCR extraction kit, end repaired, poly(A) added, and ligated to Illumina sequencing adapters. The ligation products were size selected by agarose gel electrophoresis, PCR amplified, and sequenced using Illumina HiSeq™ 4000 by Gene Denovo Biotechnology Co. Two biological replicates were used in this section.

De Novo Assembly and Gene Annotation
The raw reads (which has be uploaded in the NCBI's Sequence Read Archive, SRP095498) were first filtered by removing adaptor sequences, reads containing more than 10% of unknown nucleotides (N) and low quality reads containing more than 50% of low quality (Q-value ≤ 10) bases using Perl script. The remaining high-quality reads were assembled de novo using Trinity 2.1.1. [40]. In assembly, the reads from all three libraries were combined, and then ran Trinity. The longest transcript (from alternative splicing transcripts) was selected as the unigene in this study.

Identification and Analysis of DEGs
The unigene expression was calculated and normalized to RPKM (Reads Per kb per Million reads) [41]. To identify differentially expressed genes across samples, the edgeR package (http://www.r-project.org/) was used. We identified genes with a fold change ≥ 2 and a false discovery rate (FDR) < 0.05 in a comparison as significant DEGs. The confirmed DEGs were subjected to GO functional enrichment analysis and KEGG pathway analysis. Then, based on KEGG pathway analysis, the DEGs involved in leaf coloration were screened for up/downregulated unigenes, and used for further analysis.

Real-Time Quantitative PCR Analysis
Total RNA was extracted from leaves in S1 and S3 stages, respectively. After treatment with RNase-free DNase I (Tiangen; Beijing, China) according to the user manual, 1 µg of total RNA was reverse-transcribed to first-strand cDNA using the PrimeScript RT reagent kit (Takara; Otsu, Japan). Real-time quantitative PCR (q-PCR) experiments was performed using SYBR Premix Ex Taq (Takara) on ABI 7500 system (Applied Biosystems, Waltham, MA, USA). Primers used are listed in Table S2. Initial denaturation at 95 • C for 10 min, followed by 40 cycles of denaturation for 15 s at 95 • C, annealing for 30 s at 60 • C, and extension for 15 s at 72 • C, followed by a final extension for 15 s at 95 • C. The ubiquitin (Unigene0020903) was used to normalized the q-PCR data [42]. The relative expression levels of genes were calculated using the 2 −∆∆ct method. The statistical p-value was generated by the paired t-test. The statistical significance was defined as p < 0.05. Three biological replicates were performed for each gene.

Conclusions
Our results indicate that ANS was the key genes that control the anthocyanin level in tree peony leaf. POR was the key gene which regulated the chlorophyll content in tree peony leaf.