Comparative Transcriptome Analysis of Anthocyanin Biosynthesis in Pansy (Viola × wittrockiana Gams.)

Abstract

Pansy (Viola × wittrockiana Gams.) is an important and attractive ornamental plant with a wide variety of flower colors. To date, the molecular genetic understanding of its colorful petal pigment patterns remains largely unknown. In this study, we analyzed the bicolor petals of “Mengdie” in cytological, physiological, and transcriptomic aspects. Results showed that the difference of flower colors was mainly determined by the pigment distribution in the epidermal cells. Pigment analysis indicated that anthocyanins had strong correlations with color parameters, which acted as the main factor in flower coloration. Comparative transcriptome analysis found a total of 43,908 unigenes with the mean length of 682 bp. There were 24,323, 16,668, 8507, and 7680 unigenes annotated in the Nr, Swiss-Prot, KOG, and KEGG databases, respectively. Differential expression genes (DEGs) showed that the expression of anthocyanin late biosynthesis genes (LBGs), VwF3′H, VwF3′5′H, and VwUFGT, would be likely to play a major role in the color formation of pansy. The expression patterns of selected DEGs were verified by qRT-PCR. This study contributes an excellent insight into molecular mechanism of pigment biosynthesis and provides some useful information for flower color modification in pansy.

1. Introduction

Flower color is an important aesthetic trait that affects commercial value of ornamental plants. Anthocyanins, as a main group of flavonoid, are responsible for the majority of orange, red, purple, and blue flowers. [1]. At present, the regulation of anthocyanin biosynthesis has been extensively studied in many plants [2,3]. A number of critical biosynthetic genes have been identified in the pathway, such as chalcone synthase (CHS), chalcone isomerase (CHI), flavanone 3-hydroxylase (F3H), flavonoid 3′-hydroxylase (F3′H), flavonoid 3′,5′-hydroxylase (F3′5′H), dihydroflavonol 4-reductase (DFR), and anthocyanidin synthase/leucoanthocyanidin dioxygenase (ANS/LDOX) [4,5].

Furthermore, the biosynthetic genes were temporally and spatially regulated by transcription factors (TFs) during flower color formation [6,7]. Three main types of TFs, including MYB, bHLH, and WDR, have been characterized to regulate anthocyanin biosynthesis individually or coordinately by forming MBW complex [1,7,8]. Among these TFs, MYB proteins have been considered to be the key components to control anthocyanin biosynthesis [9], which has been supported by these findings in many horticultural plants, such as Lilium [10], Oncidium [11], Anthurium [12], and so on. The MYB TFs play equally important roles in the early and late anthocyanin biosynthesis. However, the MBW complex, as a common regulatory mechanism, usually controls the expression of late biosynthesis genes. In gentian, co-expression of GtMYB3 and GtbHLH1 could induce the accumulation of anthocyanin by upregulating the expression of GtF3′5′H and Gt5AT [13]. In-depth exploration of the regulatory mechanism of TFs could contribute to a better understanding of anthocyanin biosynthesis in higher plants.

Anthocyanin biosynthesis in different plant species share conserved pathway and regulatory mechanism to some extent, but notable differences occur among species due to multiple genes involved and the complexity of the interactions between them. Since there is a lack of complete genome sequences and limited available molecular information, rare research on anthocyanin biosynthesis has been carried out in pansy. Thanks to high efficiency and reasonable expense of transcriptome sequencing (RNA-seq), it has been widely used to obtain information regarding gene expression and regulation [14]. Currently, RNA-seq has provided unique insights into molecular mechanisms and gene expression of non-model plants without reference genomes, including Dendrobium officinale [15], Senecio cruentus [16], Bougainvillea spectabilis [17], and so on. A series of biosynthetic genes and TFs involved in flavonoid synthesis have been thoroughly analyzed [18]. These studies lay the foundation for us to explore the mysterious world of anthocyanin biosynthesis in pansy.

Pansy (Viola × wittrockiana Gams.) is a biennial horticultural plant native to Europe with a wide range of flower colors. It has various cultivars with abundant pigmentations, including white, yellow, orange, red, purple, blue, and black. According to previous studies, the flower color variations are largely determined by different types of anthocyanins [19,20]. However, little research has focused on the molecular mechanism underlying anthocyanin synthesis in pansy. Our previous study on pansy showed that the content of anthocyanin in blotched petals accumulated rapidly from stage 3 (petals with light purple blotches) and reached the maximum, stage 5 (petals with dark purple blotches) [21] (Figure S1), indicating that stage 5 might be a vital stage to explore the mechanism of anthocyanin accumulation. Therefore, we conducted comparative transcriptome studies on non-blotched-flower (NBF) and blotched-flower (BF) at stage 5, using the biosynthetic genes and TFs related to anthocyanin biosynthesis. It provided a comprehensive overview of gene expression involved in the anthocyanin biosynthesis pathway. The findings will not only improve our understanding of underlying flower pigmentation in pansy, but also lay a foundation for further research on breeding.

2. Materials and Methods

2.1. Plant Materials

Viola × wittrockiana “Mengdie” was used in this study (Figure 1) and planted in the experimental station (located at 20.03° N/110.33° E) at Hainan University, Haikou, China. All seedlings were cultivated under natural light with day and night temperatures of 25–30 °C. In October, NBF and BF of three lower petals at initial flowering stage (S5) from three-month-old seedlings were collected for further analysis, including morphological observation, pigment determination, and RNA-seq (Figure S1). For morphological observation and pigment determination, the samples were cut and used immediately. For RNA-seq, at least 20 individual seedlings were randomly selected for sample collecting with two biological replicates. The samples were snap-frozen in liquid nitrogen and stored at −80 °C until use.

2.2. Flower Color Analysis

Petal colors were determined by using the Royal Horticultural Society Color Chart (RHSCC) at noon, indoors with north light. In addition, the color parameters were measured by the spectrophotometer NF333 (Nippon Denshoku, Tokyo, Japan) based on the CIE L*a*b* scale, including lightness (L*) and two chromatic components, a* and b* [22]. Three petals were randomly selected from three individual plants as biological replicates, and the parameters were measured five times as technical replicates.

Fresh petals were cross-sectioned by razor blades. The adaxial and abaxial epidermal layers were torn off by tweezers. They were placed on a glass slide with a drop of water and observed immediately under the light microscope DM2000 (Leica Microsystems, Wetzlar, Germany).

2.3. Determination of Anthocyanin and Carotenoid

The total anthocyanin content (TAC) of different parts of petals was assayed using the method presented by Li et al. [21]. Fresh petals were picked and crushed to obtain three replicates. Anthocyanins were extracted from petals (0.1 g) in methanol (containing 0.1% hydrochloric acid) at 4 °C overnight. The extracts were centrifuged at 10,000× g for 15 min. A total of 1 mL supernatants were diluted with 9 mL pH 1.0 buffer containing 25 mM KCl and 9 mL pH 4.5 buffer containing 400 mM sodium acetate, separately. The mixtures were balanced for 40 min at room temperature before absorbance detection at 530 nm and 700 nm by spectrophotometer. The total anthocyanin contents were calculated using the following equation:

$$\mathrm{Anthocyanin} \mathrm{content}(\mathrm{mg}/\mathrm{g})=\frac{\mathsf{\Delta }\mathrm{A}\times \mathrm{MW}\times \mathrm{DF}\times \mathrm{V}\times 1000}{\mathsf{\epsilon }\times \mathrm{L}\times \mathrm{M}}$$

(1)

ΔA = (A530 − A700) pH 1.0 − (A530 − A700) pH 4.5

where ΔA is the absorbance differences, MW is the molecular weight of cyanidin-3-glucoside (449.2 g·mol⁻¹), DF is the dilute factor (10), V is the supernatant volume (1 mL), ε is the molar absorptivity of cyanidin-3-glucoside (26,900 L·mol⁻¹·cm⁻¹), L is the path length (1 cm), and M is the sample weight (0.1 g).

The total carotenoid content (TCC) of different parts of petals was measured according to the protocol described by Li et al. with slightly modification [23]. Weights of 0.2 g petals were cut into pieces and soaked in 10 mL 80% acetone overnight in the dark for carotenoid extraction. The extractions were centrifuged and the supernatants were measured at 663 nm, 645 nm, and 470 nm by spectrophotometer. The total carotenoid contents were calculated using the following equation:

Carotenoid content (mg/g) = [(4.37 × A470 − 9.10 × A645 + 2.11 × A663) × V]/M × 1000

where V is the volume of supernatants (mL) and M is the usage of samples (g).

2.4. RNA Extraction and Sequencing

Total RNA was extracted from petals using Trizol reagent (SANGON Biotech, Shanghai, China) and treated with DNase I (TIANGEN Biotech, Beijing, China) according to Li et al. [21].RNA degradation and contamination was monitored on 1% agarose gels. RNA purity was checked by NanoPhotometer^® spectrophotometer (IMPLEN, Los Angeles, CA, USA). RNA concentration was measured by Qubit^® RNA Assay Kit in Qubit^® 2.0 Fluorometer (Life Technologies, Carlsbad, CA, USA). RNA integrity was assessed using the RNA Nano 6000 Assay Kit of the Agilent Bioanalyzer 2100 system (Agilent Technologies, Santa Clara, CA, USA). A total amount of 3 μg RNA per sample was used for library constructing. NEB-Next^®Ultra™ RNA Library Prep Kit for Illumina^® (NEB, Ipswich, MA, USA) was performed and index codes were added to attribute sequences to each sample. Sequencing was performed using the HiSeq™ 2000 platform at Novogene Bioinformatics Technology Co. Ltd. (Beijing, China).

2.5. Transcriptome Assembly and Function Annotation

Raw data were firstly processed through in-house perl scripts. At this step, clean data were obtained by removing reads containing adapter, reads containing ploy-N, and low-quality reads from raw data. In the meantime, Q20, Q30, GC-content, and sequence duplication level of the clean data were calculated. De novo assembly was accomplished using Trinity with min_kmer_cov set to 2 and all other parameters set by defaults [24]. Downstream annotation analysis was performed based on clean data with high quality to ensure the reliability of RNA-Seq.

The BLASTx with an expect E-value cutoff of 1 × 10⁻⁵ in the standalone NCBI-BLAST package 2.2.28+ was used to compare the unigenes with several databases, such as the Nr (NCBI non-redundant protein sequences), Nt (NCBI nucleotide sequences), and Swiss-Prot (a manually annotated and reviewed protein sequence database). The unigenes were also annotated using the Pfam (protein family) [25], and the results based on the NR and Pfam were imported into Blast2GO software v2.5 and in-house perl scripts for gene ontology (GO) term analysis [26]. The BLASTx results with a cutoff E-value of 1 × 10⁻³ were also aligned against the euKaryotic Ortholog Groups (KOG) database to predict the possible functions [27]. The KEGG pathways were assigned to the assembled unigenes using the online KEGG Automatic Annotation Server [28].

2.6. Identification of Differential Expression Genes (DEGs)

DEGs analysis was performed by DESeq2 using a model based on the negative binomial distribution [29]. The p-value was obtained to control the false discovery rate using Benjamini approach [30]. Gene expression level of each sample was calculated by the FPKM (fragments per kb per million reads) method using RSEM software [31]. The FPKM filtering cutoff of 1.0 was used to determine expressed transcripts. DEGs were screened based on the condition of FDR < 0.05 and |log2^Fold_change| ≥ 1. To better define the DEGs, we performed Gene Ontology (GO) enrichment analysis and Kyoto Encyclopedia of Genes and Genomes (KEGG) analysis [32]. We considered FDR < 0.05 to indicate significant enrichment.

2.7. Quantitative Real-Time PCR Validation

Six unigenes related to color development were selected for validation using qRT-PCR. The primers were designed using Primer Premier 5.0 software and are listed in Table S8. Total RNA was extracted using previous method and reverse transcribed into the first strand of cDNA using HiScript III All-in-one RT SuperMix Perfect for qPCR (Vazyme, Nanjing, China). The qRT-PCR was performed using ChamQ Universal SYBR qPCR Master Mix (Vazyme, Nanjing, China) on a LineGene 9600 Detection system (Bioer, Hangzhou, China). The reaction buffer (10 uL) consisted of 5 uL SYBR qPCR Master Mix, 1 uL gene-specific primers (0.5 uL for forward and reverse primer), and 25 ng cDNA template (ddH₂O as blank). β-actin was used as internal control to normalize gene expression [33]. The thermal cycles used were as follows: 95 °C for 60 s, followed by 40 cycles of 95 °C for 15 s, and 60 °C for 30 s. The melting curve analysis was carried out by adding the following procedure: 95 °C for 30 s, 60 °C for 30 s, and 95 °C for 30 s. Relative expression level of target genes was evaluated by the 2^−ΔΔCt comparative threshold cycle (Ct) method [34]. All qRT-PCR analyses were performed with three biological and technical replications.

2.8. Statistical Analysis

Different color parameters, TAC and TCC, were expressed as the mean ± standard deviation of three replicates. The data was examined via analysis of variance (one-way ANOVA). The statistics analysis was performed using SPSS software (version 25.0).

3. Results

3.1. Phenotype Differences between NBF and BF

The flower colors of pansy were defined by the Royal Horticultural Society Color Chart (RHSCC) (

Table 1. The color parameters of NBF and BF.

Sample	RHSCC Code	Color Group	CIE L*a*b*
			L*	a*	b*
NBF	9 A	Yellow	85.96 ± 1.90	−10.69 ± 0.23	105.11 ± 6.26
BF	187 A	Greyed-purple	2.56 ± 0.72 **	12.43 ± 0.97 **	−1.69 ± 0.14 **

** Indicates significant difference (p < 0.01).

). The color grades of NBF and BF were 187A and 9A, which were classified into two separate group, greyed-purple and yellow. It could be well distinguished according to the phenotype. Moreover, color parameters L*, a*, and b* were investigated, and significant differences were observed between NBF and BF (Table 1). High L* value of NBF indicated brighter color and lower value of BF indicated a darker color. Parameter a* shows red and green from positive to negative value. Therefore, in the BF, the higher a* indicated a reddish color. The parameter b* represented a negative value in BF and positive value in NBF, which represented blue color of BF and yellow of NBF. To some extent, it indicated that the positive value of a* and negative value of b* were caused by the accumulation of anthocyanin in the BF.

To explore the influence of pigment distribution on flower color, adaxial and abaxial epidermis of NBF and BF were observed under microscope (Figure 1). The adaxial epidermis cells of both NBF and BF were closely arranged and cone-shaped, and the abaxial epidermis cell shapes appeared rounded. In addition, the abaxial epidermis colorations of NBF and BF were close to yellow. However, the adaxial epidermis of BF showed deeper purple coloration than that of NBF. It was also clearly shown from the transverse sectioning.

3.2. Pigment Analysis and Correlation with Flower Color

Variation in TCC and TAC were observed between NBF and BF of pansy (

Table 2. Variance differences in TCC and TAC.

Sample	TCC Content (mg/g)	TAC Content (mg/g)
NBF	0.064 ± 0.013	0
BF	0.045 ± 0.019	2.099 ± 0.300 **

** Indicates significant difference (p < 0.01).

). The carotenoid content of BF was 0.045 ± 0.019 mg/g, which was slightly lower than that of NBF (0.064 ± 0.013 mg/g). No significant differences were found in carotenoid content in those two petals. On the contrary, obvious differences were detected in the anthocyanin content between BF (2.099 ± 0.301 mg/g) and NBF. Interestingly, no anthocyanin was found in the non-blotched petals (0).

Pearson correlation analysis was applied to investigate the relationship between flower color and pigment composition (

Table 3. Pearson correlation between color parameters and pigment content.

CIE L*a*b*	TCC	TAC
RHSCC	−0.597	0.987 **
L*	0.593	−0.988 **
a*	0.569	0.983 **
b*	0.617	−0.984 **

** Indicates significant difference (p < 0.01).

). The results showed significant correlations between TAC and RHSCC, L*, a*, and b*. The parameter of a* was positively correlated with TAC, whereas L* and b* were negatively affected by TAC. The RHSCC number increased while TAC increased. In contrast, no correlations were shown between TCC and these color parameters, indicating that the effect of TCC on variation of flower color was not important.

3.3. Characterization of Transcriptome Data and Gene Functional Annotation

Four cDNA libraries were constructed in pansy. In total, 0.43G and 0.40G were achieved in each sample. The error rate, Q20, Q30, and GC contents of two samples were within acceptable limits (Table S1). The mapped ratios of each sample were greater than 80%. The FPKM method was used to compare the expression level of genes between two samples. The FPKM distribution of all samples tended to be stable. (Table S2). The Pearson correlation coefficients (R²) between samples were both higher than 0.8 (Figure S2), which indicates that the transcriptome data was reliable for subsequent analysis. The clean data were stored in the Short Read Archive database of the NCBI (PRJNA819825).

De novo assembly was applied to construct transcripts with clean reads (Table S1). A total of 43,908 unigenes were gained after assembly by Trinity (version: v2012-10-05) software. The length of unigenes ranged from 201 bp to 8182 bp, with an average length of 682 bp. The N50 value was 1175 bp (

Table 4. Summary of RNA-sequencing and de novo assembly.

	Unigene	Transcript
Total number	43,908	118,747
Min length	201 bp	201 bp
Max length	8182 bp	8182 bp
Mean length	682 p	866 bp
N50	1175 bp	1312 bp
N90	269 bp	380 bp
Total nucleotides	29,961,359	102,873,692

). Therefore, given the accurately sequencing data, we could use it to perform functional analysis and DEGs detection.

Functional annotation of the assembled transcriptomes was performed in several public databases, including Nr, Nt, Pfam, KOG, Swiss-Prot, KEGG, and GO (Table S3). A total of 25,989 unigenes were successfully annotated. Among them, 24,323 unigenes (55.39%) could be annotated in the Nr database; 18,979 (43.22%) in the GO databases; 16,668 (37.96%) in Swiss-Prot database; 15,525 (35.35%) in the Pfam databases; 11,466 (26.11%) in the Nt database; 8507 (19.37%) in the KOG database; and 7680 (17.49%) in the KEGG database. Due to the lack of genome sequence in pansy, 17,919 unigenes did not match any of the databases, suggesting that they may be novel transcripts, small size of the sequence, noncoding RNAs, or misassembles [35].

3.4. GO and KEGG Classification of DEGs

Based on sequencing data and annotation analysis, DEGs were detected according to the standard, log2^{|Fold change|} > 1, and FDR < 0.05. In total, 79 DEGs were detected, including 40 upregulated and 39 downregulated DEGs (Table S4). The number of upregulated DEGs was similar to downregulated ones. The distribution of DEGs are shown in the volcano map (Figure 2).

To investigate the functional features of these DEGs, GO enrichment analysis was carried out (Figure 3). A total of 79 DEGs were classified into three categories, molecular function (MF), cell component (CC), and biological process (BP), among which BP was the largest group (Table S5). In the top 10 GO terms, MF and BP were the two main groups, and the DEG terms in BP were two times more than MF. The subcategories enriched in the MF group included “catalytic activity (GO:0003824)”, “binding (GO:0005488)”, and “transferase activity (GO:0016740)”. In the category of BP, “metabolic process (GO:0008152)”, “organic substance metabolic process (GO:0071704)”, and “primary metabolic process (GO:0044238)” were the top three largest subcategories, suggesting highly differential metabolic activities between NBF and BF. The upregulated DEGs in the top 10 GO terms were much higher than the downregulated ones. The types and proportions of subcategories in both upregulation and downregulation were similar to those in the GO enrichment of total DEGs. Most upregulated DEGs in the top 10 GO term were enriched in “metabolic process (GO:0008152)”. The unigenes involved in the flavonoid metabolic process (GO:0009812) were significantly upregulated, consistent with the phenotypic data.

KEGG functional analysis was conducted to obtain the difference in the metabolic pathway between NBF and BF. The results revealed that a total of 31 pathways were identified, of which 13 were downregulated and 19 were upregulated. Especially, “Flavone and flavonol biosynthesis (ko00944)” and “Flavonoid biosynthesis (ko00941)” were upregulated (Table S6). This provided reference for the comparative investigation between NBF and BF, particularly the molecular mechanism of color variation.

3.5. Genes Involved in Anthocyanin Biosynthesis

Based on the transcriptome data and DEG functional analysis, we found the key unigenes involved in the anthocyanin biosynthesis pathway and compared the expression level of those genes between NBF and BF (Figure 4). Almost all biosynthetic genes related to anthocyanin biosynthesis were analyzed, and most of them were downregulated in NBF (Table S7). From the global overview of the anthocyanin biosynthesis pathway, the expression level of genes in the early biosynthesis pathway (EBP) were commonly higher than those in the late biosynthesis pathway (LBP). These genes included phenylalanine ammonia-lyase (PAL), 4-coumarate CoA ligase (4CL), chalcone synthase (CHS), and flavanone-3-hydroxylase (F3H). However, no significant differences were shown among those genes between NBF and BF. To some extent, the synergistic expression of these biosynthetic genes promoted the accumulation of anthocyanins in BF. The expression of flavonoid-3′-hydroxylase (F3′H) was, quite, higher than flavonoid 3′,5′-hydroxylase (F3′5′H), and F3′5′H was significantly upregulated in BF. This may lead the metabolic flow to different anthocyanin synthesis. In addition, dihydroflavonol reductase (DFR) and anthocyanin synthase (ANS) were upregulated in BF, but the expression level was not high. The UDP-glucose flavonoid 3-glucosyltransferase (UFGT) was significantly increased and the expression level in BF was about 10 times higher in NBF.

To analyze the expression level of main TFs involved in anthocyanin biosynthesis, we applied the Plant Transcription Factor Database to screen out all TFs in pansy. A total of 754 TFs were gained, including 54 bHLH, 27 MYB, 47 MYB-related, and other transcription factors (Figure S3). MYB TFs play a vital role in anthocyanin biosynthesis. Among these MYB TFs, we found that one MYB (VwMYBP, comp50362_c0) was highly upregulated almost 60 times in BF. It belonged to the R2R3-MYB family, which was homologous to AtMYB82 (Figure S4). Moreover, qRT-PCR analysis was carried out to determine the expression pattern of VwMYBP in NBF and BF at six different developmental stages (Figure S5). Interestingly, no expression of VwMYBP was found at stage 1. The expression level of VwMYBP was significantly higher in BF than NBF and increased steadily from stage 2 to 5 and sharply reached the peak at stage 6.

3.6. Validation of RNA-Seq Results by Quantitative Real-Time PCR

Six unigenes were selected to verify the expression level of the transcriptome by gene-specific primers (Figure 5). In detail, one F3′5′H gene (comp54016_c0), four UFGT genes (comp71062_c0, comp60225_c0, comp37827_c1, comp37827_c0), and one MYB TF (comp50362_c0) were validated by qRT-PCR analysis. The expression levels of those six genes were consistent with RNA-Seq data, indicating that the transcriptome data were reliable.

4. Discussion

Flower color depends on the contents of the colored cells [36]. Here, in the NBF, both sides of epidermal cells showed a similar color of yellow. In contrast, the adaxial petal surface in the BF appeared to be distinctly purple, and the abaxial surface was populated with yellow cells. As expected, no anthocyanins were found in NBF compared to the significant abundant anthocyanins in BF. These findings were consistent with previous studies [21]. Previous studies have reported that anthocyanins act as main factors in flower coloration and have strong correlations with color parameters [1,37]. Significant correlations between TAC and flower color parameters were observed in pansy at the 0.01 level. With the increase in TAC, the value of a* increased and L* and b* decreased, indicating the thick red color in blotched petals. It led to the conclusion that anthocyanins play a vital role in flower pigmentation in pansy, which was in accordance with the findings of Li et al. [21]. Therefore, following research focused on the differences between anthocyanin biosynthesis in pansy.

Anthocyanin biosynthesis is relatively well understood in plants [38]. Due to the phenotypic and physiological differences, the identification of DEGs related to anthocyanin biosynthesis has indubitable superiority to illustrate the mechanism of color variation in pansy. Biosynthetic genes and TFs involved in anthocyanin biosynthesis have been identified in pansy [21]. The biosynthetic genes were divided into early biosynthesis genes (EBGs), such as CHS, CHI, and F3H, and late biosynthesis genes (LBGs), which include F3′H, F3′5′H, DFR, ANS, and UFGT [39,40]. Chalcone synthase (CHS) catalyzes the first reaction of anthocyanin biosynthesis and has been proved to be an important key enzyme in anthocyanin biosynthesis in different species [40,41]. The lack of expression of CHS directly caused the loss of anthocyanin [42,43]. In our study, we found that CHSs showed equally high expression in NBF and BF, but no significant differences between them. The expression patterns of CHI and F3H were similar to CHS. Therefore, pansies were different from other species in the formation of bicolor flowers. We speculated that the EBGs were not the main reasons for the lack of anthocyanin in NBF.

In the late biosynthesis pathway, several genes, including F3′H, F3′5′H, DFR, ANS, and UFGT, were upregulated in BF, among which F3′5′H and UFGTs were remarkably upregulated. F3′H and F3′5′H catalyze DHK to two intermediates, DHQ and DHM, respectively. DHQ and DHM are main precursors of two colored anthocyanins, cyanidin and delphinidin, separately [41,44]. In our previous studies, we identified two types of anthocyanins in BF, delphinidin and cyanidin, in which cyanidin accounted for 80% [21]. This might explain why the expression of F3′H was much higher than F3′5′H in BF. However, the intermediate dihydroflavonols were also involved in the flavone and flavonol biosynthesis [45]. Higher FLS expression in NBF might be accompanied with higher expression level of F3′H. This was consistent with the finding in alfalfa [41]. F3′5′H was known as the “blue gene” [46], and it was barely expressed in NBF. It was responsible for the dark color of blotched petals. Although no significant differences were found in DFR and ANS between two petals, it did not mean they were less important in color variation in pansy. On the contrary, they were highly expressed in the BF at stage 2–4 [21]. They directly determined the accumulation of anthocyanins and color formation at early blooming stage. Studies have shown that UFGT is a key enzyme for anthocyanin biosynthesis, and the activity is closely related with anthocyanin concentration and components [47,48]. The expression levels of UFGTs in BF were much higher than NBF (Table S7). We found that the cyanic blotches were legible and the content of anthocyanin reached the maximum at stage 5 [21]. UFGT increases the stability of anthocyanins by the glycosylation process [49]. This illustrated that UFGT might have an important role in the accumulation of anthocyanin and flower color formation at the late blooming stage. To sum up, the expression of LBGs would be likely to play a major role in color formation of pansy. This was consistent with the accumulation of anthocyanins in peony [50], pink camellia [51], and lily [39].

In addition to the biosynthetic genes, transcription factors (TFs) also play important roles in anthocyanin biosynthesis. Increasing numbers of TFs involved in anthocyanin biosynthesis have been identified, most of which belong to MYB families, and have regulated anthocyanin synthesis based on MYB-bHLH-WD40 complex [32]. This has been reported in Arabidopsis thaliana [52], petunia, and morning glory [53]. In the MYB family, R2R3-MYB TFs are the main regulators involved in anthocyanin biosynthesis [32]. In this study, through differential analysis, we found that VwMYBP (comp50362_c0) was an R2R3-MYB, which was homologous to AtMYB82. It was highly upregulated in BF at six development stages. MYB82 is a TT2-like transcription factor, which regulated the late biosynthesis genes by combining with bHLH and WD40 [32,54]. As per previous studies, the LBGs were normally regulated by MBW complex [53]. Thus, a potential relationship exists between VwMYBP and LBGs in anthocyanin biosynthesis; however, this speculation requires more investigation in further research.

5. Conclusions

Anthocyanin biosynthesis is a complex pathway in flower color regulation. In this study, phenotypic, physiological, and transcriptomic analysis of two types of petals revealed the mechanism of color variation in pansy. High content of anthocyanin contributes to dark purple color in blotched petals. It was consistent with the epidermal observation. High expression level of LBGs and TFs involved in anthocyanin biosynthesis increased the accumulation of anthocyanin in BF. This work will be helpful for further functional analysis and reveals the mechanism of color formation in pansy. It is probably more complex than these findings here since the material is limited to only one stage. Thus, it will be an interesting and challenging subject for more comprehensive future research to illustrate the flower color formation in pansy.