Comprehensive Analysis of Transcriptome and Metabolome Reveals the Flavonoid Metabolic Pathway Is Associated with Fruit Peel Coloration of Melon

The peel color is an important external quality of melon fruit. To explore the mechanisms of melon peel color formation, we performed an integrated analysis of transcriptome and metabolome with three different fruit peel samples (grey-green ‘W’, dark-green ‘B’, and yellow ‘H’). A total of 40 differentially expressed flavonoids were identified. Integrated transcriptomic and metabolomic analyses revealed that flavonoid biosynthesis was associated with the fruit peel coloration of melon. Twelve differentially expressed genes regulated flavonoids synthesis. Among them, nine (two 4CL, F3H, three F3′H, IFS, FNS, and FLS) up-regulated genes were involved in the accumulation of flavones, flavanones, flavonols, and isoflavones, and three (2 ANS and UFGT) down-regulated genes were involved in the accumulation of anthocyanins. This study laid a foundation to understand the molecular mechanisms of melon peel coloration by exploring valuable genes and metabolites.


Introduction
Melon (Cucumis melo L.) belongs to the Cucurbitaceae family, and it is one of the most important economic crops worldwide [1]. Melon peel is essentially the boundary of the fruit. It not only maintains the integrity of the fruit but also protects it from the external environment. Peel color is one of the important characteristics that affect the commercial value of melon fruit and consumer choice [2]. Peel color variety of melon fruit includes white, yellow, orange, and green [3]. Melon peels are rich in nutritional ingredients, such as carbohydrates and minerals, and they contain significant amounts of dietary fibers and antioxidants, such as polyphenols and flavonoids [4].
Flavonoids are a large class of polyphenolic secondary metabolites. There are mainly six groups of flavonoids in plant tissues, including flavones, flavanones, isoflavones, anthocyanins, flavanols, and flavonols [5]. Most flavonoids are identified as antioxidants, protecting plants against various stresses [2,6]. In addition, flavonoids are the pigments that affect the coloration of plant organs [7]. In humans, flavonoids play important roles in maintaining normal vascular permeability and protecting against diseases, such as hyperglycemia, cancer, and diabetes [8].
In recent years, high-throughput methods have been widely applied to understand pigmentation in plants. By comprehensive analysis of the metabolome and transcriptome, the relationship between the content of various secondary metabolites and the corresponding differentially expressed genes has been revealed. Some structural genes and transcription factors were found to be involved in flavonoid biosynthesis. Three UDP-glucose flavonoid

Identification of DEGs Related to Peel Color
Three cDNA libraries were constructed by sequencing RNA extracted from the peels of W, B, and H samples ( Figure 1). The transcriptome sequencing obtained a total of 68.84 Gb clean data. More than 92.08% of the clean data had scores greater than Q30. Interestingly, 87.26% to 89.41% of clean reads were mapped to the reference genome (Table S1). The correlation coefficients of three biological replicates were more than 0.97 ( Figure 2A). A total of 22,090 unique genes were detected among the three peel samples. Using cutoff criteria with a |fold change| ≥ 2 and FDR < 0.05, DEGs from each two-treatment comparison (i.e.,  Figure 2C). Additionally, these genes differed between the green and yellow peels, implying that they may be responsible for the yellow peel formation.  Figure 2C). Additionally, these genes differed bet the green and yellow peels, implying that they may be responsible for the yellow formation.

Functional Analysis of DEGs
In order to identify the function of DEGs in the formation of fruit peel color three categories were classified, including molecular function, cellular component biological process, according to GO classifications. In total, 59 GO terms were signific enriched between W and B, most of which were enriched in the catalytic activity, bin cell part, metabolic process, and other functional categories ( Figure S1, Table S5). Bet W and H, there were 179 GO-annotated DEGs, which were mainly categorized int cell part, binding, single-organism process, metabolic process, and other functional gories ( Figure 3A, Table S6). One hundred ninety-two GO terms between B and H enriched in the binding, catalytic activity, cell, cell part, cellular process, and other tional categories ( Figure 3B, Table S7). Then, we used KEGG enrichment analysis o DEGs to confirm the DEG-associated pathways in the formation of the peel color. Bet the W and B libraries, KEGG analysis revealed that beta-alanine metabolism (ko00 insulin signaling pathway (ko04910), and carotenoid biosynthesis (ko00906) were th three significantly changed pathways. Between the W and H libraries, KEGG analys vealed the Fanconi anemia pathway (ko03460), steroid hormone biosynthesis (ko00 and retinol metabolism (ko00830) were the top three significantly changed pathway tween the B and H libraries, KEGG analysis revealed Fanconi anemia pathway (ko03   Figure 2C). Additionally, these genes differed between the green and yellow peels, implying that they may be responsible for the yellow peel formation.

Functional Analysis of DEGs
In order to identify the function of DEGs in the formation of fruit peel coloration, three categories were classified, including molecular function, cellular component, and biological process, according to GO classifications. In total, 59 GO terms were significantly enriched between W and B, most of which were enriched in the catalytic activity, binding, cell part, metabolic process, and other functional categories ( Figure S1, Table S5). Between W and H, there were 179 GO-annotated DEGs, which were mainly categorized into the cell part, binding, single-organism process, metabolic process, and other functional categories ( Figure 3A, Table S6). One hundred ninety-two GO terms between B and H were enriched in the binding, catalytic activity, cell, cell part, cellular process, and other functional categories ( Figure 3B, Table S7). Then, we used KEGG enrichment analysis of the DEGs to confirm the DEG-associated pathways in the formation of the peel color. Between the W and B libraries, KEGG analysis revealed that beta-alanine metabolism (ko00410), insulin signaling pathway (ko04910), and carotenoid biosynthesis (ko00906) were the top three significantly changed pathways. Between the W and H libraries, KEGG analysis revealed the Fanconi anemia pathway (ko03460), steroid hormone biosynthesis (ko00140), and retinol metabolism (ko00830) were the top three significantly changed pathways. Between the B and H libraries, KEGG analysis revealed Fanconi anemia pathway (ko03460),

Functional Analysis of DEGs
In order to identify the function of DEGs in the formation of fruit peel coloration, three categories were classified, including molecular function, cellular component, and biological process, according to GO classifications. In total, 59 GO terms were significantly enriched between W and B, most of which were enriched in the catalytic activity, binding, cell part, metabolic process, and other functional categories ( Figure S1, Table S5). Between W and H, there were 179 GO-annotated DEGs, which were mainly categorized into the cell part, binding, single-organism process, metabolic process, and other functional categories ( Figure 3A, Table S6). One hundred ninety-two GO terms between B and H were enriched in the binding, catalytic activity, cell, cell part, cellular process, and other functional categories ( Figure 3B, Table S7). Then, we used KEGG enrichment analysis of the DEGs to confirm the DEG-associated pathways in the formation of the peel color. Between the W and B libraries, KEGG analysis revealed that beta-alanine metabolism (ko00410), insulin signaling pathway (ko04910), and carotenoid biosynthesis (ko00906) were the top three significantly changed pathways. Between the W and H libraries, KEGG analysis revealed the Fanconi anemia pathway (ko03460), steroid hormone biosynthesis (ko00140), and retinol metabolism (ko00830) were the top three significantly changed pathways. Between the B and H libraries, KEGG analysis revealed Fanconi anemia pathway (ko03460), biosynthesis of secondary metabolites (ko01110), and protein processing in endoplasmic reticulum (ko04141) were the top three significantly changed pathways (Table 1) biosynthesis of secondary metabolites (ko01110), and protein processing in endoplasmic reticulum (ko04141) were the top three significantly changed pathways (Table 1). Clearly, isoflavonoid biosynthesis (ko00943), flavone, and flavonol biosynthesis (ko00944), and carotenoid biosynthesis (ko00906) were significantly changed in W vs. H and B vs. H (Table  1). These pathways may influence fruit peel coloration.

Identified Metabolites Involved in Flavonoid Biosynthesis
We profiled the metabolome analysis of the different peel samples via LC-MS. Abundance lower than 9.00 was not detected. The DEMs were set as fold change ≥ 2 (upregulation) or ≤ 0.5 (down-regulation). The results of PCA showed that the metabolites were clearly separated among the three samples ( Figure 4A (Table S8). Eighty-eight up-regulated and 25 down-regulated DEMs were identified between W vs. H ( Figure 4C), including 53 flavones, 17 flavanones, 17 flavonols, 7 anthocyanins, 7 isoflavones, and other compounds (Table S9). In total, 101 up-regulated and 28 down-regulated DEMs were identified between B vs. H ( Figure 4D (Table S10) (Table S10). Interestingly, we observed a clear distinction between the yellow peel and green peel samples. Moreover, compared to W vs. B, there were more of the same metabolites between W vs. H and B vs. H. In particular, metabolites of flavanones, flavonols, and isoflavones and anthocyanins were present.

KEGG Enrichment Analysis of Significantly Differential Metabolites
The DEMs were mapped to KEGG metabolic pathways. KEGG enrichment analysis showed that the DEMs were mainly involved in four pathways: flavonoid biosynthesis (Ko00941), flavone and flavonol biosynthesis (ko00944), isoflavonoid biosynthesis (Ko00943), and anthocyanin biosynthesis (Ko00942). The four pathways all existed in W vs. H and B vs. H ( Figure 5B,C). Additionally, only the flavone and flavonol biosynthesis was enriched in the group W vs. B ( Figure 5A). Therefore, the differences between the green peel and yellow peel were mainly caused by the biosynthesis of flavonoid, isoflavonoid, and anthocyanin. The color and size of the dots represent p-value and the amount of enriched differential metabolites, respectively. Rich factor means the ratio of the number of differential metabolites to the total number of metabolites enriched in a specific category.  H. KEGG analysis. The x-axis represents the richness factor. The color and size of the dots represent p-value and the amount of enriched differential metabolites, respectively. Rich factor means the ratio of the number of differential metabolites to the total number of metabolites enriched in a specific category.
A total of 40 kinds of flavonoids were enriched in four KEGG pathways, including 11 flavones, 10 flavanones, 9 flavonols, 7 isoflavones, and 3 anthocyanins (Table 2). Among these flavonoids, 36 up-accumulated and 4 down-accumulated in the H (Tables S9 and  S10). Most compounds (flavones, flavanones, flavonols, and isoflavones) were significantly higher in the H sample than in the W and B samples, but the content of anthocyanins was lower in the H sample. Furthermore, the levels of flavonoids accumulation were similar in W and B, but their accumulation pattern was more complicated in H. We analyzed the flavones in the melon peel. Flavones were detected to be the maximum number of metabolites among metabolites with significant content differences among three melon peel samples. Eleven flavones were identified. All these flavones were significantly higher in W vs. H and B vs. H, ranging from 3.51-fold to 20.18-fold increments. The content of these flavones has no other significant difference in the W and B samples, except luteolin 7-O-glucoside and chrysoeriol. The fold change of luteolin 7-O-glucoside was significantly decreased (−16.67-fold) in W vs. B (Table 2).
We analyzed the flavanones in the melon peel. Ten kinds of significantly different flavanones were detected. The content of flavanones was considerably higher in H compared with W and B. Moreover, it showed a 10. 36 Table 2).
We analyzed the flavonols in the melon peel. Nine flavonols in the melon peel were identified. The content of flavonols was higher in the H sample compared with W and B samples.

Modulation of Flavonoid Biosynthesis Genes and Metabolites in Relation to Fruit Peel Coloration
In order to further analyze the relationship between genes and metabolites of flavonoids involved in the coloration of yellow fruit peel, all results of metabolites and genes were combined to establish a network ( Figure 6). The results showed that 12 DEGs, 10 DEMs, and their derivatives participated in flavonoid biosynthesis pathway in W, B, and H. The 12 genes with significant differences of flavonoid biosynthesis were as follows: two 4CL (MELO3C017009.

Discussion
Melon shows a large variation in rind color, such as white, yellow, orange, and green [3]. Throughout the development of the fruit, the complicated network of metabolites and genes is dramatically transformed [19]. Fruit peel secondary metabolites, such as pigments, aroma compounds, and tannins, affect fruit appearance [20]. In the current study, we analyzed different metabolites, and flavones, flavanones, flavonols, isoflavones, and anthocyanins were commonly responsible for rind color differences. By analyzing tran- Using Pearson's correlation coefficient [7], we performed comprehensive analyses of transcriptome and metabolome data to investigate the association between genes and metabolites involved in flavonoid biosynthetic pathway. In total, we detected 145 significant correlations (correlation coefficient, R 2 > 0.9) between six DEGs (MELO3C035535.2, MELO3C017219.2, MELO3C005571.2, MELO3C009387.2, MELO3C014584.2, MELO3C035771.2) and 26 DEMs, including five flavones, seven flavanones, six flavonols, six isoflavones, and two anthocyanins (Table S11). Each metabolite was correlated with many different genes. Interestingly, flavonol and isoflavone shared the largest number of common genes. This suggested that flavonol and isoflavone might have evolved similar accumulation mechanisms. Of these genes, the correlations between MELO3C035535.2 and 25 metabolites were greater than 0.97, and MELO3C035771.2 had the highest correlation score (1) with homoeriodictyol.

Discussion
Melon shows a large variation in rind color, such as white, yellow, orange, and green [3]. Throughout the development of the fruit, the complicated network of metabolites and genes is dramatically transformed [19]. Fruit peel secondary metabolites, such as pigments, aroma compounds, and tannins, affect fruit appearance [20]. In the current study, we analyzed different metabolites, and flavones, flavanones, flavonols, isoflavones, and anthocyanins were commonly responsible for rind color differences. By analyzing transcriptome data, we uncovered several DEGs related to flavone and flavonol biosynthesis, isoflavonoid biosynthesis, and anthocyanin biosynthesis that are probably involved in rind color development. In addition, a combined transcriptome and metabolome study was conducted to elucidate the mechanism of fruit rind coloration.
Plant metabolomic analysis enables us to study the relationship between metabolites produced by biological processes and plant characteristics [21]. Flavonoids are vital pigments for regulating the rind color of many fruits [11,19,22]. There are six main types of flavonoids in plant tissues: flavones, isoflavones, flavonols, flavanones, flavanols, and anthocyanins. Furthermore, in yellow melons, naringenin chalcone is considered the major flavonoid pigment in rinds [3]. In our study, flavones, flavonols, flavanones, and isoflavones were up-regulated, and anthocyanins were down-regulated in the H samples, compared with W and B samples. Our results perfectly corroborate the conclusions of Wang et al. [23], who proposed that the levels of flavones, flavonols, isoflavones, and anthocyanins determined the green and yellow color of citrus peel. However, they did not analyze the genes associated with peel pigmentation. We used a combination of transcriptome and metabolome data to uncover genes involved in the biosynthesis of flavones, flavonols, isoflavonoids, and anthocyanins, thus searching for valuable information to understand the changes in melon rind color.
Genes that participate in flavonoid biosynthesis and regulation have been reported in melon [15,24]. However, many genes have been discovered and identified via traditional research technology. Transcriptome analysis is considered an important method for studying the expression level, structure, and function of genes in order to uncover phenotypic traits. Thus, the combination of metabolome and transcriptome analyses has increasingly become a practical tool for the mining of genes involved in various metabolic pathways [25].
Using phenylalanine as precursor, flavonoids are synthesized via flavonoid metabolic pathway. First, phenylalanine catabolism is catalyzed by phenylalanine ammonia-lyase (PAL), cinnamic acid 4-hydroxylase (C4H), and 4 coumarate CoA ligase (4CL) to form the starting substrate for flavonoid synthesis, p-coumaroyl-CoA. Then, the chalcone precursors for all flavonoids are the condensation of p-coumaroyl-COA and three molecules of malonyl-COA by chalcone synthase (CHS) [15,26]. Subsequently, chalcone is converted to naringenin by chalcone isomerase (CHI). The following reaction catalyzed by flavanone 3-hydroxylase (F3H) produces dihydrokaempferol, which is further transformed into dihydroquercetin and dihydromyricetin under the actions of flavonoid 3 -hydroxylase (F3 H) and flavonoid 3 ,5 -hydroxylase (F3 5 H) [27]. Finally, consecutive reactions catalyzed by dihydroflavonol reductase (DFR) and anthocyanidin synthase (ANS), which can be further methylated, glycosylated, and acylated, lead to the formation of a wide variety of species-specific anthocyanins [28]. Flavone synthase (FNS) and flavonol synthase (FLS) are responsible for the formation of flavones and flavonols, respectively [29]. Previous studies identified some key structural genes from the flavonoid pathways that regulate coloration. Structural gene 4CL was down-regulated, which largely explained the reduced accumulation of naringenin chalcone and naringenin [17]. Here, we found that MELO3C017009.2 and MELO3C035535.2 may also be involved in the high accumulation of the naringenin chalcone and naringenin content in melon. When the expression of F3H gene decreased, the dihydrokaempferol component did not accumulate and the color of petals and fruits had been changed [30,31]. Additionally, the high expression of F3 H provides sufficient content of dihydroquercetin [32]. We further identified that the genes F3H (MELO3C035771.  (Tables S3  and S4). Moreover, they resulted in dihydrokaempferol and dihydroquercetin accumulation. Overexpression of IFS resulted in the accumulation of isoflavones formononetin, daidzein, and genistein glycosides in transgenic plants [27,33,34]. High expression of FNS mostly accumulates apigenin and luteolin. The high expression levels of FLS observed here are consistent with the essential role in the biosynthesis of the flavonols quercetin and kaempferol, as well as with the critical role of CnFLS1 in yellow pigmentation in C. nitidissima [35]. In agreement with our findings, the genes IFS (MELO3C010951.2), FNS (MELO3C005570.2), and FLS (MELO3C035771.2) were proposed as important structural genes for regulating the pigmentation in melon. Moreover, the genes ANS and UFGT were found to differentially regulate the anthocyanin accumulation in different melon cultivars [9,26], which was also revealed in this study. Obviously, 12 DEGs associated with isoflavonoid biosynthesis (Ko00943), flavone, and flavonol biosynthesis (Ko00944), and biosynthesis of secondary metabolites (Ko01110) was also related to the formation of peel color. The difference in the color of fruit is mainly due to the differential accumulation of flavonoids and anthocyanins [22]. The present study confirmed these conclusions, and the transcriptome and metabolome analyses highlighted major changes in flavonoid biosynthetic pathway-related genes and metabolites, which may correspond with the changes in peel coloration of melon.

Plant Materials
Three different melon cultivars were used in this study. The fruits were classified into three types according to the peel color: grey-green (W), dark-green (B), and yellow (H) (Figure 1). The fresh fruits were harvested at the maturity stage (35 days after anthesis for W and B, 45 days after anthesis for H) at the Institute of Vegetables and Flowers of the Chinese Academy of Agricultural Sciences in Beijing. Fruits were grown in greenhouse conditions during the spring season of 2018. The fruits were selected for uniformity in color, shape, and size, and they were quickly transported to the laboratory [36]. Fruit peels (0.1 cm thick) were carefully excised with scalpels, collected, immediately frozen in liquid nitrogen, and then stored at −80 • C for subsequent transcriptome sequencing and metabolite extraction. Three independent biological replicates were collected per treatment.

Transcriptome Sequencing
Nine libraries representing the three peel samples and the three replicates were constructed for RNA-Seq. Total RNAs were extracted from W, B, and H samples using an RNA extraction kit (DP441, TIANGEN, Beijing, China). The quantity and quality of mRNAs were measured using a Nano-Drop ultraviolet spectrophotometer (Nanodrop Technologies, Thermo, Waltham, MA, USA) and a Bioanalyzer 2100 System (Agilent Technologies, Santa Clara, CA, USA), respectively. RNA integrity of the samples was determined by 1% agarose gel electrophoresis.
The sequencing libraries were constructed using 3 µg RNA per sample and sequenced using 150 bp paired-end Illumina Hiseq 4000 platform. Briefly, the mRNA was purified using magnetic beads with Oligo (dT). We synthesized first-strand cDNA with random hexamer primer and M-MuLV reverse transcriptase and then second-strand cDNA was synthesized using RNase H and DNA polymerase I. DNA fragments were first repaired at the end. A tail was added, the sequencing connector was connected, and the fragments were purified using the AMPure XP beads to select cDNA fragments of 240 bp in length. PCR was performed, and a cDNA library was constructed by purifying PCR products with the AMPure XP beads. The libraries were sequenced on an Illumina HiSeq platform, and 150 bp paired-end reads were produced.

Transcriptome Data Analysis
Low-quality sequencing reads containing adapter and poly-N were removed before downstream analyses. The filtered reads were mapped to the reference genome sequence using HISAT2 software [37]. The gene expression level was calculated by the FPKM (fragments per kilobase of transcript per million fragments mapped) method. Differential expression analyses among the three groups (W vs. B, W vs. H, B vs. H, with three biological replicates per treatment) of colored samples were conducted using the DESeq R package (1.10.1) [38], with adjusted p-values. Genes with the following parameters were considered differentially expressed genes (DEGs): |Log 2 fold change| ≥ 1 and false discovery rate (FDR) <0.05. Gene Ontology (GO) annotation and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analyses of DEGs were applied, employing clusterProfiler software [39].

Metabolite Qualification and Quantification
To monitor the repeatability of the analytical process, a quality control sample consisting of a mixture sample extract of every three analytical samples was used during instrumental analysis. Then, raw mass spectrometry (MS) data were imported into Analyst 1.6.3. Metabolites were annotated by public databases, including KNAPSAcK, MassBank, MoToDB, METLIN, and HMDB. The metabolites were quantified using MRM. We used MultiaQuant for peak areas, peak integration, and peak alignment calculations. The differentially expressed metabolites (DEMs) among three peel samples were evaluated using univariate analysis methods, such as fold change (FC) analysis, and multivariate analysis methods, such as partial least squares-discrimination analysis (PLS-DA). The different metabolites were identified with the following parameters: metabolites with VIP ≥1 and fold change ≥2 or fold change ≤0.5 [40]. The different metabolites were subjected to data normalization.

Integrated Metabolome and Transcriptome Analysis
We used principal component analysis (PCA) to compare the trends of metabolites. Enrichment analysis of functions and signaling pathways of the differentially expressed genes was conducted by KEGG database. A p-value less than 0.05 was defined as DEGs. Pearson's correlation coefficients were calculated between the genes and metabolites. On the basis of calculation results, DEGs and DEMs with a coefficient of R 2 > 0.8 were selected. Then, gene and metabolite datasets were analyzed by canonical correlation analysis (CCA) to reveal melon fruit peel color metabolites and gene molecular interactions, as well as to construct a network.

Conclusions
In the present study, the regulatory network of fruit peel coloration in melon was carried out by transcriptome and metabolome for the first time. We identified differentially expressed genes associated with fruit rind coloration between yellow and green peel color, among which 4CL, F3H, F3 H, IFS, FNS, FLS, ANS, and UFGT play important roles in flavonoid biosynthetic pathways. Metabolomic analysis indicated that crucial flavone, flavonol, flavanone, isoflavone, and anthocyanins responsible for fruit peel pigmentation showed significantly different expression between yellow and green rind types. In summary, the flavonoids and their correlated genes detected in this study provide a vital metabolic and functional role for flavonoid biosynthetic pathways in melon peel coloration. The findings from our study will benefit molecular breeders in improving the quality of melon fruit.
Supplementary Materials: The following are available online, Table S1: Statistics of the reads' quality in the RNA-Seq study; Table S2: FPKMs and function annotation of genes in W vs. B.; Table S3: FPKMs and function annotation of genes in W vs. H.; Table S4: FPKMs and function annotation of genes in B vs. H.; Table S5: GO annotation of differentially expressed genes in W vs. B.; Table S6: GO annotation of differentially expressed genes in W vs. H.; Table S7: GO annotation of differentially expressed genes in B vs. H.; Table S8: Differentially expressed metabolites in W vs. B.; Table S9: Differentially expressed metabolites in W vs. H.; Table S10: Differentially expressed metabolites in B vs. H.; Table S11: Pearson correlation (r 2 ) between metabolites and genes. Figure S1: GO analysis of DEGs in W vs. B.
Author Contributions: H.W. and Q.F. designed studies and contributed to the original concept of the project. A.Z., J.Z., X.C., X.S., and Q.F. managed the project. A.Z., J.Z., and Q.F. prepared the figures, tables and drafted the manuscript. All authors have read and agreed to the published version of the manuscript.

Data Availability Statement:
The data presented in this study are available on request from the corresponding author.