{'fraction_non_alphanumeric': 0.046831955922865015, 'fraction_numerical': 0.0380222681359045, 'mean_word_length': 5.110643520953884, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 7, 'https://': 1, 'lorem ipsum': 0, 'www.': 3, 'xml': 0}, 'pii_count': 2, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Integrated transcriptome and metabolome profiling of Camellia reticulata reveal mechanisms of flower color differentiation

Camellia reticulata (Lindl.) is an important ornamental plant in China. Long-term natural or artificial selections have resulted in diverse phenotypes, especially for flower colors. Modulating flower colors can enhance the visual appeal and economic value in ornamental plants. In this study, we investigated the molecular mechanisms underlying flower color differentiation in C. reticulata. We performed a combined transcriptome and metabolome analysis of the petals of a popular variety C. reticulata (HHYC) (red), and its two cultivars “Xuejiao” (XJ) (pink) and “Tongzimian” (TZM) (white). Targeted metabolome profiling identified 310 flavonoid compounds of which 18 anthocyanins were differentially accumulated among the three samples with an accumulation pattern of HHYC > XJ > TZM. Likewise, transcriptome analysis showed that carotenoid and anthocyanin biosynthetic structural genes were mostly expressed in order of HHYC > XJ > TZM. Two genes (gene-LOC114287745765 and gene-LOC114289234) encoding for anthocyanidin 3-O-glucosyltransferase are predicted to be responsible for red coloration in HHYC and XJ. We also detected 42 MYB and 29 bHLH transcription factors as key regulators of anthocyanin-structural genes. Overall, this work showed that flavonoids, particularly anthocyanins contents are the major determinants of flower color differentiation among the 3 C. reticulata samples. In addition, the main regulatory and structural genes modulating anthocyanin contents in C. reticulata have been unveiled. Our results will help in the development of Camellia varieties with specific flower color and quality.


Introduction
Camellia reticulata (Lindl.), also known as Yunnan camellia is among the decorative trees belonging to the Camellia genus, which are members of the Theaceae family (Xia, 2003;Xin et al., 2015). There are 97 different species of Camellia, 76 of which are endemic to China (Yu, 1985;Xia, 2003;Zhang et al., 2019). China is considered the origin and distribution hub of Camellia. Generally, members of the genus Camellia are classified into one of five groups. Yunnan camellia (C. reticulata) and its near relatives constitute a significant group that is primarily found in Yunnan Province. In Yunnan province, Camellia is regarded as an exceptional ornamental plant with high economic values. Some of them are very important flowering ornamentals that also produce oil, and have a wide variety of cultivars (Silan et al., 2013;Xin et al., 2015). In addition to its cultural and medicinal significance, C. reticulata is also a valuable source of food products (Xia, 2003;Silan et al., 2013;Xin et al., 2015). C. reticulata is mostly red-colored, and there are few accounts of color variations in the species. As a result, producers and researchers have constantly sought to breed C. reticulata of diverse colors to enhance its aesthetic value (Xin et al., 2015).
Flower color is a key feature of horticultural crops, directly linked to pigments' content and distribution (Zhao and Tao, 2015). Numerous studies have revealed that carotenoids, flavonoids, and alkaloids are among the most important pigments that contribute to the development of flower colors (Wang et al., 2001;Ma et al., 2018;Zhu et al., 2019). Particularly, carotenoids and flavonoids have garnered increasing interest as they are pigments extensively dispersed in plants and are responsible for the spectacular hues found in petals (Mekapogu et al., 2020;Cui et al., 2021;Gao et al., 2021;Mei et al., 2021). Carotenoids can cause flowers to display a wide spectrum of colors, from vivid red to orange to yellow (Clotault et al., 2008;Tanaka et al., 2008;Zhao and Tao, 2015). Likewise, flavonoids, which are the most significant secondary metabolites in plants, can create colors ranging from light yellow to purple, depending on the type of flavonoids. Camellia blooms have been the subject of recent studies due to the plant's status as an ornamental plant. Carotenoids and flavonol glycosides were revealed to be the primary pigments in the golden yellow blooms of the C. nitidissima by Zhou et al. (2017). The pigment composition of red-flowered Camellia species was investigated by Jianbin Li et al. (2013), who concluded that cyanidin-core structure pigments (such as cyanidin 3,5-di-O-glucoside) have the potential to produce the most dominant phenotype among wild red-flowered Camellia species in China. According to Norihiko et al. (2001), the red blossoms of Benibana-cha also contain anthocyanins, albeit in somewhat different proportions. In spite of findings from a recent study that anthocyanins are not present in many white and yellow flowers of ornamental plants , practically all of the core metabolites needed for anthocyanin production are found in the white grape hyacinth petal (Lou et al., 2014). In a similar vein, the sepals of developing white tea flowers (C. sinensis) contain maximum levels of proanthocyanidins between the stages of complete separation and full bloom (Joshi and Gulati, 2011;Yu et al., 2020).
Anthocyanins are one of the most significant secondary metabolites that plants produce, and research has shown that they possess anti-cancer, anti-oxidation, and antiatherosclerosis characteristics (Escribano-Bailon and Santos-Buelga, 2012;Jiang et al., 2021). Anthocyanins are responsible for the orange, red, magenta, violet, and blue flower colorations. The biochemical route that leads to petal pigment accumulation has been thoroughly studied, and the genes that code for the necessary enzymes and transcriptional factors have been identified (Escribano-Bailon and Santos-Buelga, 2012;Feng et al., 2013). Extensive research on the molecular mechanisms involved in the anthocyanin biosynthesis pathway (ABP) has been conducted. However, the majority of previous studies focused mostly on the coloration of fruits and leaves, while anthocyanin production in colored-petals is yet to be studied in C. reticulata. Anthocyanin biosynthesis and accumulation is a highly complicated process that is controlled by a number of enzymes, transcription factors, and influenced by varied external environmental factors such as light, water stress, and temperature (He and Giusti, 2010;Hu et al., 2015;Gao et al., 2021). As a result, additional research is required to get complete insights into the mechanism that underlies petal coloration in C. reticulata.
Early biosynthetic genes (EBGs) and late biosynthetic genes (LBGs) are the two categories of ABP genes that are found in dicots, along with phenylalanine ammonia-lyase (PAL) genes (Medda et al., 2021). The genes, chalcone synthase (CHS), chalcone isomerase (CHI), and flavanone 3hydroxylase (F3H) that are part of the common flavonoid biosynthesis pathway are included in the EBGs. These genes have an effect on downstream flavonoids; while LBGs consist of F3′H, flavonoid 3′, 5′ hydroxylase (F3′5′H), dihydroflavonol 4-reductase (DFR), and anthocyanidin synthase (ANS) as well as UDP-glucose:flavonoid 3-O-glucosyltransferase (UFGTs) (Ren et al., 2021). Flowers come in a variety of colors, and these colors are determined by a number of factors, including co-pigments, cell structure, pH, and other factors (Stavenga et al., 2021). There is a correlation between the MYB-bHLH-WD40 (MBW) complex and the transcriptional regulation of ABP. MYB (Allan et al., 2008;Fraser et al., 2013;Naik et al., 2021) and bHLH (Li et al., 2020a) are two of the components of this complex that act as positive regulators of ABP genes. In most cases, the expressions of these two components are unique to pigmented tissues. On the other hand, WD40 is found to perform analogous functions in pigmented and nonpigmented tissues alike. Despite this, WD40 is able to maintain the stability of the MBW complex (Hichri et al., 2011;Li et al., 2020b).
Frontiers in Genetics frontiersin.org Transcriptome and metabolome high-throughput sequencing technology have been widely used in recent years to explore diverse plants' color mechanisms. We performed transcriptome and metabolome profiling, annotation, and bioinformatics analysis to understand the molecular mechanism of C. reticulata petal color formation and differentiation. The results of this study will aid in the breeding of cultivars with specific-flower colors.
2 Materials and methods

Plant materials
Materials tested were C. reticulata flowers harvested from three different genotypes that had distinct blossom colors: A popular C. reticulata stock (red flower, 'HHYC') and its two cultivars C. reticulata 'Xuejiao' (pink flower, XJ), and C. reticulata 'Tongzimian' (white flower, TZM). Plant materials were collected in Zixi Mountain of Chuxiong (H: 2133m, E: 101°24′1″, N: 25°0′7″), Yunnan, China. They were planted at the experimental station of the College of Landscape and Horticulture, Southwest Forestry University (H:1950m, E: 102°45′53″, N:25°4′0″), Kunming, Yunnan, China. Flowers were sampled during blooming in March 2021. There was a total of 9 C. reticulata flowers (petaLs) used in this experiment. As a result, all of the data were acquired based on the three distinct biological replicates that were gathered for each sample. All of the fresh C. reticulata petals were flash frozen in liquid nitrogen, and then placed in an ice chest and stored in −80°C refrigerator until further use.

RNA sequencing
The methods for isolating RNA, purifying, and monitoring it, as well as building a cDNA library and sequencing, were adapted from Sterling et al. (2015). The Qubit ® RNA Assay Kit in Qubit ® 2.0 Flurometer (Life Technologies, Carlsbad, CA, United States), the RNA Nano 6000 Assay Kit of the Agilent Bioanalyzer 2100 system (Agilent Technologies, Santa Clara, CA, United States), and the NanoPhotometer ® spectrophotometer (IMPLEN, Westlake Village, CA, United States) were used to check, measure, and evaluate the purity, concentration, and integrity. The sequencing libraries were prepared with the NEBNext ® UltraTM RNA Library Prep Kit for Illumina ® (NEB, Ipswich, MA, United States) in accordance with the recommendations by the manufacturer, and they were sequenced using an Illumina Hiseq platform. The NEBNext ® UltraTMRNA Library Prep Kit for Illumina ® (NEB, Ipswich, MA, United States) was used to prepare sequencing libraries, which were then sequenced using an Illumina Hiseq platform. This was done in accordance with the instructions provided by the manufacturer.

Identification of differentially expressed genes and transcription factors
With the help of the iTAK software (Zheng et al., 2016), transcription factors were identified. Previous methods by Zhao et al. (2012); Fraser et al. (2013) and Li et al. (2018) were adapted for both the identification and classification of transcription factors. Using a model of negative binomial distribution as a foundation, the DESeq R package (1.10.1) (Anders, 2010) was utilized to perform differential expression analysis. The false discovery rate was purposefully kept under control by the altered p values, which were changed using the methodology developed by Benjamini and Hochberg (1995). Genes were considered to have differential expression if their adjusted p-value was less than 0.05.

Real-time quantitative PCR
The RevertAidTM First Strand cDNA synthesis kit was utilized to analyze approximately 1 μg of total RNA in preparation for cDNA synthesis (Thermo Fisher Scientific, Waltham, MA, United States). An amplified reaction was carried out in a volume of L using the LightCycler ® 480 realtime PCR equipment in conjunction with a 96-well plate. Each reaction mix contained 1.0 μl of cDNAs, SYBR Premix ExTaq ™ 10 μl, PCR forward primer (10 μmol·L-1) 0.5 μl, PCR reverse primer (10 μmol·L − 1 ) 0.5 μl, ddH 2 O 8.0 μL, to a final volume of 20 μL. The amplified reaction began with 95°C, maintained for 5 min, and was then followed by 45 cycles of 10 s at 95°C, 20 s at 60°C, and 20 s at 72°C. After each experiment was finished, a melt-curve analysis was performed with the default parameters (5 s at 95°C and 1 min at 65°C). Normalization was accomplished with β-actin gene as control (Selvey et al., 2001;Tang et al., 2022). Every analysis was carried out in triplicate with separate biological replicates serving as controls. Table 1contains a listing of the primer sequences.

Extraction of metabolites
Metabolites were extracted using a 50% methanol buffer after the samples were thawed on ice (50% solution of methanol in distilled water). Following incubation at room temperature for 10 min and vortexing for 1 min to extract the 20 μl of sample in 120 μl of precooled 50% methanol, the extraction liquid was stored overnight at 60°C. In 96-well plates, the supernatants were transferred following a 20-min centrifugation at 4000 rpm. The extraction and subsequent analysis of data was conducted in three separate instances. An ultra-performance liquid chromatography (UPLC) system was used for all chromatographic separations (SCIEX, Cheshire, United Kingdom). The reversed phase separation was performed using an ACQUITY UPLC BEH Amide Frontiers in Genetics frontiersin.org column (100 mm × 2.1 mm, 1.7 m, Waters, United Kingdom). The temperature of the column oven was kept at 35°C. Solution A (25 mM ammonium acetate +25 mM NH 4 H 2 O) and Solution B were used as the mobile phase. The flow rate was 0.4 ml/min. The following were the gradient elution conditions: 95% B for the first 0.5 min; 95%-65% B for the next 0.5-9.5 min; 65%-40% B for the next 9.5-10.5 min; 40% B for the next 10.5-12 min; 40% B for the next 12-12.2 min; and 95% B for the next 12.2-15 min. Each sample received a 4 L injection. The metabolites eluted from the column were detected using a highresolution tandem mass spectrometer TripleTOF5600plus (SCIEX, Cheshire, United Kingdom). Positive and negative ion modes of the Q-TOF were both utilized. The chromatograph and mass spectrometer were controlled by XCMS software 3.2.0 (UC, Berkeley, CA, United States).

Detection and evaluation of metabolite concentrations
The MSConvert software was utilized to convert raw LC-MS data into the mzXML format, which was further processed by the XCMS, CAMERA, and metaX toolboxes, which were all integrated within the R program (Li et al., 2018;Stanstrup et al., 2019). To identify each ion, the combined retention time (RT) and m/z data was utilized. The Plant Metabolic Pathway Databases (PLANTCYC; http://www.plantcyc.org/), Kyoto Encyclopedia of Genes and Genomes (KEGG; http:// www.kegg.jp/), and Saskatoon Public Library in-House

Databases
(http://spldatabase.saskatoonlibrary.ca/), and Human Metabolome Databases (HMDB; http://www.hmdb.ca/ ) were used to perform level-one and level-two identification and annotation to explain the physico-chemical properties as well as the biological functions of metabolites. The metaX software (http://metax.genomics.cn/), was employed for screening and quantitative analysis of differential metabolites.

Statistical analysis
The relative gene expressions were computed utilizing the 2 −ΔΔCT method (Livak and Schmittgen, 2001), and data visualization was performed utilizing GraphPad Prism 5 (GraphPad Software Inc., San Diego, California, United States). Both the heatmap and the cluster analyses were performed with the help of R-3.4.2 and MEGA6, respectively. For the analysis of significant differences, IBM SPSS (IBM SPSS Software Inc., San Diego, California, United States) was used.

Flower color variation among Camellia genotypes
In attempt to understand the molecular mechanism underlying the color formation and differentiation in C.
3.2 Transcriptome profiling of petals from three contrasting flowers of C. reticulata The transcriptome profiling was performed on nine libraries (3 genotypes) with contrasting petals/flowers × 3 biological repeats) leading to an average of 47, 49 and 48 Gb data for HHYC, XJ and TZM, respectively (Table 2). Of these, an average of 45,500,491 bp were clean reads with an error rate of 0.02% (Table 2). Furthermore, the Q30 (%) and guanine-cytosine content ranged from 94.55 to 95.15 and 44.51-45.30%, respectively (Table 2). These results suggest that the transcriptome results qualified for further downstream analyses.
In all, 61,277 genes comprising 44,751 known genes and 16,526 novel genes were detected from the nine libraries (Supplementary Table S1). From these, 53,861, 53,739 and 54,631 genes were expressed in HHYC, XJ and TZM, respectively. The fragments per kilobase of exon per million fragments mapped (FPKM) distribution and principal component analysis (PCA) of the expressed genes are shown in Figures 2A,B. The Pearson correlation coefficients (r ≥ 0.81) and PCA based on FPKM showed high levels of reproducibility of petal samples from the same genotype (Supplementary Figure  S1). Most of the expressed genes were successfully annotated to at least one of the nine public databases, i.e., Kyoto Encyclopedia of Genes and Genomes (KEGG), KEGG pathway, NCBI nonredundant (NR), Swiss protein (Swissprot), KEGG orthologous (KOG), Gene ontology (GO), Protein family (Pfam) and Transcription factor (TF) family (Supplementary Table S1). The successful annotation of most of the identified genes coupled with high reproducibility and repeatability of petal samples from the same genotype give credence to our transcriptome results for the identification of differentially expressed genes (DEGs) and subsequent analyses.

Differentially expressed genes and kyoto encyclopedia of genes and genomes pathway analyses
A total of 11,746 DEGs were detected among the HHYC, XJ and TZM genotypes in pairwise comparisons upon application of strict criteria of log2 fold change (log2FC) ≥ 1 and false discovery rate (FDR) with adjusted p-value (padj) < 0.05. We subsequently performed hierarchical clustering heatmap of biological replicates from the three genotypes based on the FPKM values of the DEGs, the nine samples were divided into two major groups with each biological repeat in the sub-group. The two major groups comprised HHYC in group I, while XJ and TZM formed the group II ( Figure 3A). The grouping pattern suggests that XJ and TZM have some level of similarity in terms of flower coloration mechanism (Figure 1). Among the pairwise groups, 6169, 7512 and 5134 DEGs were detected in HHYC_vs_XJ, HHYC_vs_TZM and XJ_vs_TZM, respectively (Supplementary Tables S2A-C; Frontiers in Genetics frontiersin.org Figure 3B) with nearly equal proportion of up-and downregulated genes. The least number of DEGs in XJ_vs_TZM further pinpoints the similarity in flower coloration mechanism. We also compared the results of DEGs in each pairwise group and a total of 449 core conserved DEGs were detected among the three pairwise groups ( Figure 3C). These results highlight that XJ and TZM have similar flower coloration mechanism than HHYC.  In order to verify the transcriptome results, 18 genes were randomly selected for qRT-PCR to establish the correlation between gene expression profiles from the RNA-seq and qRT-PCR. The relative expression levels of the selected genes were mostly consistent with that of RNA-seq with correlation score of 84% (Supplementary Figures S2A,B). This further confirms the reliability of transcript quantification from our RNA-seq.

Analysis of key transcription factor families
It is well established that MYB-bHLH-WD40 (MBW) complex is involved in modulating coloration in plants (Fraser et al., 2013;Naik et al., 2021). However, in this study, no gene with WD40 was detected. Hence, we mined for MYB and bHLH transcription factors (TF) families among the DEGs and their log2 transformed FPKM were heatmapped. A total of 276 and 174 genes with MYB and bHLH TFs were identified, respectively, from which 42 MYB and 29 bHLH were detected among the DEGs. We further studied the expression profiles of these genes (Supplementary Figure S3A,B). For instance, the expression of gene-LOC114321352, gene-LOC114268549, gene-LOC114264400 and gene-LOC114281024 followed in HHYC or TZM were greater than XJ (Supplementary Figure S3A). In similar vein with rearrangement of the three genotypes by the 29 DEG encoding bHLH, gene-LOC114233324, gene-LOC114258501, gene-LOC114298318, gene-LOC114287545 and gene-LOC114319413 expressed lower in XJ relatively to either HHYC or TZM (Supplementary Figure S3A). The identified MYB-bHLH genes could be useful as genomic resources to deepen our understanding of transcriptional regulation flower coloration in C. reticulata.
3.6 Analysis of carotenoid biosynthetic pathway genes among three contrasting C. reticulata genotypes Thirty-eight carotenoid biosynthetic structural genes were differentially expressed among the three pairwise groups. Their expression profiles were subjected to hierarchical clustering heatmap. The expression profiles of the 38 genes grouped the 3 genotypes of C. reticulata into 2 groups (Figure 4), just as the Frontiers in Genetics frontiersin.org hierarchical clustering heatmap of the 11,746 DEGs in Figure 3A. It is reported that carotenoid gives yellow-orange-red color to fruits/flowers in plants (Mortensen, 2006;Tanaka et al., 2008). Two ortholog genes (gene-LOC114265684 and novel.6342) of At4g38540 encoding for zeaxanthin epoxidase (ZEP) on the carotenoid biosynthetic pathway expressed more than 10-fold higher in HHYC than either XJ or TZM (Figure 4), suggesting these genes could be candidate genes for red pigmentation in HHYC. Several studies have demonstrated that ZEP is the major contributor to carotenoid composition (Gonzalez-Jorge et al., 2016;Lee et al., 2021). From the Figure 4, another gene worth elaborating is gene-LOC114269817 encoding for abscisate betaglucosyltransferase (ABA-UGT) [EC:2.4.1.263]. It was nearly not expressed in either XJ or TZM. This suggests that gene-LOC114269817 may contribute to red coloration in HHYC as a recent study by Li et al. (2022) revealed that increase in ABA-UGT gene contributed to red coloration in strawberry fruit, but its reduction contributed to white coloration in an unripe fruit.
3.7 Analysis of phenylpropanoid/ flavonoid/anthocyanin biosynthetic pathways genes among the three contrasting C. reticulata genotypes Aside carotenoid biosynthesis, it is well known that flavonoids are major molecules involved in plant pigmentation specifically anthocyanins (Lai et al., 2014). Phenylpropanoids are precursors of flavonoid/flavone and flavonol which lead to anthocyanins biosynthesis and accumulation. Therefore, we explored DEGs in the pairwise groups based on the KEGG pathway enrichment analyses (Supplementary Table S3). Sixteen structural genes involved in flavone and flavonol biosynthesis (Ko0094) were differentially expressed between HHYC and XJ (Supplementary Table S2A). Out of these 16 structural genes, gene-LOC114309343 (flavonoid 3′-monooxygenase [EC:1.14.14.82]) and other 2 genes (gene-LOC114273964 and gene-LOC114273955) encoding flavonol-3-O-glucoside/galactoside glucosyltransferase and anthocyanidin 3-O-glucoside 2″-O-glucosyltransferase-like, were not expressed in XJ genotype, making them potential candidate genes for the white pigmentation in XJ.

Flavonoid-based metabolome profiling of flowers of the three contrasting C. reticulata
To further understand the metabolites involved in the reddish color formation in HHYC compared to the two other Heatmap clustering of thirty-eight structural genes involved in carotenoid biosynthesis among the three contrasting flowers of C. reticulata (HHYC), C. reticulata 'Xuejiao' (XJ) and C. reticulata 'Tongzimian' (TZM). The expression levels of the genes were log2 transformed from averages of the three biological repeats from each genotype and the transformed data were used to draw the heatmap. Color diagram at the right side shows the intensity of the gene expression.  Table S4). We performed PCA and hierarchical heatmap clustering based on the ion intensities of the 310 compounds. We observed that the grouping/topology ( Figures 5A,B) followed the FKPM of genes obtained from the RNA-seq ( Figure 3A). To identify the differentially accumulated metabolites (DAMs), we performed partial least squares discriminant analysis (PLS-DA) with threshold of log2FC ≥1 and variable importance in projection (VIP) ≥1. From these stringent criteria, we detected 138, 149 and 133 DAMs in HHYC_vs_XJ, HHYC_vs_TZM and XJ_vs_TZM, respectively (Supplementary Figures S5A,B; Supplementary Table S5).
From above, 18 anthocyanin-based DAMs were detected with an accumulation pattern: HHYC > XJ > TZM (Table 3). All the 18 DAMs were detected in HHYC, while only nine and seven DAMs were detected in XJ and TZM genotypes, respectively. These trends suggest that anthocyanin is the major pigment responsible for the flower coloration among the three contrasting genotypes of C. reticulata. peonidin-3-O-glucoside (pmf0203) accumulated higher in HHYC genotype but were either completely absent or lower in either XJ or TZM (Table 3).

FIGURE 5
Profiles of metabolites accumulated among the 3 from petals of contrasting flowers of C. reticulata (HHYC; red colored), C. reticulata 'Xuejiao' (XJ; pink colored) and C. reticulata 'Tongzimian' (TZM; white colored). (A). Principal component analysis based on metabolites normalized ion intensities. Mix samples represent equal proportion of samples from HHYC, XJ and TZM. Metabolome profiling was done in triplicates for each genotype (HHYC-1 to -3, XJ-1 to -3, TZM-1 to -3 and mix-1 to -3). (B). Hierarchical clustering heatmap based on metabolites normalized ion intensities. The red, blue and green colors on column cluster (shown righthand side) represent HHYC, XJ and TZM, respectively. The ion intensity of each metabolite are shown green color represents lowly-accumulated, and red color represents highly-accumulated metabolites.
Frontiers in Genetics frontiersin.org

Conjoint analysis of transcriptome and metabolome results
We further performed a conjoint analysis of the transcriptome and metabolome data to assess the relationship between the genes and metabolites. This helps to statistically prioritize genes and metabolites that alter flower colors in C. reticulata. It was observed that anthocyanin is the key pigmentation that differentiates flowers of HHYC_vs_XJ and HHYC_vs_TZM (p = 0.004), however XJ_vs_TZM were not differentiated, possibly as a result of high similarity between the flowers of XJ and TZM as observed in the heatmap clustering ( Figure 3A: Figure 5B; Figure 6).

Discussion
Flower color is one of the most important features of horticultural crops, including C. reticulata (Xin et al., 2015). It does not only affect economic value, but eye appeal, aesthetics and quality. With more than 500 cultivars of C. reticulata with diverse flower color phenotypes (Xin et al., 2015), the regulatory mechanism underlying flower coloration among the cultivars remain unknown.
Flavonoid/anthocyanins, betalains, carotenoids and chlorophylls are well known to modulate color pigmentation in plants (Mortensen, 2006;Tanaka et al., 2008;Lai et al., 2014;Zhao et al., 2017;Shen et al., 2019). With exception of chlorophylls, the other three pigmentations have been elucidated to regulate color variation depending on their synthesis, structures and subcellular localization (Zhao and Tao, 2015). Carotenoids are well documented to confer Frontiers in Genetics frontiersin.org yellow, orange and red colorations , while anthocyanins are responsible for blue, purple, red or white colorations (Yan and Kerr, 2013). The transcriptome profile from the flowers with red, white and intermediate (pink) genotypes (HHYC, XJ and TZM, respectively) revealed that carotenoids and anthocyanins (Supplementary Table S3) coexist which were pronounced in the red flower genotype (HHYC) (Figure 1). For instance, two ZEP genes (gene-LOC114265684 and novel.6342) and one ABA-UGT (gene-LOC114269817) (Figure 4) involved in carotenoid biosynthesis as well as flavonoids/anthocyanins structural genes (Supplementary Table S2) including three anthocyanidin 3-Oglucosyltransferase (gene-LOC114289234, gene-LOC114285765, gene-LOC114285774 and gene-LOC114288492) expressed higher in the HHYC genotype than XJ and TZM genotypes. These results highlight the co-regulation of red flower color in HHYC genotype by anthocyanins and carotenoids, with a least contribution of the latter. Our findings are consistent with the findings of Li et al. (2022); Xue et al. (2016). In addition, downregulation/absence of anthocyanidin 3-O-glucosyltransferase is responsible for inhibition of the synthesis of delphinidin 3glucoside (Takeda et al., 1990;Yabuya et al., 2002;Muhammad et al., 2022) and Yamazaki et al. (1999) andNakatsuka et al. (2008) reported that a reduced level of delphinidin limits red coloration. The absence/downregulation of anthocyanidin 3-O-glucosyltransferase genes in XJ or TZM implies that late stage of anthocyanins biosynthesis may have been blocked (Lai et al., 2013). The two ZEP genes (gene-LOC114265684 and novel.6342), one ABA-UGT (gene-LOC114269817) and three anthocyanidin 3-O-glucosyltransferase (gene-LOC114289234, gene-LOC114285765, gene-LOC114285774 and gene-LOC114288492) could be targeted to study the structural variation which is responsible for color variation between XJ and TZM. There are reports that nucleotide variation in structural genes involved in anthocyanin biosynthesis could cause variations in coloration in some plants (Ben-Simhon et al., 2015).
Transcription factors are proteins involved in converting, or transcribing DNA into RNA. It is well established the involvement of MYB-bHLH-WD40 (MBW) complex in modulating coloration (B). Heatmap of seven differentially accumulated metabolites. Metabolites with hollow circles were not detected in this study, while those in purple circles were detected in this study. Heatmaps were produced with log2 transformed data. Color diagram at the right side shows the intensity of the metabolite accumulation or gene expression.
Frontiers in Genetics frontiersin.org in plants (Li et al., 2020a;Naik et al., 2021). Among the 61,277 expressed genes including 44,751 known genes and 16,526 novel genes were detected in the present study, no WD40 structural gene was found (Supplementary Table S1). Four MYB transcription factor genes (gene-LOC114321352, gene-LOC114268549, gene-LOC114264400 and gene-LOC114281024) and four bHLH genes (gene-LOC114233324, gene-LOC114258501, gene-LOC114298318, gene-LOC114287545 and gene-LOC114319413) expressed higher in genotypes with red color (either HHYC or TZM) than in the white pigmented genotype, XJ (Supplementary Figure S3). These suggest that MYB and bHLH transcription factors are potential positive regulators of carotenoid and anthocyanin biosynthesis in C. reticulata. Zhou et al. (2019) revealed that CpbHLH1/2 modulates carotenoid biosynthesis-related genes in papaya, while He et al. (2022) demonstrated that UpMYB44 controls carotenoid biosynthesis-related genes in Ulva prolifera.
In comparison with other crops, MYB-bHLH role in regulating carotenoid and anthocyanin biosynthesis needs to be unearthed with genes reported in this study (Ramsay and Glover, 2005;Gonzalez et al., 2008;Zhao et al., 2012). It is important to link the identified key TFs to the target structural genes involved in the synthesis of anthocyanins through gene co-expression analysis. Such analysis will require more transcriptome data, which could be the subject of a further study. We further performed targeted metabolome profiling with the petals of three contrasting genotypes of C. reticulata with the aim of identifying specific flavonoid/anthocyanin compound(s) that may account for white coloration observed in XJ and TZM genotypes relative to HHYC genotype. Targeted metabolome profiling offers absolute quantification, higher accuracy, and greater selectivity (Beckles and Roessner, 2012;Cao et al., 2020;Lelli et al., 2021). From this, 310 metabolites comprising 282 and 28 flavonoid/ anthocyanin and tannin compounds, were detected, respectively (Supplementary Table S4). These compounds mostly accumulated higher in either HHYC or TZM genotype than in XJ genotype (Supplementary Table S4). For example, Luteolin-7-O-gentiobioside which is flavone compound accumulated in order of HHYC (61,965,667) > TZM (1,657,300), but was completely absent in XJ. This compound is reported to cause pink coloration in Myoporum bontioides (Myoporaceae) (Iwashina and Kokubugata, Frontiers in Genetics frontiersin.org 2014; Cao et al., 2020). In all, eighteen anthocyanin compounds were detected among the three contrasting genotypes with wide variation in flower colors (Table 3). Of these, only 9 and 7 compounds were detected in XJ and TZM genotypes, respectively. The absence/ presence of some anthocyanin compounds could be the basis for color differentiation among the three genotypes of C. reticulata. Cyanidin, a major anthocyanin compound which accounts for red coloration had 8 derivatives compounds (cyanidin-3-O-(6″-Omalonyl) glucoside, cyanidin-3-O-(6″-O-p-coumaroyl-2″-Oxylosyl) glucoside, cyanidin-3-O-arabinoside, cyanidin-3-Ogalactoside, cyanidin-3-O-glucoside (kuromanin), and cyanidin-3-O-rutinoside (keracyanin). Of these, cyanidin-3-O-glucoside, and cyanidin-3-O-rutinoside were exclusively detected in HHYC genotype which is consistent with the findings of Kim and Kokubugata (2010) who found these compounds in red pigmented flower cultivar of buckwheat (Gan-Chao), but were absent in the white-flowered cultivar (Tanno). In another study, Li et al. (2020b) found delphinidin 3-O-rutinoside and petunidin 3-O-glucoside to be absent in white fruit of Lycium ruthenicum Murray. Similar observations were made among the white, violet and red flowers of Rhododendron schlippenbachii Maxim (Kim et al., 2013). All these indicate that the absence and reduction of anthocyanin compounds could account for the white and pink (red-white) flower coloration in XJ and TZM, respectively.
In summary, this study provides important insights into flower coloration in C. reticulata by unearthing key structural genes and transcription factors. This forms a foundation for biotechnological applications to explore their regulatory mechanisms in flower coloration (Kim et al., 2013;Varshney et al., 2020Varshney et al., , 2021. Also, the flavonoid/anthocyanin compounds identified in this study could be validated and used to develop biomarkers for metabolome-assisted breeding programs (Varshney et al., 2020;Jamaloddin et al., 2021).

Data availability statement
The data presented in this study are deposited in the National Genomic Data Center under the accession number PRJCA012977 (https://ngdc.cncb.ac.cn/search/?dbId=&q=PRJCA012977).