Integrative Metabolome and Transcriptome Analyses Reveal the Pericarp Coloration Mechanisms in Bitter Melon (Momordica charantia L.)

Abstract

Pericarp colors are critical agronomic traits that affect the quality and economic values of fruits. Although a diversity of bitter melon pericarp (BMP) colors is available, the fruit pigmentation mechanisms remain elusive. Hence, this study aimed to unveil the key metabolites and molecular mechanisms underlying variation in BMP coloration through integrative metabolomics and transcriptomics analyses of four differently colored genotypes, including K1102 (grayish orange), 262 (grayish yellow), 1392 (very soft green), and K115 (dark grayish cyan). The four BMPs exhibited significant metabolite profile and transcriptional differences, as over 112 and 1865 DAMs (differentially accumulated metabolites) and DEGs (differentially expressed genes), respectively, were identified. The variation in the content of six anthocyanins, including malvidin 3-O-glucoside, petunidin 3-O-glucoside, rosinidin O-hexoside, cyanidin, cyanidin 3-p-hydroxybenzoylsophoroside-5-glucoside, and pelargonidin 3-O-beta-D-glucoside, might be the major driving factor of BMP color changes. Notably, malvidin 3-O-glucoside, rosinidin O-hexoside, and petunidin 3-O-glucoside are the dominant pigments in K115, while carotenoids and other flavonoids may contribute to other colors. Candidate flavonoid structural and regulatory (MYBs, NACs, MADSs, bHLHs, and bZIPs) genes were identified. Of them, gene13201 (anthocyanin reductase), gene8173 (polyphenol oxidase), gene2136 (NAC43), gene19593 (NAC104), and gene15171 (tetrapyrrole-binding protein) might play essential roles in K115 pericarp color development. Our findings deepen our understanding of BMP pigmentation and provide fundamental resources for higher-valued bitter melon breeding perspectives.

1. Introduction

Bitter melon or bitter gourd (Momordica charantia L.) is a highly economic value and promising industrial vegetable crop of the Cucurbitaceae family [1]. It is cultivated in tropical regions, mainly in Asian and African countries, where its fruits are consumed as an essential vegetable for many culinary uses and commonly exploited for medicinal purposes [1,2]. For instance, pharmacological investigations have proved that bitter melon possesses diverse health-promoting abilities, such as antidiabetic, anti-cancer, anti-obesity, anti-neurodegenerative, antioxidant, and anti-dyslipidemia properties [2,3,4,5,6,7,8]. It has also shown potential for treating SARS-CoV-2 [9]. Its fruits are mainly consumed immature due to their typical bitter flavor induced by the high accumulation of cucurbitacins (cucuriate triterpenoids) [10]. Preferences for bitter melon fruits are attributable to many traits, among which fruit color and fruit ridges are the most critical [10]. As per cucumber, bitter melon fruit displays a variety of colors ranging from white to dark green. QTLs and candidate genes underlying the variation in bitter melon fruit ridges and yield-related traits have been explored [10,11]. However, the molecular determinant of the fruit pericarp pigmentation remains elusive. Fruit color determines essential quality traits and has been the interest of tremendous studies in many horticultural crops, such as apple [12], cucumber [13,14], mango [15], raspberry [16], papaya [17], chokecherry [18], banana [19], etc. These studies have provided vital molecular resources for breeding preferential fruits in these crops. Hence, the molecular dissection of BMP pigmentation mechanisms will generate key targetable metabolic resources for improving the quality and commercial values of bitter gourd fruits.

Fruit color is primarily determined by anthocyanins, a subclass of flavonoid compounds widely distributed in higher plants [20,21]. The most commonly found anthocyanidins are cyanidin, pelargonidin, delphinidin, petunidin, peonidin, and malvidin, and these compounds display different colors varying from red and blue to purple depending on the pH, light, and temperature in plant cells [20]. For example, glycosylated cyanidins are the dominant pigments in apple, berry, and dark red Chinese cherry [12,20,22]. Petunidins and malvidins are the major pigments of black currants and matured darker red wine, respectively [20]. In black goji, petunidin was also identified as the major pigment [23]. Apart from the primary anthocyanin pigments, other phenolic compounds, including carotenoids, chlorophyll, and other flavonoids (flavone, flavanone, flavanols, etc.), influence the final color appearance of fruits [14,24,25]. In cucumber, studies have shown that variations in the composition and amounts of carotenoids and chlorophyll contribute to the yellow, orange, red, light green, green, dark green, and white peel and flesh colors [13,14,26]. We thus hypothesized that color formation in bitter melon fruit might be a complex mechanism involving several metabolites and genes. Although flavonoid profiling of different colored BMPs was performed, no anthocyanin was detected, limiting the exploration of the pigmentation metabolism [27].

Anthocyanins are biosynthesized via the flavonoid pathway from phenylalanine through subsequent catalyzing action of multiple structural and regulatory genes. The main enzyme-encoding structural genes include PAL (phenylalanine ammonia-lyase), 4CL (4-coumarate-CoA ligase), C4H (cinnamate 4-hydroxylase), CHS (chalcone synthase), CHI (chalcone isomerase), F3′H (flavonoid 3′-hydroxylase), F3H (flavanone 3-hydroxylase), F3′5′H (flavonoid 3′,5′-hydroxylase), DFR (dihydroflavonol-4-reductase), ANS/LDOX (anthocyanidin synthase/leucoanthocyanidin dioxygenase), UFGT (UDP-glucose: flavonoid-3-Oglucosyl-transferase), and LAR (leucoanthocyanidin reductase), ANR (anthocyanidin reductase), and GST (glutathione S-transferase) [20,24,28]. The expressions of these genes are regulated by a complex of TFs (transcription factors) at the transcriptional levels. This complex is composed of R2R3-MYBs, bHLHs (basic helix–loop–helix), and WD40 proteins (MBW complex) [20,24,28]. Of them, MYBs play critical roles in regulating the composition and concentration of anthocyanins in plant organs [29,30,31]. In addition to these TFs, many other TFs, such as NAC, bZIP, WRKY, SPL, GATA, etc., network with the MBW complex to modulate the biosynthesis and storage of anthocyanins [24,28]. The anthocyanin component and regulation mechanisms vary depending on the species. Therefore, it is of significant interest to explore the potential candidate genes and reveal the molecular mechanism governing variation in anthocyanin levels of different colored bitter melon fruits. The integrative analysis of metabolome and transcriptome is the most common approach used to effectively break down the metabolic pathways and regulatory genes of diverse fruits’ pigmentation [18,22,32,33,34].

In the present study, we combined metabolomics and transcriptomics analyses of four differently colored bitter melon genotypes, including K1102 (grayish orange), 262 (grayish yellow), 1392 (very soft green), and K115 (dark grayish cyan), to unveil the key anthocyanins and potential candidate genes underlying variation in BMP colors. In addition, we investigated the pericarp pigmentation-associated molecular mechanisms and the variability of metabolites in the different colored BMPs. Our results deepen our understanding of bitter melon fruit color variation and lay a molecular foundation for suitable color breeding in the future.

2. Materials and Methods

2.1. Plant Material and Experimental Procedures

Four bitter melon varieties with different fruit colors, including K1102 (grayish orange), 262 (grayish yellow), 1392 (very soft green), and K115 (dark grayish cyan) (Figure 1A), were used in the present study. They were planted in the Modern Agricultural Science and Technology Innovation Demonstration Park of Sichuan Academy of Agricultural Sciences (30.77° N, 104.21° E) on 11 March 2022. The photon flux density ranged from 600–900 W m². The temperature and relative humidity during the plants’ growth varied: 15–30 °C and 50–85%, respectively. The soil in this region is predominantly loam and slightly acidic with a pH of 6.6. Eighteen days after simultaneous pollination, fifteen fruits with a length of 10–30 cm and 4.2–6.7 cm in diameter were selected for each genotype. Then, the flesh in the middle of the fruit was immediately collected, mixed, and divided into two parts for transcriptome sequencing and metabolomics profiling. All samples were collected in triplicate. The collected samples were in situ, frozen in liquid nitrogen, and stored at −80 °C up until the analyses.

2.2. LC–MS-Based Widely Targeted Metabolomics Analysis

Freeze-dried pericarp samples’ powder (100 mg) was extracted at 4 °C (overnight) with 1.2 mL of 70% methanol. After centrifugation (15 min at 12,000× g), the supernatants were collected and filtrated with a 0.22 μm micropore membrane (SCAA-104, ANPEL, Shanghai, China). All extracts were stored at −20 °C up until the UPLC–ESI–QqQLIT–MS/MS analysis at MWDB (Metware Biotechnology Co., Ltd., Wuhan, China) [35,36,37]. All sample extracts were mixed equally to constitute quality control (QC) samples. The metabolomics was achieved as per previously described methods, and detailed information about the liquid phase and MS conditions is presented in Table S5 [36,37].

The spectrum information, retention times, and mass spectra were integrated for the qualitative identification of the metabolites. Specifically, the values of Q1 (precursor ions) and Q2 (product ion), retention times, and fragmentation patterns were compared with standards when available (Sigma-Aldrich, St. Louis, MO, USA). When no standards were available, the metabolites were identified via the database of MWDB and verified in public databases (MassBank, KNApSAcK, MoTo DB, HMDB, and METLIN) [35,36,38]. The relative contents of metabolites were computed via the MRM modes (triple quadrupole (QqQ) MS analysis).

The multivariate analyses were achieved after data quality validation and standardization. The PCA (principal component analysis), HCA (hierarchical clustering analysis), and OPLS-DA (orthogonal partial least squares discriminant analysis) analyses were carried out in R with the packages prcomp, pheatmap, and MetaboAnalystR (www.r-project.org), respectively. Significant DAMs (differentially accumulated metabolites) were screened out considering fold change (FC ≥ 1), VIP (variable important in projection) ≥ 1, and p-value < 0.05. Finally, DAMs functional annotation was conducted by mapping to the KEGG (Kyoto Encyclopedia of Genes and Genomes) database and subsequent significant enrichment analyses by MSEA (metabolite sets enrichment analysis).

2.3. RNA Extraction, Library Construction, and Sequencing

Total RNA from each flesh sample was extracted with a Trizol reagent kit (Invitrogen, Carlsbad, CA, USA) per the manufacturer’s instructions. The genomic DNA was discarded using DNase I (TaKaRa, Beijing, China). RNA quality was investigated on an Agilent 2100 Bioanalyzer (Agilent Technologies, Palo Alto, CA, USA) and quantified using the ND-2000 (NanoDrop Technologies). Only high-quality RNA (OD260/280 = 1.8~2.2, OD260/230 ≥ 2.0, RIN ≥ 6.5, 28S:18S ≥ 1.0, >1 μg) samples were used for sequencing library construction. The libraries were prepared with 1 μg of total RNA using the TruSeqTM Kit (Illumina, San Diego, CA). Next, fragmentation was followed by cDNA construction using the NEB Next Ultra RNA Library Prep Kit (NEB #7530, New England Biolabs, Ipswich, MA, USA). After adapters ligation, the resulting cDNA libraries were subjected to Illumina sequencing (HiSeq xten/NovaSeq6000 sequencer) at Gene Denovo Biotechnology Co. Ltd. (Guangzhou, China). The raw paired-end reads’ quality was assessed with SeqPrep (https://github.com/jstjohn/SeqPrep, accessed on 11 September 2022) and Sickle (https://github.com/najoshi/sickle, accessed on 11 September 2022). Thereafter, clean reads were aligned to the M. charantia L. reference genome (https://www.ncbi.nlm.nih.gov/genome/?term=bitter+Melon, accessed on 11 September 2022) by HISAT2 (http://ccb.jhu.edu/software/hisat2/index.shtml, accessed on 11 September 2022) software [39]. Finally, the mapped reads of each sample were assembled using StringTie (http://www.string-db.org/, accessed on 11 September 2022) [40].

2.4. Differentially Expressed Genes (DEGs) and Functional Enrichment Analysis

All transcript expression levels were computed via the TPM (transcripts per million reads) method. Next, RSEM (http://deweylab.biostat.wisc.edu/rsem/, accessed on 11 September 2022) was used for gene abundance quantification [41]. DEGs analysis was conducted with DESeq2 software (v 1.30.0) [42]. Significant DEGs were detected at FDR (false discovery rate) < 0.05 and |fold change| ≥ 1. GO (Genes Ontology, http://geneontology.org/, accessed on 11 September 2022) and KEGG (Kyoto Encyclopedia of Genes and Genomes, http://www.genome.jp/kegg/kaas, accessed on 11 September 2022) analyses were achieved with GO seq and KOBAS software (v2.0), respectively. Significant enrichment terms were identified at p-value < 0.05.

2.5. Quantitative RT-PCR Analysis

RNA was extracted from leaf samples (RNA Extraction Kit, Code No. 9769S, Takara, Beijing, China). Reverse transcription (RT) was conducted using the PrimeScriptTM 2II1st Strand cDNA Synthesis Kit (Code No. 6210A, Takara, Beijing, China). Quantitative RT-PCR was achieved as per the guidance in TB Green^® Premix Ex TaqTM II (Code No. RR820A, Takara, Beijing, China). All samples had three biological and technical replicates. The PP2A (Serine/Threonine-protein phosphatase) gene was used as an internal control [43] to normalize the expression levels of target genes via the 2^−ΔΔCT method [44]. Primers were designed with primer5 software (Table S4).

2.6. Data Analysis

Excel 2016 software and GraphPad Prism (v9.0.0121, La Jolla, CA, USA) were used for data processing and graph construction. TBtools software was used for heatmaps’ construction [45].

3. Results

3.1. Metabolite Profiles of White, Grey White, Light Green, and Dark Green Bitter Melon Pericarps

To identify key metabolites associated with bitter melon pericarp (BMP) color variation, we conducted a widely targeted metabolomics analysis of K1102 (grayish orange), 262 (grayish yellow), 1392 (very soft green), and K115 (dark grayish cyan) pericarps (Figure 1A). We accurately identified a total of 573 metabolites (Table S1). The dominant metabolites were flavonoids (16.58%), organic acids and derivatives (15.71%), amino acids and derivatives (12.04%), lipids (11.34%), nucleotides and derivatives (8.20%), phenylpropanoids (6.98%), and alkaloids (5.93%) (Figure 1B). Of the flavonoids, flavones and flavonols were major, and only seven anthocyanins were identified (Figure 1C).

To explore metabolites’ differentiation and variability between the BMP of different colors, we carried out PCA, HCA, and correlations analyses (Figure 1D,E, and Figure S1). The PCA results showed that the metabolite profiles of different colored BMPs were very different as they clustered apart and could be discriminated by PC1 (29.10%) and PC2 (23.53%) (Figure 1D). All the QC samples were gathered together next to the PCA plot center, indicating the repeatability and reliability of the metabolomics analysis. As shown in Figure 1E and Figure S1, the results from HCA and correlation analyses confirmed the dissimilarity of the metabolite profiles of the K1102, 262, 1392, and K115 BMPs. Many metabolites exhibited higher relative content in 1392 (very soft green) compared to other colored pericarps (Figure 1E). For further confirmation of the observed differences in the metabolite profiles, we conducted OPLS-DA analyses. We obtained supportive results with high predictability and strong goodness of fit (R²Y = 1, Q² ˃ 0.973) (Figure S2).

3.2. Variation in Important Metabolic Classes and Differentially Accumulated Metabolites (DAMs)

To explore the distribution characteristics of important metabolites in the different colored BMPs, we calculated the sum of the relative contents of metabolites within the same class and compared their levels in the four sample groups (Figure 2A–K). We observed no major differences in the relative content of terpenes, vitamins and derivatives, and alkaloids (Figure 2B,E,F). BMPs 1392 and 262 exhibited the highest relative content of polyphenols and lipids, respectively (Figure 2C,I). Compared to 262 and K115, K1102 and 1392 had higher relative content of flavonoids, amino acids and derivatives, nucleotides and derivatives, and organic acids and derivatives (Figure 2A,G,H,J). In contrast, 262 and K115 exhibited the highest relative content of carbohydrates (Figure 2K). BMP 1392 had the highest relative content of phenylpropanoids, followed by K115, K1102, and 262 (Figure 2D).

We performed DAMs analysis to uncover DAMs between the different colored BMPs and insights into the differentially regulated pathways. Significant DAMs in pairwise comparison between groups were identified as shown in the volcano plots in Figure S3. In total, 154 (50 up-regulated in 262), 112 (77 up-regulated in 1392), 149 (71 up-regulated in 262), 183 (70 up-regulated in K115), 169 (132 up-regulated in 1392), and 176 (116 up-regulated in 262) significant DAMs were identified in pairwise comparison between 1392_vs_262, K1102_vs_1392, K1102_vs_262, K1102_vs_K115, K115_vs_1392, and K115_vs_262, respectively (Figure 2L). To identify key DAMs, we took K1102 as the control group, and we constructed a Venn diagram (Figure 2M). The results revealed thirty metabolites belonging to diverse classes overlapped (Figure 2M, Table S2). The classification of DAMs in pairwise comparisons of other pericarp color groups against K1102 revealed that flavonoids and organic acids and derivatives were the most differentially accumulated across the different colored BMPs (Figure 2N). The KEGG functional analyses assigned the DAMs mainly to the biosynthesis of secondary metabolites, phenylpropanoids, and amino acids (Figure S4).

3.3. Major Flavonoid Pigments in the Different Colored Bitter Melon Pericarps

Flavonoids represent the most diversified class of phenolic compounds and are primarily associated with the variation in colors of plant organs [24,27]. Therefore, we examined the distribution characteristics of differentially accumulated flavonoids in the four color groups of BMPs (Figure 3). Most differentially accumulated flavones and flavonols had the highest relative content in 1392 (very soft green), followed by 262 (grayish yellow), K1102 (grayish orange), and K115 (dark grayish cyan) (Figure 3A,B). Among flavonoids, astibin and tiliroside exhibited the highest relative content in 262, while isorhamnetin 3-O-glucoside and narcissoside showed the highest content in 1392 (Figure 3C). K115 exhibited the lowest relative content of phloretin and hesperetin 5-O-glucoside, while K1102 had the lowest relative content of hesperetin O-malonylhexoside and xanthohumol (Figure 3D). Anthocyanins are the main pigments of colored plant organs [20,21]. Our analysis revealed that K115 contained higher levels of four anthocyanins, including malvidin 3-O-glucoside (Oenin), rosinidin O-hexoside, cyanidin 3-p-hydroxybenzoylsophoroside-5-glucoside, and petunidin 3-O-glucoside (Figure 3E). We then inferred that they might represent the key pigments of dark grayish cyan BMP. Similarly, cyanidin and petunidin 3-O-glucoside might be the major pigment in soft green BMP, while cyanidin 3-p-hydroxybenzoylsophoroside-5-glucoside and pelargonidin 3-O-beta-D-glucoside might be the key pigment in grayish yellow BMP (Figure 3E). Cyanidin, cyanidin 3-p-hydroxybenzoylsophoroside-5-glucoside, and pelargonidin 3-O-beta-D-glucoside showed higher relative content in K1102 (grayish orange) (Figure 3E). Variations in the accumulation patterns of these six anthocyanins may represent the major underlying factor for BMP color changes (Figure 3E).

3.4. Comparative Transcriptome Sequencing and Differentially Expressed Genes (DEGs)

To provide more insights into the BMP pigmentation mechanisms, we subjected the four different colored pericarp samples to transcriptome sequencing. The RNA sequencing yielded 40.63–49.18 Mb of raw reads, with higher quality clean reads of 40.1–48.54 Mb (Table S3). The Q20, Q30, and GC contents varied from 95.66 to 96.65%, from 89.38 to 91.28%, and from 45.42 to 45.73%, respectively (Table S3). The total and unique mapping rates against the reference genome were 82.61–86.23% and 77.75–82.09%, respectively (Table S3). Samples of the same groups displayed significantly higher correlations, indicating the repeatability of the transcriptomics analysis (Figure S5B). The 1392 and 262 samples clustered closely on the PCA pot and exhibited high correlations (r = 0.92), suggesting no major transcriptional differences in grayish orange and soft green BMPs (Figure S5A,B). To verify the reliability of the RNA-seq data, we randomly selected ten genes for qRT-PCR analysis. As shown in Figure S6, the expressions of these genes by RNA-seq and qRT-PCR were very consistent (R² = 0.91), confirming that the transcriptome data are reliable.

To explore transcriptional changes, all DEGs in pairwise comparisons were filtered out. A total of 3324 (1641 up-regulated), 3037 (2007 up-regulated), and 4488 (2450 up-regulated) DEGs in pairwise comparisons between K1102_vs_262, K1102_vs_1392, and K1102_vs_K115, respectively, were identified (Figure 4A). The volcano plots are shown in Figure S7. A Venn diagram shows that 1363 DEGs were common in the three pairwise comparisons (Figure 4B). GO annotation and enrichment analysis assigned the DEGs to diverse cellular processes, such as tetrapyrrole binding, phenylpropanoid biosynthesis, secondary metabolism, etc. (Figure S8). The KEGG analysis results were supportive of those of GO (Figure 4C,D, and Figure S9). However, the flavonoid biosynthesis and tyrosine metabolism were also enriched (Figure 4C,D).

3.5. Expression Patterns of Flavonoid, Carotenoid, and Other Pigment Biosynthesis-Related DEGs

To identify key genes involved in BMP pigmentation, we examined the expression patterns of DEGs related to flavonoid biosynthesis (ko00941, ko00942, ko00943, ko00944), carotenoid biosynthesis (ko00906), and pigment biosynthesis (GO:0046148). Most of the flavonoid biosynthesis-related genes, such as phenylalanine ammonia-lyases, chalcone synthase (gene9135), chalcone isomerase (gene4079), flavonoid 3′-monoxygenase (gene4850), flavonol synthase (gene9134), vinorine synthase, anthocyanidin 3-O-glucosyltransferase (gene5091 and gene5092), etc., were highly induced in K1102, followed by 1392 and 262 (Figure 5A,B). gene18877 (tyrosine decarboxylase 1), gene219 (isoflavone 2′-hydroxylase), gene5598 (isoflavone reductase homolog PCBER), gene13201 (anthocyanin reductase, ANR), and gene11072 (isoflavone reductase homolog TP7) were significantly up-regulated in K115 compared to others (Figure 5A). Meanwhile, gene10002 (dihydroflavonol-4-reductase isoform X8) and gene8428 (coumaroyl-CoA:anthocyanidin 3-O-glucoside-6″-O-coumaroyltransferase 1) were also significantly up-regulated in 262 and 1392 compared to others (Figure 5A). Genes encoding phenolic glucoside malonyltransferase 2 were also induced in all colored BMPs except K1102 (Figure 5B).

As shown in Figure 6, most carotenoid catabolism-related genes, such as abscisic acid 8’-hydroxylase (gene11152, gene734, gene18471, gene5711, and gene5809), carotene epsilon-monooxygenase (gene10840), 9-cis-epoxycarotenoid dioxygenase (gene8322, and gene10129), zeaxanthin epoxidase (gene14628), beta-carotene hydroxylase 2 (gene8971), protein LUTEIN DEFICIENT 5 (gene6079), etc., were highly induced in other colored BMPs compared to the grayish orange (K1102). Most other pigment biosynthesis-related genes, including gene8173 (polyphenol oxidase), gene15171 (tetrapyrrole-binding protein), gene12160 (protein neoxanthin-deficiency), etc., were highly up-regulated in K115 compared to others (Figure S10).

3.6. Potential Candidate Transcription Factors (TFs) Regulating Bitter Melon Pericarp Pigmentation

TFs are key regulators of plant developmental and biochemical processes. To identify potential TF family genes that may govern BMP pigmentation, we screened out all DEGs encoding TFs. In total, we identified 337, 286, and 390 differentially expressed TFs in pairwise comparisons between K1102_vs_262, K1102_vs_1392, and K1102_vs_K115, respectively. MYB, AP2/ERF, bHLH, WRKY, and NAC were the major differentially expressed TFs (Figure 7A). We mainly focused on significantly up-regulated (Log2FC ≥ 3) MYBs, NACs, EP2/ERF, bHLHs, bZIPs, and MADSs with regard to their major roles in regulating flavonoid biosynthesis [20,24,28]. Most MYB (MYB30, MYB303, MYB36, MYB46, MYB20, etc.) and NAC (NAC7, NAC31, NAC26, NAC86, NAC90, and NAC98) genes were highly induced in other colored BMPs compared to K1102, particularly in 1392 (soft green) (Figure 7A,B). Three MYB (MYB2, MYB306, and MYB93) and two NAC (NAC104 and NAC43) genes exhibited high expression in K115 (dark grayish cyan) pericarp (Figure 7A,B). Meanwhile, MYB26 and NAC26 were the most induced in 262 (grayish yellow pericarp) (Figure 7A,B). Regarding other TFs, four ERFs (gene4902, ERF054, ERF2, and ERF106), one bHLH (bHLH97), and three MADSs (AGL12, AGL8, and EJ2) were highly induced in K115 (Figure S9). Two bHLHs (MYC4 and bHLH94) and two bZIPs (bZIP44 and bZIP61) were highly up-regulated in 262 (Figure S11). In 1392, two bHLHs (AIG1 and bHLH126) were also highly induced (Figure S11).

To facilitate an overview of the main anthocyanins and flavonoid structural and regulatory genes involved in BMP color variation, we constructed a diagram that integrates the main results by metabolomics and transcriptomics. These key findings are presented in Figure 8.

4. Discussion

Fruit color is one of the most critical traits affecting fruits’ commercial and quality values. Hence, deciphering and understanding the molecular regulation of color formation in fruit crops is essential to uncover key metabolic and genetic resources for breeding novel fruit varieties of high consumer preference and economic value. Accordingly, this study took advantage of advanced metabolomics and transcriptomics tools to provide insights into the mechanisms underlying the pericarp color of four bitter melon genotypes, including K1102 (grayish orange), 262 (grayish yellow), 1392 (very soft green), and K115 (dark grayish cyan).

The metabolomics analysis revealed a great variation in the metabolite profiles of the four different colored BMPs. K1102, 262, 1392, and K115 samples clustered separately on the PCA and HCA plots and accumulated different levels of important metabolites, such as alkaloids, flavonoids, polyphenols, lipids, amino acids and derivatives, carbohydrates, etc. These results indicate that fruit color may significantly influence bitter gourd fruits’ quality and medicinal values. Fruit of different colors may display considerable variations in pharmaceutical and therapeutical attributes. Supportively, Oude Griep et al. have demonstrated that the ability of fruits and vegetables to treat coronary heart diseases varies depending on their color [46]. Furthermore, these findings corroborate with reports in other fruit crops [12,14,47]. All DAMs in pairwise comparisons between groups were identified, providing fundamental data for in-depth exploration of the variation in bioactive compounds among the four different colored BMPs. Thirty key overlapped DAMs were filtered out and could represent discriminatory metabolic markers to differentiate between BMPs of different colors. Metabolites’ diversity and variability are driven by different expression of genes regulating metabolic pathways. Concordantly, transcriptome sequencing and analysis revealed high numbers of significant DEGs across the four different colored BMPs. These results show that metabolism regulation in the four different colored BMPs is different, and further interest should be given to pathways related to the fruit quality traits.

Anthocyanins are the main metabolic compounds that determine the color of fruits [20,21]. Our analysis unveiled that variation in the levels of six anthocyanins, including malvidin 3-O-glucoside, petunidin 3-O-glucoside, rosinidin O-hexoside, cyanidin, cyanidin 3-p-hydroxybenzoylsophoroside-5-glucoside, and pelargonidin 3-O-beta-D-glucoside, is the major driving factor of fruit color variation in bitter melon. Compared to K115, the other BMPs exhibited higher relative content of most differentially accumulated flavones, flavonols, flavanones, and flavonoids. Accordingly, we found that many flavonoid pathway structural genes were differentially regulated in the four different colored BMPs. Besides anthocyanins, carotenoids, xanthophylls, chlorophyll [25], and other flavonoids [24] contribute to the final color appearance of fruits [14]. For instance, studies have revealed that variations in the composition and amounts of carotenoids and chlorophyll contribute to the yellow, orange, red [26], light green, green, dark green, and white peel and flesh colors in cucumber [13,14]. The examination of DEGs related to carotenoid metabolism showed that the carotenoid catabolism-related genes were not induced in K1102, 262, and 1392. Moreover, tetrapyrrole-binding protein (related to chlorophyll accumulation) genes were significantly down-regulated in K1102 and 262. Taken together, these results show that the grayish orange and grayish yellow colors of K1102 and 262 fruits, respectively, may be due to co-pigmentation between anthocyanins, other flavonoids, carotenoids, and xanthophylls, with cyanidins and pelargonidins being the major pigments. Pelargonidin is known for its natural orange color [20]. The very soft green color of 1392 may result from a co-pigmentation between anthocyanin (mainly cyanidin and petunidin 3-O-glucoside), other flavonoids, chlorophyll, xanthophylls, and carotenoids. Further quantification of chlorophyll and carotenoid contents of different colored BMPs is required to confirm these statements.

K115 displayed a higher content of malvidin 3-O-glucoside, rosinidin O-hexoside, petunidin 3-O-glucoside (dominant pigments), cyanidin 3-p-hydroxybenzoylsophoroside-5-glucoside, and low relatively contents of other flavonoids. Tetrapyrrole-binding protein, carotenoid catabolism, and ABA metabolism-related genes were highly induced in K115 compared to other BMPs. These results indicate that the dark grayish cyan color of K115 may be due to a co-pigmentation between these four anthocyanins and chlorophyll, with malvidin 3-O-glucoside, rosinidin O-hexoside, and petunidin 3-O-glucoside being the dominant pigments. Petunidins have been identified as the major pigments in black currants and black goji [23], while malvidins were found to be the dominant pigments in darker red wine [20]. Furthermore, we found that the genes gene8173 (polyphenol oxidase) and gene13201 (ANR) were up-regulated in K115 compared to other BMPs. These results infer that the browning reaction (polymerization of proanthocyanins) may contribute to the fruit color development of K115. The enzymatic browning reaction catalyzed by polyphenol oxidases (PPOs) influences the color quality of vegetables and fruits [48]. The significant induction of ABA (abscisic acid) metabolism-related genes in K115 suggests that ABA may promote anthocyanin biosynthesis and accumulation, leading to their predominance in the fruit pericarp. In fig (Ficus carica L.), it is shown that ABA stimulates anthocyanin synthesis and accumulation to enhance fruit color during repining [49]. Further studies on ABA-promoting anthocyanin accumulation in BMPs are needed for biochemical control of better gourd fruit color perspectives.

TFs, particularly R2R3-MYBs, bHLHs, and WD40 proteins, play critical roles in regulating the expressions of flavonoid structural genes at the transcriptional level [20,24,28]. We have identified potential candidate MYB, NAC, EP2/ERF, bHLH, bZIP, and MADS family genes that might regulate the levels of anthocyanin compounds in different colored BMPs. These genes should be subjected to functional genomics studies in order to dissect the regulatory network of flavonoid biosynthesis and accumulation in bitter melon fruits. Notably, particular attention should be given to the candidate MYB family genes due to their recognized major regulation function of anthocyanins in plant organs [29,30,31,50,51]. Attention should also be paid to NAC family candidate genes, particularly NAC104 and NAC43, which showed significantly higher expression in K115 compared to other BMPs. NAC TFs ANAC078 and MdNAC52 were found to regulate anthocyanin and other flavonoid synthesis and accumulation in Arabidopsis and apple, respectively [52,53]. The flavonoid pathway-related differentially expressed structural genes identified may also represent an important resource for modulating the composition and concentration of flavonoids in bitter melon fruits.

5. Conclusions

The present study identified key metabolites and insights into the molecular mechanisms underlying the pericarp color variation in bitter melon. The variation in the content of six anthocyanins, including malvidin 3-O-glucoside, petunidin 3-O-glucoside, rosinidin O-hexoside, cyanidin, cyanidin 3-p-hydroxybenzoylsophoroside-5-glucoside, and pelargonidin 3-O-beta-D-glucoside, is the major driving factor of BMP color differences. The grayish orange pericarp color may be a co-pigmentation between cyanidin, cyanidin 3-p-hydroxybenzoylsophoroside-5-glucoside, pelargonidin 3-O-beta-D-glucoside (major pigments), other flavonoids, and carotenoids. Meanwhile, co-pigmentation between cyanidin 3-p-hydroxybenzoylsophoroside-5-glucoside, pelargonidin 3-O-beta-D-glucoside (major pigments), other flavonoids, and carotenoids may be responsible for grayish yellow pericarp color. The soft green pericarp color may result from a co-pigmentation between petunidin 3-O-glucoside, cyanidin (major pigments), other flavonoids, chlorophyll, xanthophylls, and carotenoids. Malvidin 3-O-glucoside, rosinidin O-hexoside, and petunidin 3-O-glucoside are the dominant pigments in K115 (dark grayish cyan), but cyanidin 3-p-hydroxybenzoylsophoroside-5-glucoside and chlorophyll may participate in the color. MYBs, NACs, MADSs, bHLHs, EP2/ERFs, and bZIPs TF family genes may play critical roles in regulating BMP pigmentation. Candidate biosynthesis and regulatory genes were identified for future studies. Among them, gene13201 (anthocyanin reductase), gene8173 (polyphenol oxidase), gene2136 (NAC43), gene19593 (NAC104), and gene15171 (tetrapyrrole-binding protein) need to be paid attention to when dissecting dark grayish cyan color development. Our findings will facilitate the molecular dissection of the regulatory network of pericarp pigmentation in bitter melon and contribute to breeding varieties with high nutritional value.