Integrated analysis of the metabolome and transcriptome in different developmental stages of Cerasus humilis peel coloration

Background: Cerasus humilis is a unique dwarf shrub and fruit color is an important trait in the species. In this study, we evaluated the transcriptomic and metabolomic proles of the plant at different developmental stages to elucidate the mechanism underlying color formation. Results: In a metabolomics analysis, 16 anthocyanin components were identied at different developmental stages, and high levels of cyanidin O-syringic acid and pelargonidin 3-O-beta-d-glucoside (callistephin chloride) were correlated with the reddening of the fruit peel. Additionally, transcriptome analysis showed that most anthocyanin biosynthetic genes and two MYB transcription factors were signicantly up-regulated. qRT-PCR results for these differentially expressed genes were generally consistent with the high-throughput sequencing. A co-expression analysis revealed that ANS and UFGT play keys role in pigmentation. Conclusions: These ndings provide insight into the mechanisms underlying anthocyanin accumulation and coloration during fruit peel development, providing a basis for the breeding of anthocyanin-rich C. humilis cultivars.

anthocyanins. Furthermore, MYB transcription factors (TFs) are key factors in regulatory networks controlling growth and development, metabolism, and responses to biotic and abiotic stresses [23]. For example, anthocyanin accumulation is regulated by the MYB transcription factor [7,[24][25][26]. Anthocyanin biosynthesis, controlled by a diverse array of exogenous and endogenous factors, is the most wellstudied secondary metabolic pathway in plants. Gene expression controlling the accumulation of metabolites under different stages of development is a direct determinant of fruit color changes.
However, the molecular and metabolic pathways involved in C. humilis fruit coloration are still elusive. Therefore, the genetic basis of anthocyanin biosynthesis and related metabolic pathways requires further exploration. Knowledge of the genetic basis of characters related to fruit coloration will facilitate the manipulation of this trait to obtain more attractive and healthier fruits for consumers [27].
Transcriptomics and metabolomics can provide insight into transcriptional regulation and metabolic ow.
Bai et al. [28] performed a spatiotemporal transcriptomic analysis of the red Chinese sand pear peels, revealing that a light-responsive pathway functions in anthocyanin accumulation with different temporal expression patterns, they suggested that PpMYB10 as well as other light-responsive transcription factors are involved in the regulation of anthocyanin biosynthesis. Using a combination of transcriptomic and metabolite pro ling, Zhuang et al. [29] found that purple turnip essentially diverts dihydro avonols to the biosynthesis of anthocyanins over avonols by strongly down-regulating one avonol synthase gene, while strikingly up-regulating dihydro avonol 4-reductase, anthocyanidin synthase, and UDP-glucose: avonoid-3-O-glucosyltransferase compared to levels in green turnip. In the same way, Meng et al. [30] demonstrated the molecular mechanism underlying the WP1-mediated regulation of oral carotenoid pigmentation by network analyses of metabolomics and transcriptomics data. Thus, multi-omics joint analysis is a powerful means to elucidate the complex process of plant development [7,31].
In the present study, based on UPLC-MS/MS, metabolite changes in the process of peel color change in C. humilis were qualitatively and quantitatively analyzed, and the unique pattern of anthocyanin-related metabolites during peel color change was established for the rst time. Additionally, a correlation analysis of the transcriptomic pro les and anthocyanin contents was performed to elucidate the mechanism underlying variation in fruit peel pigmentation in C. humilis, followed by further analyses of the regulatory mechanism based on gene expression assays. Our transcriptome analysis and metabolic pro ling of four stages of fruit development, with a focus on anthocyanin biosynthetic pathways, will facilitate breeding and the functional characterization of genes of interest in C. humilis as well as several other species in Cerasus Mill.

Results
Principal component analysis (PCA) reveals differences in anthocyanin metabolite pro les Traditionally and practically, the fruiting of C. humilis is classi ed into four developmental stages. Metabolite content data processing was performed to improve normalization using R (www.r-project.org/) to compare the accumulation of metabolites in different samples by HCA. Based on relative differences in accumulation patterns among different stages, the HCA of anthocyanins resulted in four main clusters in C. humilis (Fig. 1a). Anthocyanins belonging to clusters 1 and 2 accumulated at the highest levels in S4, while anthocyanins within clusters 3 and 4 did not differ signi cantly among stages.

Overview of transcriptome sequencing
To assess the gene expression pro le of the fruit peel at S1-S4, RNAs from 12 C. humilis fruit samples were used for RNA-seq with three replicates per stage. Before using Illumina platform to read the cDNA libraries, the low-quality reads and adapter sequences were deleted. The clean reads for each sample ranged from 40,704,832 to 47,000,394, and a total of 541,781,582 clean reads were obtained. The Q30 values > 92%, meeting the criterion for gene discovery. Based on ltered clean data, the full-length transcript sequence was assembled using Trinity, and the total number of assembled transcripts was 176,109, with a mean length of 1340.05 nt and N50 length of 2,227 nt. These were assembled into 73,035 unigenes with a mean length of 1,018 nt and N50 length of 1,917 nt (Additional le 2: Table S2).

Annotation and identi cation of unigenes
All assembled unigenes and ORFs were annotated using Trinotate against the public databases NT, NR, UniProt, RNAMMER, eggNOG, GO, and KEGG. Using the Trinotate annotation results, genes were subjected to a GO analysis. According to the GO database, 73,035 unigenes are divided into 60 functional groups,, in the three main functional categories (biological process, cellular component, and molecular function) (Additional le 6: Figure S1). In the biological process category, 42,236 (57.83%) unigenes were associated with metabolic processes.
The assembled unigenes were searched through the COG (Cluster of Orthologous Group, KOG) protein database to evaluate potential functions. By a KOG analysis, 25,127 unigenes were clustered into 25 functional categories (Additional le 7: Figure S2). Among them, general function prediction only (Class R, 5,338 unigenes) was most common, accounting for 21.24% of unigenes. Class R is a category for which only a general prediction regarding protein function is feasible (e.g., biochemical activity). Most of the COGs could be assigned to a well-de ned functional category. However, the fact that the largest category remains unspeci c in nature indicates insu cient information within the COG database. The next largest categories were signal transduction mechanisms (Class T, 3,071 unigenes, 12.22%), posttranslational modi cation, protein turnover, chaperones (Class O, 2,593 unigenes, 10.32%), and translation, ribosomal structure and biogenesis (Class J, 1,111 unigenes, 7.02%). The smallest group was cell motility (Class N, 39 unigenes, 0.02%). The COG and GO annotations showed that unigenes expressed in the C. humilis fruit peel encode various proteins related to metabolism.
To gain more insight into the transcriptomic differences, we compared DEGs in C. humilis fruit peel among different developmental stages. Using the DESeq2 package, normalized RPKM (Reads Per Kilobase Million Mapped Reads) values were obtained from RNA-seq data to quantify transcript expression. We identi ed 7,930 DEGs in comparisons among different fruit developmental stages.
Among the three comparisons (S2_S1, S3_S2, and S4_S3), the largest number of DEGs was found between S3 and S2, with 4,333 DEGs, including 1,540 up-regulated and 2,793 down-regulated DEGs. In addition, there were 430 up-regulated DEGs in S2_S1, and 1,179 up-regulated DEGs in S4_S3 (Fig. 2). These results indicated that DEGs were most active during the second developmental stage of the fruit peel and were perhaps largely responsible for anthocyanin biosynthesis. Therefore, we speculate that the transition from S2 to S3 is a key stage in fruit peel development.

GO enrichment and KEGG pathway analyses of DEGs
We functionally categorized the DEGs based on GO terms. A total of 7,927 unigenes were annotated, including 429 up-regulated and 201 down-regulated unigenes in the S2_S1-enriched GO term libraries and 1,539 up-regulated and 2,794 down-regulated DEGs in the S3_S2-enriched GO term libraries. The subcategories with the highest degree of enrichment were all cell part, followed by binding and cellular process, as shown in a histogram in Additional le 8: Figure S3. Apart from this, 1,178 up-regulated and 1,786 down-regulated unigenes in the S4_S3 comparison were enriched for GO terms, particularly cell part, followed by cellular process, and binding.
To classify the DEGs based on related pathways, a KEGG enrichment analysis was performed. In total, 125 biosynthetic and metabolic pathways were enriched in the pairwise comparisons of fruit peel developmental stages. The avonoid biosynthetic pathway began to be enriched (ko00941, 10 genes) in S3_S2, and the anthocyanin biosynthetic pathway began to be enriched (ko00942, 2 genes) in S4_S3. We obtained 12 DEGs across the three comparisons related to pigment synthesis, suggesting that these genes are associated with fruit peel coloration (Additional le 9: Figure S4).

Page 6/23
The changes of TFs in the process of peel coloring Next, we focused on differentially expressed TFs based on our transcriptome data. The MYB and bHLH TF families were the most abundant TF families during the fruit coloring process, followed by the B3, AP2/ERF, and NAC families ( Table 1), many of these were expressed before samples began to show color.
Additionally, the TFs were mostly up-regulated before the color change and down-regulated after the color change. Notably, the MYB and bHLH TF families are important transcription factors for anthocyanin transcriptional regulation, and 68% of bHLHs and 70% of MYB genes were down-regulated. Interestingly, we found that 2 MYB genes (TRINITY_DN26515_c0_g1 and TRINITY_DN21536_c0_g1) were signi cantly up-regulated after the fruit peel color changes, and their expression levels were consistent with the pattern of fruit peel pigment accumulation, indicating that they might be key regulators of fruit peel color.

Analysis of unigenes related to anthocyanidin biosynthetic pathways in the fruit peel
Previous studies have suggested that the anthocyanin biosynthetic pathway is an important metabolic branch of the avonoid pathway, responsible for the production of anthocyanins in different plant tissues. Thus, we evaluated the mechanisms underlying red pigmentation in fruit peel based on a comparative transcriptome analysis. The compositions of compounds in anthocyanin biosynthesis pathways differed signi cantly depending on the growth stage (Fig. 3). The PAL gene regulates the synthesis of cinnamic acid and consistently showed a gradual decrease over the four growth stages. Seven 4CL genes exhibited signi cant differences in expression levels, including 2 up-regulated and 5 down-regulated genes. The CHS gene is the rst key enzyme in the avonoid pathway, which regulates the synthesis of chalcone; it consistently showed a gradual increase over the four growth stages. The expression levels of CHI, F3H, DFR, and ANS were signi cantly up-regulated from S1 to S4, and the expression levels were positively correlated with the accumulation of anthocyanin. UDP-glucose: avonoid-3-O-glycosyltranferase (UFGT), which mainly catalyzes the transformation of unstable anthocyanin glycosides into stable anthocyanin. By using a transcriptomic analysis, two UFGT genes (TRINITY_DN23815_c1_g1 and TRINITY_DN19893_c1_g5) involved in the anthocyanin biosynthesis pathway were screened (Fig. 3). In order to test the accuracy of the transcriptome data, we selected 12 anthocyanin synthesis genes and 2 MYB genes for qRT-PCR veri cation. The results showed that the expression trend of qRT-PCR was consistent with the transcriptome data, indicating that the transcriptome data has high reliability (Fig. 4).

Integrated analysis of the transcriptome and metabolome
By transcriptome and metabolome analyses, we constructed a co-expression network to identify interactions (Fig. 5). We found that 11 anthocyanin biosynthesis DEGs and 11 differential metabolites exhibited interactions. Highly positive correlations were obtained among most differential metabolites, and only pmb0542 (cyanidin 3-O-malonylhexoside) was minimally relevant to all other differential metabolites. The correlation analysis between DEGs and differential metabolites also showed that levels of 5 DEGs (TRINITY_DN15939_c0_g1, TRINITY_DN20784_c0_g2, TRINITY_DN21962_c0_g2, TRINITY_DN23815_c1_g1, and TRINITY_DN26276_c0_g1) were positively correlated with all differential metabolites, while the other 6 DEGs (TRINITY_DN17945_c0_g2, TRINITY_DN21141_c0_g1, TRINITY_DN25788_c2_g8, TRINITY_DN22885_c0_g4, TRINITY_DN25099_c1_g5, and TRINITY_DN24108_c0_g2) were negatively correlated with all differential metabolites (Fig. 5a, Additional le 4: Table S4). In addition, the combined metabolome and transcriptome analysis revealed a signi cant positive correlation between both ANS (TRINITY_DN15939_c0_g1) and UFGT (TRINITY_DN23815_c1_g1) (PCC > 0.82) and most metabolites (Fig. 5b). These data indicated that expression patterns of ANS and UFGT have decisive roles in color development, and their regulatory relationships need to be further studied.

Discussion
Fruit peel color is an important trait that effects the edible value of many fruits. Several recent reports have described key genes and enzymes for fruit peel color formation and anthocyanin accumulation in fruit trees, including pear and apple [13,26]. However, the mechanism underlying C. humilis anthocyanin biosynthesis is not well understood. In addition, as far as we have known, the differences of metabolites of C. humilis peel color changes at different developmental stages have not been studied. We rst completed a transcriptome analysis and metabolic pro ling for four developmental stages of fruit development to focus on anthocyanin biosynthetic pathways.
Based on UPLC-MS/MS, changes in metabolites during peel color changes in C. humilis were qualitatively and quantitatively analyzed. The unique patterns of anthocyanin-related metabolite changes in the process of peel color change were preliminarily established for the rst time. Using a targeted metabolomics approach, we detected and annotated 11 metabolites with signi cant differences among developmental stages of fruit peel, all of which were up-regulated. Notably, the levels of two anthocyanins, cyanidin O-syringic acid and pelargonidin 3-O-beta-d-glucoside (callistephin chloride), increased signi cantly during fruit peel development in C. humilis. Zhou et al. [33] found that cyanidin 3-O-glucoside and pelargonidin 3-O-beta-d-glucoside are unique anthocyanins in pink tea owers. Accordingly, we suggested that the accumulation of these anthocyanins accounted for the red pigmentation of fruit peel during fruit development.
It has been deeply studied in antho biosynthesis patyway, and the structure genes who play important roles in this pathway has also been well know [33]. We found that most of upstream structural genes (i.e. PAL, C4H, and 4CL) in the early committed steps were highly active in young fruits but were gradually silenced during fruit development (Fig. 4). In contrast, UFGT genes (TRINITY_DN23815_c1_g1 and TRINITY_DN19893_c1_g5), which are involved in the synthesis and accumulation of anthocyanins, were gradually activated during ripening process of C. humilis. UFGT played an important role in anthocyanin metabolism, which catalyzes the glycosylation of anthocyanins, and makes anthocyanins exist in plant vacuoles in the form of stable anthocyanins. Our results suggested that the accumulation of anthocyanin metabolites is positively correlated with UFGT expression levels during the reddening of C. humilis fruit peel. The accumulated metabolites are converted by UFGT into the stable anthocyanins involved in fruit peel reddening with fruit ripening [31]. Similarly, Steyn et al. [34] reported that UFGT activity increases in pear fruit development. Therefore, we proposed that UFGT during fruit ripening could potentially promote C. humilis peel coloring.
Studies of C. humilis fruit color have revealed that in addition to DEGs, TFs play a key role in the regulation of color changes and the activity of pigment biosynthetic pathways [35]. MYB TFs contribute to the regulation of plant secondary metabolism, including anthocyanin biosynthesis [26,36]. For example, AcMYB110a regulates anthocyanin production in petals by promoting the expression of F3GT1 in Actinidia [37]. McMYB10 plays an important role in ever-red leaf coloration by positively regulating McF3'H in crabapple [38]. Consistent with these results, we identi ed several candidate MYB genes (TRINITY_DN26515_c0_g1 and TRINITY_DN21536_c0_g1) with the potential to positively regulate the majority of avonoid biosynthetic genes, their expression abundance was consistent with anthocyanin accumulation. Therefore, the results suggested that these 2 MYB genes might play an important role in the temporal and spatial regulation of anthocyanin synthesis.
Metabolites can provide direct insight into biological processes and mechanisms. In this study, we detected 16 metabolites, including 11 differential metabolites, providing reference values for the isolation and identi cation of functional compounds in fruit peel. Our combined metabolome and transcriptome analysis uncovered a signi cant positive correlation between the synthesis of two genes (ANS, TRINITY_DN15939_c0_g1; UFGT, TRINITY_DN23815_c1_g1) in the anthocyanin pathway (ko00942) and differential metabolites (Fig. 5b), establishing an new interesting connection between relative secondary metabolites and genes. The high expression of ANS and UFGT, as major determinants of anthocyanin biosynthesis in fruit peel, might account for C. humilis fruit coloration.

Conclusions
Collectively, the mechanisms underlying anthocyanin accumulation in C. humilis fruit peel were analyzed by MRM, transcriptomics, and qRT-PCR. At least 16 anthocyanins exhibited differences among developmental stages, and cyanidin O-syringic acid and pelargonidin 3-O-beta-d-glucoside (callistephin chloride) were identi ed as the major anthocyanins. FurthermoreMoreover, The anthocyanin biosynthesis genes were identi ed and analyzed by transcriptome data and qRT-PCR. These results indicated that ANS and UFGT were the key genes involving in anthocyanin accumulation. Meanwhile, two MYB genes (TRINITY_DN26515_c0_g1 and TRINITY_DN21536_c0_g1) participate anthocyanin synthesis pathway and potentially regulated the expression of related genes in C. humilis. This study improves our understanding of anthocyanin accumulation and the molecular mechanism underlying anthocyanin biosynthesis in C. humilis fruit peel. However, the precise mechanism requires further studies.

Metabolome analysis
The sample preparation, analysis, metabolite qualitative and quanti cation were carried out by MWDB (Metware Biotechnology Co., Ltd., Wuhan, China) according to their standard procedures [39][40][41][42]. The instrument system of metabolite data acquisition was mainly ultra-high performance liquid chromatography (HPLC). The data of metabolite content were normalized by range method and R software was used (www.r-project.org/), and the accumulation patterns of metabolites among different samples were analyzed by hierarchical cluster analysis (HCA). The analysis of metabonomics data combines the methods of univariate and multivariate statistical analysis, and according to the characteristics of data, multi-perspective analysis is also carried out, and nally the differential metabolites will be accurately mined.

RNA isolation and Illumina sequencing
Total RNA was extracted from four developmental stages of the fruit peel using the Omega RNA Extraction Kit (Biel, Switzerland). The 1% agarose gel electroassay was used to detect degradation and impurities. RNA quantity and quality were determined using a Keio K5500 spectrophotometer (Kaio, Beijing, China), Agilent Bioanalyzer 2100 system, and Agilent RNA 6000 Nano Kit (Agilent Technologies, Palo Alto, CA, USA). A total of 4 g of total RNA per sample (12 samples = 4 stages × 3 biological replicates) was used as input material for RNA-seq library preparation (Annoroad Biotechnology, Zhejiang, China). Trinity was for de novo assembly. Raw data (in FASTq format) were processed using Perl scripts to measure Q30 values and GC contents. After removing low-quality reads, the remaining clean data were used for all downstream analyses. Heat map and PCA of the relative differences in anthocyanins in different fruit development stages in Cerasus humilis. a: Heat map visualization. Each sample is visualized in a single column and each metabolite is represented by a single row. Red indicates high abundance and green indicates low abundance. b: Score plots for PCA. young fruit (S1), green fruit (S2), slightly red fruit (S3), and red fruit (S4)