Transcriptomic and Chemical Analyses Reveal the Hub Regulators of Flower Color Variation from Camellia japonica Bud Sport

: Camellia japonica is a woody ornamental plant with multiple flower color variations caused by bud sport; however, the molecular mechanism remains unclear. Here, chemical and transcriptomic analyses of C. japonica were performed with white, pink, red, and dark red flowers caused by bud sport. Seven anthocyanins were detected in these samples, except in C. japonica ’YuDan’ (white petals). The total anthocyanin content of C. japonica ’JinBiHuiHuang’ was the highest, and cyanidin 3- O - β glucoside (Cy3G) was the main anthocyanin affecting the redness of petals. Furthermore, the ratio of Cy3G and cyanidin-3- O -(6- O -( E )- p -coumaroyl)- B -glucoside) was significantly correlated with the red petal phenotype. In total, 5673 genes were identified as differentially expressed genes (DEGs). The potential co-expression modules related to anthocyanin accumulation were established, which fea-tured transcription factors, anthocyanin biosynthesis, and plant hormone signal transduction. Thir-teen structural genes in the anthocyanin biosynthetic pathway were identified as DEGs, most of them were upregulated with deepening of flower redness. An integrated promoter and cluster analysis suggested that CjMYB62 , CjMYB52 , and CjGATA may play important roles in anthocyanin accumulation. These results provide insight and candidate genes for the transcriptional mechanism responsible for the bud sport phenotype.


Introduction
Bud sport is a common somatic mutation that occurs in leaves, flowers, and fruits [1] and leads to qualitative and quantitative changes in plant phenotypes [2] but retains the original qualities of the parent plant. Thus, many new plant cultivars have been obtained by breeding mutant buds with good characters. The mutants caused by bud sport often display changes in the anthocyanin content of different tissues. The bud sport induces changes in phenotypic colors in periclinal chimeras of carnation [3].
Bud sport contributes to different expression levels of genes in the anthocyanin synthesis pathway between red and white flowers of peach trees [4]. In addition, the induction of anthocyanin biosynthesis in grapevine fruit [5] and the mutations in some branches lead to a loss of color in leaves and stems in pear [6]. However, the transcriptional mechanism remains unclear.
Pigments are the most critical factor affecting the floral color. Pigments are divided into three groups, including flavonoids, betalains, and carotenoids. Anthocyanins are mainly the red-colored flavonoids in petals. The floral color turns white when anthocyanins degrade, and the floral color turns red when anthocyanins accumulate [7].
It has been reported that anthocyanin contents play the main role in the floral color of herbaceous peony cultivars [8]. In addition, seven anthocyanins (cyanidin 3-O-β-glucopyranoside, and others) have been isolated from the red petals of C. saluenensis [9]. The anthocyanin synthesis pathway involves the regulation of several genes [10], including PAL, C4H, 4CL, CHS, and CHI, which promote naringenin synthesis.
The formation of various flavonoids is regulated by the FLS, FNS, F3H, F3′H, and F3′5′H genes, and DFR, ANS, and LAR promote anthocyanidin biosynthesis. Moreover, UDP-glucose transferase further stabilizes anthocyanin by glycosylation. Interestingly, some transcription factors (TFs) (MYB, bHLH, and MADS-box) regulate the structural genes in the anthocyanin biosynthesis pathway by binding to the cis-acting element [11]. Although studies on flower color have been carried out, few are available on changes in anthocyanin content caused by bud sport.
Camellia japonica is an important ornamental flowering plant, and the flower color is a major characteristic of ornamental value. Due to a lack of molecular resources, the breeding of new cultivars mainly depends on bud sport and cross-breeding. Bud sport causes rich flower color in C. japonica [9]. C. japonica is considered an excellent model plant for research on flower color, as flowers range from white to dark red in C. japonica 'ChiDan'.
Initially, two novel C. japonica cultivars with pink or dark red petals were obtained by budding sport (C. japonica 'FenDan' and C. japonica 'JinBiHuiHuang'). Then, a novel C. japonica cultivar with white petals was obtained by budding sport (C. japonica 'YuDan') (Figure 1A). This greatly improved the ornamental value of C. japonica. The bud sport transcriptional regulatory mechanism is of great significance for understanding the different C. japonica color phenotypes; however, the molecular mechanisms have not been adequately explored.
To understand the molecular mechanism of flower color variation from C. japonica bud sport, seven anthocyanins were identified from four C. japonica cultivars formed by bud sport using high-performance liquid chromatography (HPLC). Correlation analysis of petal color and pigment content was performed. Moreover, the traditional divergence of anthocyanin accumulation in C. japonica was inferred by RNA-seq, which provided candidate genes involved in the bud sport floral color phenotypes and valuable information for the color diversity of C. japonica. The total anthocyanin content of the samples; the data are the mean ± SD of three independent replicates, different letters indicate a significant difference between groups using Tukey's Honest Significant Difference test at p < 0.05 after one-way ANOVA.

Plant Materials and Chroma Measurements
All camellias used in this experiment were obtained from the Germplasm Resource Center of the Institute of Subtropical Forestry, Chinese Academy of Forestry (Daqiao Road, Hangzhou City, Zhejiang Province). The annual rainfall is approximately 1500 mm, and the frost-free period is approximately 320 days. The soil at the test site is sandy loam soil, and the soil pH is 5.5-6.5. The irrigation and drainage conditions are good.
C. japonica were used as the control (CJ), C. japonica 'ChiDan' with red petals (CD), C. japonica 'JinBiHuiHuang' with dark red petals (JH), C. japonica 'FenDan' with pink petals (FD), and C. japonica 'YuDan' with white petals (YD) were selected as the experimental subjects, and their petals were collected. The camellias used in the experiment were all 6 years old. The petals were divided into five groups according to color and species, with three biological replicates in each group.
The chroma values were measured using a spectroscopic colorimeter (NF555, Spectrophotometer, TRISUM) according to the CIE L* a* b* method. Position was measured at the central part of the upper epidermis of the petals, and the light source was C/2°. Each measurement had three biological replicates.

Extraction and Quantitative and Qualitative Analyses of the Anthocyanins
Fresh petals (2 g) were ground in liquid N2. Flowers were extracted in buffer (methanol: water: formic acid: trifluoroacetic acid = 70:27:2:1, v/v) for 24 h under a shaded condition. Then, the supernatants were filtered through a 0.22 µm membrane. Three biological replicates were collected and stored at −20 °C.

Library Construction for Next-Generation Sequencing (NGS) and RNA Extraction
The 15 samples of total RNA from the five groups (CJ, CD, JH, FD, and YD) were extracted using the DP441 plant kit (TIAGEN, Beijing, China), following the manufacturer's instructions. Each treatment was comprised of three replications. RNA quality was determined using the NanoDrop1000 (ThermoFisher, Waltham, MA, USA) and Agilent 2100 instruments (Agilent Technologies, Palo Alto, CA, USA). A NGS RNA library was constructed. Oligo (dT) magnetic beads were used for mRNA enrichment by AT complementary pairing with the poly-A tail of the mRNA. Polymerase chain reaction (PCR) enrichment was performed to obtain the cDNA library [15], and then libraries were sequenced on an Illumina HiSeq 2000 platform. All the clean reads were assembled using a de novo assembly program Trinity.

Identification of DEGs and the Weighted Gene Co-Expression Network Analysis (WGCNA)
The expression levels of the transcripts were quantified based on the location of the mapped reads in the transcripts [24] and were calculated using the Fragments Per Kilobase of transcriptome per Million mapped reads method. DESeq [25] was used for the differential expression analyses between the sample groups. The DEG screening conditions were fold change ≥2 and a false discovery rate (FDR) < 0.01. Genes with the highest contribution to 25% of the variance were set as targets, and the R package WGCNA (The R Foundation for Statistical Computing, Vienna, Austria) was used to identify co-expression modules related to anthocyanin accumulation.

Cloning and Cis-Element Analysis of the Promoter
DNA was extracted using the CTAB method [26]. According to the genome sequences of related species, PCR-specific primers were designed with Primer 5.0 software (Supplementary Table S1). The sequence was connected to the 18 T vector for Sanger sequencing. The cis-acting elements were identified using PlantCARE [27].

qRT-PCR Analysis
Primer Premier v6.0 was used to design the gene-specific primers (Supplementary Table S1). The RNA extracted from the sequencing samples was reverse transcribed to cDNA using the PrimeScript RT Reagent Kit + gDNA Eraser (Takara, Shiga, Japan). The qRT-PCR reactions were carried out according to the instructions for TB Green ® Premix Ex Taq II (Takara, Shiga, Japan) using the ABI 7500 Fast Real-Time PCR System (ABI, Foster City, CA, USA) instrument. Gene expression was quantified by the 2 −△△CT method, and 18 s rRNA was used as the internal reference gene.

Statistical Analysis
All data were analyzed with three biological replicates. The statistical analysis was conducted using SPSS (v22.0) software (SPSS Inc., Chicago, IL, USA). A p-value < 0.05 was considered significant between mean values based on one-way analysis of variance and Tukey's test. The data are presented as the mean ± standard deviation of three independent experiments.

The Chroma Values of the Different Petals
The chroma value is used to scientifically represent the floral color phenotype. The a* value indicates redness and greenness, ranging from −0.74 (YD) to 53.63 (CJ), whereas the b* value indicates blueness and yellowness, ranging from 1.66 (FD) to 18.28 (CJ). The L* value indicates lightness, ranging from 35.21 (JH) to 91.49 (YD) (Supplementary Table  S2). The color index of these samples was significantly different, which was consistent with the PCA results for the chroma values ( Figure 1B).

Qualitative and Quantitative Analyses of the Anthocyanins
Seven anthocyanin peaks were detected in the samples, except in YD ( Figure 2A). The anthocyanins were identified by referencing published data of retention time, elution order, UV-vis spectrum, and the MS fragmentation pattern. The absorption peak of the UV-vis spectrum (513-520 nm) and the fragment ions (287 u) showed that the seven anthocyanins were cyanidin-based [28]. The A440/Avis-max ranged from 32% to 34%. All anthocyanins were the 3-O-glycoside type [29]. The molecular ions P1 and P2 were 449 u, and the difference between m/z 449 u to m/z 287 u was 162 u, e.g., glucoside ( Figure 2B).
The elution time of anthocyanin galactose is less than that of anthocyanin glucoside. Overall, P1 and P2 were tentatively identified as cyanidin 3-O-β-glucoside and cyanidin 3-O-β-galactopyranoside, respectively. Similarly, five anthocyanins were established, in- The total anthocyanin content was the highest in JH, which was more than 29 times that in FD and six times that in CD. The anthocyanin content in YD was not measured ( Figure 1C and Supplementary Table S3), which was consistent with the chroma measurements. Interestingly, despite the total anthocyanin content, Cy3G was the main contributor to red petals ( Figure 2C).
Furthermore, the correlation analysis of seven anthocyanins was performed, indicating that the Cy3G and Cy3GaEpC contents were positively related (R = 0.74) ( Figure 2D). In addition, the proportion of different anthocyanins is another important factor affecting flower color; the redness of the petals deepened along with the increase in the proportion of Cy3Ga and Cy3G. Moreover, the proportion of Cy3GaEpC and Cy3GEpC decreased along with the deepening of the red color and was the lowest in JH ( Figure 2E).

Overview of the Transcriptome Data
To further study the molecular mechanism of anthocyanin biosynthesis in camellia bud sport cultivars, a transcriptome analysis was performed based on the Illumina sequencing platform. Fifteen samples were divided into five groups, including CJ, YD, FD, CD, and JH. The number of raw reads in each library ranged from 22 to 28 million. After the low-quality reads were removed, and the high quantity sequences were trimmed, 107.23 Gb of clean data were generated with Q30 ≥ 92.50% and a GC content of 45.15-46.13%.
The average reads of all samples were 7.14 Gb of data, the average number of reads was 25,657,417 (Supplementary Table S4). The sequencing data were sufficient for further analysis. In total, 102,021 transcripts and 44,356 unigenes were generated using the Trinity sequence assembler. The average length of a transcript was 1136 bases, and the mapped ratio was 59.29%. The N50 length was 1638 bp.

Identification of DEGs among Camellias of Different Flower Colors
The screening conditions were set to fold change ≥2 and an FDR < 0.01 to identify the DEGs. The CJ group was the control group, and 5673 genes were identified as DEGs. We identified 464 DEGs specific to the CJ vs. CD comparison; 327 DEGs were specific to the CJ vs. FD comparison; 352 DEGs were specific to the CJ vs. JH comparison, and 510 DEGs were specific to the CJ vs. YD comparison.
A total of 1653 DEGs were specific in the four groups. In addition, 1369 DEGs were shared by the four groups. The maximum number of DEGs (1531) was found in the JH vs. YD comparison ( Figure 2B), suggesting that the number of DEGs was positively correlated with the color difference. The heatmap of all DEGs is shown in Figure 2C. These genes may play an important role in petal color mutations, indicating that the transcriptional mechanism of flower color change is a complex integrated network in C. japonica.
To explore the molecular mechanism for the formation of different flower colors (dark red, red, pink, and white) in C. japonica, the functional categories of the DEGs were classified into three groups (cellular components, molecular functions, and biological processes) in GO analysis. The first 17 enriched GO terms of the CD vs. JH, FD vs. YD, and CD vs. FD comparisons are shown in Figure 3B and suggest that the DEGs were primarily enriched in "catalytic activity", "protein phosphorylation", "regulation of transcription", "transferase activity", "integral components of the membrane", and "binding" ( Figure 3A and Supplementary Table S5).
These results indicate that TFs play an important role in the flower color mutation. A KEGG enrichment analysis was performed, and the top 19 pathways are shown in Figure  3B (Supplementary Table S6). The results show that the phenylalanine metabolism and flavonoid biosynthesis were significantly enriched, which are key pathways related to the formation of red petals.
Plant hormone signal transduction and glycolysis/gluconeogenesis metabolism were significantly enriched, indicating that hormones and glycosylation play a key role in white petals turning red. Glutathione metabolism and oxidative phosphorylation were other terms in the analysis. Overall, the results indicate that the formation of red petals involves complex transcription and metabolic pathways.

Identification and Cluster Analysis of the Transcription Factors
TFs play a key role in regulating gene expression [30]. Here, 774 transcripts were annotated as TFs. MYB was the most abundant TF family (190 TFs), followed by the bHLH family (80 TFs) (Supplementary Table S7). A total of 83 differentially expressed transcription factors (DETs) were identified in the RNA-seq data. The number of MYB-like TFs was 21. These results indicate that MYB TFs play an important role in the flower color mutation. A cluster analysis of the 87 DETs was performed using the R package 'Mfuzz'. The number of clusters was set to six ( Figure 3E), and the results showed that the number of genes in cluster 1 was the highest.
The expression patterns of clusters 1 and 6 were significantly correlated with flower color. Twenty-one TFs in cluster 1 and 12 TFs in cluster 6 were upregulated and downregulated along with anthocyanin accumulation, respectively ( Figure 3E). The top three gene families in Cluster 1 were "MYB", "bHLH", and "HSF". The top three gene families in Cluster 6 were "MYB", "Homeobox", and "HSF".

Identification of Candidate Genes Involved in Anthocyanin Biosynthesis
The anthocyanin biosynthetic pathway mainly involves the flavonoid synthetic pathway and phenylalanine metabolism. Our results show that the DEGs participating in the anthocyanin biosynthetic pathway changed as shown in the pattern diagram ( Figure 4). The RNA-seq results revealed that 13 structural genes were identified as DEGs. Based on the chemical data, these genes were related to cyanidin biosynthesis, and four CjPAL (c114370, c110860, c108794, and c126434) genes, Cj4CL (c148171), two CjDFR (c142804 and c134588) genes, and one CjLAR (c125393) gene were significantly differentially expressed, which decreased along with the fading of flower color.
In addition, CjF3H (c136386) was upregulated in the CD vs. JH comparison, whereas CjLAR (c125393) was downregulated in the CD vs. YD comparison. These genes might be important structural genes that affect anthocyanin biosynthesis. Moreover, 11 CjGSTs and five CjABCCs were identified. These genes participate in anthocyanin transport, and most were highly expressed in JH. The degradation of anthocyanins plays an important role in flower color variations, and the related genes were identified, including CjLAC and CjPPO. Most of these genes decreased as anthocyanin accumulated. Overall, these genes correlated with the flower color variations.
To further explore the candidate genes related to anthocyanin synthesis, we cloned and analyzed the promoters of the above genes, and the MYB, MYC, and GATA binding sites were identified (Supplementary Table S8). CjMYBs (c118169.graph_c0, c144149.graph_c3, and c144055.graph_c0) and Cjzf-GATA (c122187.graph_c0) were identified as candidate TFs involved in the anthocyanin pathway according to the integrative promoter and correlation analyses, indicating that TFs plays an important role in anthocyanin biosynthesis (Supplementary Figure S1). Furthermore, methyl jasmonate (MeJA)-responsiveness, the auxin-responsive element, and the anaerobic induction element were identified, providing evidence for hormones participating in anthocyanin biosynthesis. For the biosynthesis of salicylic acid (SA), two CjTGA (c128839.graph_c1 and c89768.graph_c0) genes and one CjPR (c124349.graph_c0) gene increased with deepening of the flower color.

The Co-Expression Modules Related to High-Anthocyanin Accumulation in Petals
A WGCNA was performed to identify the hub network related to high anthocyanin accumulation. A total of 1418 genes made the largest contributions to 25% of the variance and were screened from 5673 DEGs to construct the modules. Power 18 was settled to obtain relatively balanced scale independence and connectivity of the results. Eight modules (brown, turquoise, red, green, yellow, black, blue, and grey) were identified ( Figure  5A).
The TOM plot suggested almost no correlation among these modules but a high correlation within modules ( Figure 5B). Thus, our results are reliable. The yellow module was significantly positively correlated with JH (the highest anthocyanin content in the five samples) (r = 0.97, p = 4 ×10 −8 ; Figure 5C). In contrast, the negative correlation between the yellow module and petal redness increased as anthocyanin content decreased, indicating that the yellow module contained a hub genetic component related to anthocyanin accumulation. In addition, one gene (c136386.graph_c0) from the flavonoid biosynthetic pathway and six genes (c135611.graph_c0, c144556.graph_c0, c128992.graph_c0, c134612.graph_c0, c129739.graph_c0, c120831.graph_c0, and c137350.graph_c0) involved in the phenylpropanoid biosynthetic pathway were highly enriched in the yellow module, indicating that these genes are possible regulators of high anthocyanin accumulation.
To further explore the traditional mechanism of anthocyanin accumulation, the promoters of some genes were cloned to predict the cis-acting element and identify the binding site of traditional factors. The results suggested some of the binding sites, such as the MYB and MYC binding sites (Supplementary Table S8). CjMYB62 (c118169.graph_c0), CjMYB1R1 (c118516.graph_c1), CjMYBs (c122246.graph_c0, c132870.graph_c0, and c144055.graph_c0), CjMYB52 (c144149.graph_c3), and CjMYB142 (c117274.graph_c0) were tentatively screened by integrative binding site and expression pattern analyses. Finally, the potential co-expression network related to anthocyanin accumulation was constructed using genes from the yellow module ( Figure 6C). Pearson's correlation analysis of the candidate genes and TFs (D) A potential DEG regulatory network with TFs related to anthocyanin accumulation in the yellow module. Co−expressed genes with weights >0.5 were screened as candidate genes, and TFs were screened via integrative cis-acting element and expression pattern analyses.

Quantitative Real-Time PCR Verification of the Transcriptome Data
To check the validity of the RNA-seq results, the expression patterns of six candidate genes were determined in all samples by qRT-PCR (Supplementary Figure S2). The expression levels of starch synthase 4 (c137365) and UDP-glucosyl transferase 74B1 (c139663) decreased dramatically along with the reduction of redness in the petals, indicating that starch synthase 4 and UDP-glucosyl transferase 74B1 regulate the biosynthesis of anthocyanins by affecting glycosylation. Similarly, aquaporin protein 9 (c119563) was expressed at the highest level in C. japonica 'JinBiHuiHuang'. In addition, the expression level of GATA TF 24 increased along with anthocyanin accumulation. The qRT-PCR results were consistent with the RNA-seq data, demonstrating that our transcriptome data of C. japonica petals are reliable.

Discussion
Bud sport is a common reason for color variations in plants. Most azalea varieties are color sports, which may be the result of transposon-mediated changes in pigment genes [3]. A series of bud sport cultivars from C. japonica, e.g., JH, CD, FD, and YD, with colors ranging from white to dark red, provided the material to explore the bud sport phenotype. Flower color is an important ornamental evaluation index. The contents and components of pigments are key factors that affect the brightness of petals [31].
Flavones and flavonols provide yellowness [32], while anthocyanins are correlated with redness [33], such as the red leaves in Eucommia ulmoides [34], purple-fleshed sweet potato [35], and grape berries [36]. In our results, red petals formed due to anthocyanin accumulation. The anthocyanin content was positively correlated with the a* value within a certain range (0-40.6).
Interestingly, the continual increase of anthocyanin content in JH significantly raised the h* value, but the a* value did not change significantly. These results are consistent with a previous study reporting that reducing anthocyanin content decreases the a* value [37], and the h* value describes seven color tones (red, orange, yellow, green, cyan, blue, and purple) [38].
Six common anthocyanins are pelargonidin, cyanidin, delphinidin, peonidin, petunidin, and malvidin [39]. Cyanidin correlates with redness, and delphinidin contributes to a blue-purple color [40]. Our biochemistry data suggested that only cyanidins were detected in the C. japonica samples, corroborating our previous result [13], e.g., the red petals of C. japonica mainly contain cyanidins. Furthermore, the stability of anthocyanins increases by binding to glycosides [41]. Here, seven cyanidin molecular structures were identified, which were involved in the binding of glucoside, glucopyranoside, caffeoyl, and coumaroyl.
Cy3G contributed the most to the pink, red, and crimson phenotypes in all C. japonica used in the experiment. The proportion of Cy3Ga in CD was not significantly different from that in FD but was significantly enhanced in JH, indicating that Cy3Ga plays an important role in the formation of the crimson phenotype. The proportion of Cy3GaEpc and Cy3GEpc decreased along with the deepening of petal redness. We speculate that the importance of Cy3GaEpc and Cy3GEpc is lower than that of Cy3Ga and Cy3G in the crimson phenotype. This result was consistent with a previous study reporting that Cy3G is the main anthocyanin in the red petals of peach [42].
The anthocyanin biosynthetic pathway was reported [43], and a transcriptomic analysis showed that 13 candidate structural genes involved in anthocyanin biosynthesis were identified as DEGs. Most of these genes were downregulated along with a reduction in anthocyanin content. Cyanidin-like glucoside biosynthesis in camellias occurs as follows. The first step is phenylalanine ammonia-lyase (PAL), cinnamate-4-hydroxylase (C4H), and 4-coumarate-CoA ligase (4Cl) promote the conversion of phenylalanine to coumarate-CoA [44].
In our results, four PAL and one 4CL gene were identified as DEGs. Unlike previous studies, C4H is not differentially expressed in Malus halliana [45], indicating that CjPAL and Cj4CL are the main DEGs affecting phenylalanine biosynthesis in C. japonica. The second step is that chalcone synthase (CHS) and chalcone isomerase (CHI) promote the conversion of 4-coumaroyl-CoA to naringenin and flavanone 3-hydroxylase (F3H) promotes the synthesis of dihydrokaempferol. The transcriptome analysis of Paeonia petals revealed that the PsCHS and PsF3H genes are differentially expressed and related to flower color [46].
These results combined with our chemical data indicate that the expression pattern of these DEGs was consistent with the seven types of cyanidin synthesis, particularly that of Cy3G. Furthermore, CjANR and CjLAR were significantly downregulated in C. japonica 'YuDan'. Catechin and epicatechin are often considered colorless substances, while previous studies have shown that catechins affect the transcriptional balance of F3′H and F3′5′H in the white leaves of Camellia sinensis [32], and higher concentrations of catechins accelerate the anthocyanin coloring reaction [47].
In addition, catechins have strong antioxidant properties [48]. Our GO enrichment analysis of the DEGs revealed that oxidation reactions play an important role in flower color variations. GST may be involved in anthocyanin transport in Petunia hybrida [49].
We found that the expression levels of seven CjGSTs and two CjABCCs were significantly higher in JH, indicating that these genes promote the transport of anthocyanins.
The degradation of anthocyanins affects their content [50], which involves laccase (LAC) [51]. In our results, the expression patterns of three LAC and four PPO coding genes were negatively correlated with the deepening of the red flower color, as these genes are critical regulators of anthocyanin degradation.
Although hormones play an important role in flower development, their effect on flower color is largely unknown. In our results, some genes involved in ABA, indole acetic acid (IAA), and SA signal transduction were identified as DEGs, including CjSAUR, CjAUX, CjTGA, CjGH3, and CjPP2C. High expression levels of auxin-responsive genes promote the biosynthesis of anthocyanins [52].
The expression level of CjAUX increased in the crimson petal phenotype. Treating green fruit with exogenous ABA promotes anthocyanin accumulation [53]. Similarly, we observed that one CjPP2C gene was downregulated as the flowers faded. The SA signaling pathway mediates anthocyanin synthesis in A. thaliana [54]. The changes in genes in the hormone signal transduction pathways may regulate the variations in flower color.
TFs play important roles in regulating anthocyanin biosynthesis [55]. Many studies have shown that MYB-like transcription factors directly regulate anthocyanin biosynthesis [56]. AtMYB12 is a flavanol-specific regulator of phenylpropanoid biosynthesis in A. thaliana [57]. The expression level of MhMYB10 is downregulated during the flower color fading period and is related to the expression of anthocyanin biosynthetic genes in the petals of Malus halliana [45], thus, indicating its important role in the anthocyanin biosynthetic pathway by regulating the expression of MpDFR and MpANS in Malus pumila [58].
Our results show that 21 MYB TFs were identified as DETs, and the expression levels of six CjMYBs were upregulated with deepening flower color (CjMYB62, CjMYB1R1, and CjMYB85). Furthermore, a previous study showed that bHLH TFs regulate ANS expression by forming the MBW complex [59]. Here, cluster analysis showed that the expression patterns of four CjbHLHs (CjbHLH51, CjBEE3, CjbHLH75, and CjbHLH117) were significantly correlated with anthocyanin accumulation. bHLH1 and bHLH2 control different pigments pattern by regulating structural genes of anthocyanin biosynthesis [60]. These results reveal that CjMYB and CjbHLH TFs play important roles in regulating anthocyanin biosynthesis in C. japonica.
High anthocyanin accumulation involves different metabolic pathways, such as glycosylation and oxidation [61]. In this study, the regulatory network for the formation of the crimson phenotype in JH was predicted. MYB and GATA interacted with the binding sites on the promoters of genes related to phenylpropanoid synthesis. Furthermore, ARE, MeJA, light, and IAA were also identified. Auxin inhibits the biosynthesis of anthocyanins in apple [52], and anthocyanin content increases after a MeJA treatment [62]. Overall, anthocyanin accumulation may be regulated by hormones and environmental factors.

Conclusions
In summary, the chemical and transcriptomic analyses provided evidence for the pigment composition and genes related to anthocyanin biosynthesis in C. japonica. Among the seven anthocyanins identified, Cy3G was the main pigment. The 13 structural genes (e.g., CjPAL, Cj4CL, CjF3′H, CjDFR, CjANR, and CjLAR) and 33 TFs (MYB, HSF, bHLH, GATA, and others) were hub regulators affecting the variations in flower color.
In addition, a predicted potential module (Figure 7) explained the high anthocyanin accumulation in C. japonica 'JinBiHuiHuang', indicating that increased phenylpropanoid promoted the deepening of redness in petals and that a high expression level of beta-glucosidase promoted anthocyanin glycosylation. This genetic information will be helpful for breeding studies of C. japonica floral color. Supplementary Materials: The following supporting information can be downloaded at: www.mdpi.com/article/10.3390/horticulturae8020129/s1. Figure S1. The correlation analysis of TFs and candidate genes, Figure S2. The expression level analysis of RT-qPCR and transcriptome, Table  S1. The primer of promoter cloned and RT-qPCR, Table S2. The raw data of chroma values, Table  S3. The raw data of anthocyanin content, Table S4. Summary of sequencing quality, Table S5. The GO enrichment analysis of DEGs, Table S6. The KEGG enrichment analysis of DEGs, Table S7. The classification of transcription factors, Table S8. The identified cis-acting element.

Acknowledgments:
The authors are grateful to the editors and referees for their valuable comments to improve our manuscript. In addition, Menglong Fan wants to thank the invaluable support received from Rui Zhou.

Conflicts of Interest:
The authors declare no conflict of interest.