Flesh color assessment and carotenoid contents variation during fruit development of five watermelon genotypes
Watermelon flesh features at different developmental stages have been shown in Fig. 1. We determined the color space values to confirm flesh color variations. At the 10 days after pollination (DAP), all fruits were white flesh, and there is no significant difference in color space parameters between different genotypes (Additional file1: Table S1). At the 20 DAP, The fruit's flesh presents varying degrees of white, pink or yellow owe to carotenoids accumulation. 20 DAP is a critical period for the rapid accumulation of pigments. At the 34 DAP, the fruits were matured and the flesh has vivid colors except for the white flesh genotype. Significant differences of L*, a*, b*, and Chroma (C) can be observed in different flesh-colored fruits at 20 DAP and 34 DAP. The difference in flesh color appeared at the 20 DAP and was more pronounced at 34 DAP in this study.
Transcript sequencing of watermelon flesh with different colors
To explore the potential molecular mechanisms underlying the flesh coloration during the fruit development of 5 watermelon genotypes, RNA-Sequencing analysis was conducted on fruit flesh to generate transcriptome profiles. Samples of fruit flesh at three critical stages (10 DAP, 20 DAP, and 34 DAP) were obtained from five genotypes (Fig. 1). All samples were analyzed as three independent biological replicates.
In total, 45 libraries were constructed and analyzed. After removing low quality reads, the average reads number of per library was 51.9 million, with an average GC content of 44.15% (Additional file 1: Table S3). The RNA-Seq reads were aligned with the reference map of the watermelon (97103) genome (http://cucurbitgenomics.org/organism/1), using HISAT(version 2.0.4) [27]. More than 97% of the total clean reads had Phred-like quality scores at the Q20 level (Additional file 1: Table S3). Ultimately, 24,794 genes (including 1354 novel genes) were identified by Cufflinks v2.1.1 [28]. The numbers of transcripts identified in each sample were expressed in FPKMs. Approximately 39.08% of expressed genes were in the 0-1 FPKM range, and 13.48% of expressed genes showed high expression levels (above 60 FPKM) (Additional file 1: Table S4). Genes with normalized reads lower than 1 FPKM were removed from the subsequent analysis. The gene expression levels among different experimental groups were compared in Fig. 3a. The expression patterns among biologically repeated samples were highly consistent (Fig. 3b) and the correlation coefficient was close to 1 (Additional file 1: Table S5). Therefore, this high-quality RNA-Seq data provided a solid foundation for identifying key genes participating in carotenoid syntheses during watermelon fruits development.
Identification of differentially expressed genes in five genotypes
We conducted a pairwise comparison at three developmental stages of five genotypes to identify the genes correlating with the flesh color. The DEGs were screened with FDR<0.05, |log2(FoldChange)| >1 as a threshold, the numbers and lists of significantly differentially expressed genes (up-regulated and down-regulated) of each pairwise comparison were shown in Additional file 1: Table S6 and Additional file 3-6: Dataset 1-4. 16781 genes were differentially expressed between at least one comparison.
The global hierarchical clustering (Fig. 4a) and principal component analysis (Fig. 4b) were performed based on the FPKM values for all the DEGs. The results revealed that 45 samples could be generally assigned into three main groups corresponding to development stages based on gene expression patterns. The samples from 10 DAP were distinctly clustered as one group, the samples from 20 DAP and 34 DAP were clustered as another group except for the samples at 34 DAP of red flesh (Fig. 4a), suggesting that the expression patterns of most DEGs were consistent during fruit development of different genotypes. In particular, we can see the differences in gene expression patterns between genotypes become clearer at 34 DAP compared to that of 10 DAP and 20 DAP in the PCA analysis (Fig. 4b). At 20 DAP and 34 DAP, the white, pink and orange genotypes were clustered together, while the yellow and red genotypes were separated from each other (Fig. 4b). Three biological repeats for red color at 34 DAP were not clustered together (Fig. 4b), which may be caused by environmental differences during cultivation. The difference in the overall gene expression patterns indicates that there must be a set of differentially expressed genes associated with the difference of flesh coloration in watermelon.
At the early development stage (10 DAP), 5318 significantly differentially expressed genes were identified (Fig. 5a, Additional file 3: Dataset 1). Specifically, 510,262,588,349 differentially expressed genes were identified in the red flesh genotype as compared to the pink, orange, yellow, and white flesh genotypes. The numbers of other pairwise comparisons were listed in Additional file 1: Table S6. The candidate gene linked to fruit shape, ClFS1 (Cla011257), was differentially expressed in different genotypes at 10 DAP, consistent with the previous research [29]. Cla019403 encode xyloglucan endotransglucosylase is related to plant cell growth [30] and highly expressed at this stage (Additional file 2: Fig. S2a, Additional file 1: Table S7). An auxin response factor (ARF, Cla009800), a growth-regulating factor 5 (GRF, Cla006802), an auxin-induced SAUR-like protein (Cla016617), were associated with the fruit development and expansion [31] and highly expressed at early developmental stages (Additional file 2: Fig. S2a, Additional file 1: Table S7). There are less DEGs at 10 DAP compared to later stages, maybe due to little differences among five genotypes at this stage.
At the pigment accumulation stage (20 DAP), 11814 significantly differentially expressed genes were identified (Fig. 5b, Additional file 4: Dataset 2). Specifically, 2498,4830,3123,4876 differentially expressed genes were identified in the red flesh genotype compare to the pink, orange, yellow, and white flesh genotypes. The numbers of other pairwise comparisons were listed in Additional file 1: Table S6. The geranylgeranyl pyrophosphate synthase (Cla015797), phytoene synthase protein (Cla005425), phytoene desaturase (Cla010898), carotenoid isomerase (Cla017593), lycopene cyclase (Cla016840), violaxanthin de-epoxidase-related protein (Cla007759), 9-cis-epoxycarotenoid dioxygenase (Cla015245) were involved in carotenoid biosynthesis and differentially expressed among 5 genotypes (Additional file 2: Fig. S2b, Additional file 1: Table S7). Two AP2-EREBPs (Cla000701, Cla017389) and two bHLHs (Cla020193, Cla022119) were differentially expressed in 5 genotypes (Additional file 2: Fig. S2b, Additional file 1: Table S7). The expression level of Cla017389 was positively related to the contents of lycopene (Pearson's r =0.85) and γ-Carotene (r =0.69) in 15 experimental groups. The expression level of Cla020193 was negatively correlated with the contents of phytofluene (r =-0.60) and phytoene (r =-0.57), the expression level of Cla022119 was negatively correlated with the contents of phytofluene (r =-0.59) and phytoene (r =-0.58). These genes maybe the potential color regulators in watermelon according to the studies in other crops [2, 32].
At the maturity stage (34 DAP), 10779 significantly differentially expressed genes were identified (Fig. 5c, Additional file 5: Dataset 3). Specifically, 2097,2572,2429,3316 differentially expressed genes were identified in red flesh genotype compared to pink, orange, yellow, and white flesh genotypes. The numbers of other pairwise comparisons were listed in Additional file 1: Table S6. Most of the carotenoids pathway genes and many TFs are differentially expressed among 5 genotypes at this stage. The geranylgeranyl reductase (Cla003139, Cla019109), geranylgeranyl pyrophosphate synthase (Cla015797, Cla020121), phytoene synthase protein (Cla005425, Cla009122), phytoene desaturase (Cla010898), carotenoid isomerase (Cla017593, Cla011810), lycopene cyclase (Cla005011, Cla017416, Cla016840), 9-cis-epoxycarotenoid dioxygenase (Cla015245, Cla009779, Cla005404, Cla005453, Cla019578), beta-carotene hydroxylase (Cla011420, Cla006149), zeta-carotene desaturase (Cla003751), zeaxanthin epoxidase (Cla020214), and many TFs (AP2-EREBPs, MADSs, MYBs, G2-likes, NACs, AUXs), were differentially expressed at 34 DAP (Additional file 2: Fig. S2c, Additional file 1: Table S7). Cla015245 and Cla005404 (9-cis-epoxycarotenoid dioxygenase) were highly expressed in white flesh may lead to the degradation of xanthophyll and colorless flesh. Still, more notably, Cla015245 and Cla005404 have the highest expression level in the mature pink fruit (Additional file 2: Fig. S2c, Additional file 1: Table S7) which possess a vivid flesh color. The results indicate that different genotypes have different color formation mechanisms. Cla017389 and Cla015515 (AP2-ERFBP) were highly expressed in red and pink fruit, respectively. Cla017389 and Cla015515 were homologous to ethylene-responsive transcription factor RAP2-2 (E-value: 1.4e-28 and 3.0e-82) that are involved in the regulation of carotenoid biosynthesis in Arabidopsis thaliana [33]. More importantly, there was a positive correlation between the expression level of Cla017389 and the contents of lycopene (r =0.85), as mentioned above. The transcription factor bHLH is related to the carotenoid metabolism in tomato [34], papaya [35], and citrus [36]. In this study, the expression of Cla006599 and Cla022119 (bHLH) were decreased at 20 and 34 DAP as compared to 10 DAP, and their expression patterns were similar to CpbHLH1/2 that regulate carotenoid biosynthesis in papaya [35] (Additional file 2: Fig. S2c, Additional file 1: Table S7). Further analysis showed that the expression level of Cla006599 was negatively correlated with the contents of phytofluene (r = -0.56) and phytoene (r = -0.62). The expression level of Cla022119 was also negatively related to the contents of phytofluene and phytoene as described above. Here we also note that a zinc finger CCCH domain-containing protein (Cla007686) showed a significant increase in expression level (~3 times) in red fruit at the ripening stage than that of early stage, (Additional file 1: Table S7), and the expression level was significantly associated with the lycopene content in 15 groups (r = 0.81). Five MYB related genes (Cla020633, Cla007790, Cla009263, Cla017995, and Cla019223) were differentially expressed in 5 genotypes at 34 DAP (Additional file 2: Fig. S10a, Additional file 1: Table S11). The content of phytoene was positively correlated to the gene expression levels of Cla009263 (r = 0.67) and Cla017995 (r = 0.62).
For each genotype, 6430, 9019, 9605, 9147, and 10881 developmental DEGs were obtained in red-, pink-, orange-, yellow-, and white-fleshed watermelon genotypes, respectively (Fig. 5d, Additional file 6: Dataset 4). These results indicate that a large number of genes were involved in the regulation of fruit development. More genes are differentially expressed in the white flesh watermelon, suggesting a very complicated regulatory network of gene expression in this genotype. We also identified some differentially expressed genes related to fruit development using comparative transcriptome analysis. A cytokinin dehydrogenase gene (Cla022463) was highly expressed at 10 DAP (Additional file 1: Table S7), may be contribute to the early fruit development. The pyrabactin resistance 1-like protein (PYL8) can regulate the plant growth and stress responses by mediate ABA signaling in Arabidopsis [37]. Here, we found the expressions of four abscisic acid receptor PYL8 (Cla004235, Cla004904, Cla015009, Cla021167) were significantly different in 5 watermelon genotypes (Additional file 2: Fig. S2d, Additional file 1: Table S7). The gene expression level of Cla004904 had a negative correlation with the contents of antheraxanthin (r = - 0.61), violaxanthin (r = - 0.59) and neoxanthin (r = - 0.58). This gene may be involved in carotenoid degradation and abscisic acid metabolism in fruit development and ripening.
Clustering of DEGs into six groups based on gene expression patterns
Based on the pattern of gene expression, the 16781 DEGs are grouped into 6 different subclusters using the h-cluster clustering approach (Fig. 6a, Additional file 1: Table S8). The genes grouped in the same cluster shared a similar expression pattern and have a similar function or participate in the same biological process. In subcluster 1, 8931 genes were displaying relatively stable expression levels across different stages and genotypes. In subclusters 2 and 3, genes expression exhibited a general downward trend according to the developmental stages. Genes in subclusters 4, 5, and 6 exhibited a general upward trend according to the developmental stages (Fig. 6a).
To functionally characterize the biological roles of these DEGs in subcluster 1, GO enrichment analyses were performed. Almost the same proportion of genes were enriched into biological process, cellular component, and molecular function, GO terms related to various basic life activities, such as binding, cellular macromolecule metabolic process, intracellular, and cell (Additional file 2: Fig. S3a, Additional file 7: Dataset 5). In addition, DEGs are enriched to the spliceosome, RNA transport, and ribosome biogenesis in eukaryotes pathways by Kyoto Encyclopedia of Genes and Genomes (KEGG) analysis (Fig. 6b, Additional file 7: Dataset 5). Previous studies showed that HY5 (ZIP) is involved in chloroplast biogenesis in Arabidopsis and tomato [38, 39]. Here, the transcription factor Cla016581, Cla017361, Cla002873, and Cla021184 (ZIP) were differentially expressed in different samples (Additional file 1: Table S7). The gene expression of Cla002873 was positively correlated to the contents of γ-Carotene (r = 0.77) and lycopene (r = 0.85). The gene expression of Cla017361 was also positively correlated to the contents of γ-Carotene (r = 0.66) and lycopene (r = 0.72). Also, two GLK2 TFs (Cla010265, Cla020369), were differentially expressed in different samples (Additional file 1: Table S7). The gene expression of Cla010265 was positively correlated to the contents of lycopene (r = 0.77) and total carotenoids (r = 0.74). The gene expression of Cla020369 was also related to the content of lycopene (r = 0.74). The transcription factor Cla010815 (MADS) is homologous to SlMADS1, which plays an important role in fruit ripening as a repressive modulator in tomato [40] (E-value: 3e-76, identity: 77.24%). Its expression level was highest in red flesh fruit at 20 DAP and 34 DAP (Additional file 1: Table S7) and related to the lycopene content (r = 0.55). The transcription factor Cla009725 (MADS) is homologous to CsMADS6, which was coordinately expressed with fruit development and coloration in citrus [41] (E-value: 4e-108, identity: 68.62%), its expression level was negatively correlated with the content of phytoene (r = - 0.86), γ-carotene (r = - 0.82), lycopene content (r = - 0.73), and total carotenoids (r = - 0.85). The zinc finger CCCH domain-containing protein (Cla007686) mentioned above was also in subcluster 1.
3282 genes have a slightly higher expression at 20 DAP and 34 DAP than 10 DAP in subcluster 2 (Fig. 6a). GO terms such as single-organism metabolic process, small molecular metabolic process, and organonitrogen compound metabolic process were enriched (Additional file 2: Fig. S3b, Additional file 7: Dataset 5). KEGG enrichment analysis showed that the most significantly enriched pathways are proteasome, biosynthesis of amino acids, carbon metabolism pathway, TCA-cycle, glycolysis/gluconeogenesis proteasome, and pyruvate metabolism pathway. The proteasome pathway is containing 26S protease regulatory subunit genes and proteasome subunit type genes. Acetyl-CoA carboxylase biotin carboxylase, pyruvate kinase, and malate dehydrogenase are in the pyruvate metabolism pathway (Fig. 6b, Additional file 7: Dataset 5). The transcription factor Cla000691 is homologous to SlMADS1 that plays as a repressive modulator in tomato fruit ripening [40] (E-value: 2e-87, identity: 64.93%) and highly expressed in the pink-fleshed fruits at later development stages (Additional file 1: Table S7).
Subcluster 3 represents genes that were highly expressed at 34 DAP, and the range of change was more obvious than that of subcluster 1 (Fig. 6a). 754 DEGs in this cluster mainly allocated into molecular function and biological process according to GO term analysis, with 310 and 287 DEGs were classified into the metabolic process and catalytic activity (Additional file 2: Fig. S3c, Additional file 7: Dataset 5). Notably, the DEGs were significantly involved in pathways associated with photosynthesis-antenna proteins biosynthesis, photosynthesis, protein processing in endoplasmic reticulum, and biosynthesis of unsaturated fatty acids (Fig. 6b, Additional file 7: Dataset 5). Cla006149 and Cla011420 (beta-carotene hydroxylase), Cla009122 (phytoene synthase), and Cla009779 (9-cis-epoxycarotenoid dioxygenase) are related to the carotenoid pathway and differentially expressed in 5 genotypes (Additional file 2: Figs. S2b and S2c, Additional file 1: Table S7). Some MYBs, AP2-ERFBPs, bHLHs, NACs, and WRKYs may be important regulators in fruit development and ripening [2]. In subcluster 3, we observed some differentially expressed transcription factors (Additional file 7: Dataset 5). Two MYBs (Cla018631 and Cla006739) and two WRKYs (Cla002243 and Cla002084) were highly expressed in yellow and pink flesh, respectively (Additional file 2: Fig. S2e, Additional file 1: Table S7). The gene expression levels of Cla018631 and Cla006739 were correlated with the contents of antheraxanthin and violaxanthin (r: Cla018631 – antheraxanthin = 0.88, Cla018631 - violaxanthin = 0.69, Cla006739 – antheraxanthin = 0.85, Cla006739 - violaxanthin = 0.69). The gene expression of Cla002243 was negatively correlated with the content of violaxanthin (r = - 0.58), but there was a positive correlation between the expression of Cla002084 and the content of phytoene (r = 0.64).
There were 2274, 1174, and 366 genes in subclusters 4, 5, and 6, respectively. These genes were highly expressed at 10 DAP and decreased to low gene expression levels at the later stages with different change magnitudes (Fig. 6a). GO enrichment analysis of subcluster 4 indicated that biological processes were most enriched. (Additional file 2: Fig. S3d, Additional file 7: Dataset 5). KEGG analysis showed that the genes significantly involved in the pathway of plant hormone signal transduction, such as signal transduction histidine kinase (Cla000685, Cla005808), auxin responsive protein (Cla003635), Ein3-binding f-box protein (Cla020970), and so on (Fig. 6b, Additional file 7: Dataset 5). The transcription factor Cla019630 (MADS) in subcluster 4 was homologous to CsMADS6 that coordinately expressed with citrus fruit development and coloration [41] (E-value: 2e-98; Identity: 75.27%), its gene expression was negatively related to content of violaxanthin (r = - 0.57). GO enrichment of subcluster 5 genes assigned to the biological process and molecular function, such as protein phosphorylation, protein kinase activity (Additional file 2: Fig. S3e, Additional file 7: Dataset 5). Genes are enriched into plant hormone signal transduction, alanine, aspartate, glutamate metabolism, and others by KEGG (Fig. 6b, Additional file 7: Dataset 5). In subcluster 6, the enriched Go terms were predominantly related to molecular function and biological processes, such as enzyme inhibitor activity and endopeptidase regulator activity (Additional file 2: Fig. S3f, Additional file 7: Dataset 5). KEGG enrichment was mostly related to the pathways of plant hormone signal transduction, phenylpropanoid biosynthesis, and phenylalanine metabolism (Fig. 6b, Additional file 7: Dataset 5). Cla019806, Cla0004102, Cla002975, and Cla016617 were involved in hormone synthesis, and highly expressed at early stages of fruit development as compared to later fruit developmental stages (Additional file 2: Fig. S2f, Additional file 7: Dataset 5, Additional file 1: Table S 7).
Co-expression network analysis identified carotenoid-related DEGs
To identify the potential genes (structural genes and putative transcription factors) highly associated with different kinds of carotenoid accumulation. The carotenoids content in each sample was used as phenotypic data and 16781 DEGs were used to perform the weighted gene co-expression network analysis (WGCNA).
A sample dendrogram and trait heatmap was constructed to illustrate the expression of each phenotypic parameter at different developmental stages (Additional file 2: Fig. S4a). The best parameter value determination for module construction is 7.7 for this dataset (Additional file 2: Fig. S4b). A total of 40 distinct co-expression modules were formed according to the pairwise correlations of gene expression across all samples and co-expression patterns of individual genes, as shown in the cluster dendrogram (Additional file 2: Fig. S5a). Moreover, a network heatmap of all the DEGs in gene-modules was also drawn to exhibit the correlation between modules (Additional file 2: Fig. S5b). Notably, 6 co-expression modules (indicated with red underlines) have a high positive correlation with most carotenoids (Fig. 7a), meaning that the genes in these modules play an important role in the carotenoid accumulation.
'Yellow' module contains 846 Genes (including 34 TFs) (Additional file 1: Table S9) exhibited a stronger positive relationship with zeaxanthin (correlation coefficient, r=0.91), neoxanthin(r=0.91), antheraxanthin (r=0.84), violaxanthin (r=0.79), phytofluene (r=0.81), apocarotenal (r=0.76), β-cryptoxanthin (r=0.68), lutein (r=0.69), and α-carotene (r=0.62) (Fig. 7a). In this module, a set of genes were related to the cellular metabolic process, and involved in pyruvate metabolism and proteasome pathway were identified (Additional file 2: Fig. S6). Pyruvate is an important mediator of carbohydrate, fat, and protein metabolism, and participates in several important metabolic functions in vivo. The proteasome is related to the regulation of carotenoid content in tomato [42]. This module contains differentially expressed genes between the yellow-/orange-fleshed genotypes and the red-/pink-/white-fleshed genotypes (Additional file 2: Fig. S7a), and maybe the important factors involved in yellow pigments accumulation. According to gene function annotation, Cla005011 is lycopene beta-cyclase in watermelon [23]. Cla003751 encodes a zeta-carotene desaturase, Cla020214 encodes zeaxanthin epoxidase. Hence, they were involved in the carotenoid pathway (Additional file 2: Fig. S2c, Additional file 1: Table S7 and S9). Cla018406 (a chaperone protein dnaJ-like protein) was highly expressed in orange and yellow flesh and its gene expression level was related to the β-Carotene content (r = 0.71) and neoxanthin content (r = 0.93). The Cla014416 (plastid-lipid-associated protein, ClPAP) is homologous to SlPAP(NP_001234183.1) that affect carotenoid content in tomato [43] (E-value: 1e-145, Identity: 68.67%). The expression level of Cla014416 was higher in yellow-/orange-flesh than in red-/pink-/white-flesh genotypes used in this study (Additional file 2: Fig. S2g, Additional file 1: Table S7 and S9). This result is different from the previous report that ClPAP highly expressed in red/orange than yellow /white color genotypes [44], possibly due to different genotypes used in these two studies. Cla000655 (encodes a cytochrome P450) is homologous to the protein lutein deficient 5 (CYP97A3), which involved in the biosynthesis of xanthophylls in Arabidopsis thaliana [45] (E-value: 6.5e-267, Identity: 80.47%). Cla010839 is homologous to 15-cis-zeta-carotenoid isomerase in Arabidopsis thaliana [46] (E-value: 1.1e-135, Identity: 67.43%). Cla018347 (encodes a cytochrome P450) is related to carotenoid epsilon-monooxygenase (CYP97C1) in Arabidopsis thaliana [47] (E-value: 1.2e-233, Identity: 76.20%). Cla000655, Cla010839, and Cla018347 were highly expressed in yellow and orange color fruits (Additional file 2: Fig. S2g, Additional file 1: Table S7 and S9) and maybe involved in the biosynthesis of xanthophylls in watermelon. The MYB transcription factor can regulate the carotenoid contents in Mimulus lewisii flowers [48], Cla013280 and Cla010722 belong to the MYB family and highly expressed in yellow-fleshed fruits (Additional file 2: Fig. S2g, Additional file 1: Table S7 and S9). The gene expression of Cla013280 was positively related to the contents of antheraxanthin (r = 0.85) and violaxanthin (r = 0.81). The gene expression of Cla010722 was also related to the contents of antheraxanthin (r = 0.75) and violaxanthin (r = 0.75). The hub genes linked to this module were further analyzed using Cytoscape cytoHubba (Fig. 7b), the ATP synthase protein I, (Cla013542), cysteine desulfuration protein SufE (Cla003340), membrane protein (Cla003760), and others were identified as hub genes responsible for yellow color formation in watermelon (Additional file 2: Fig. S8, Additional file 1: Table S9 and S10).
The 'dark-red' module containing 111 genes was positively associated with the contents of antheraxanthin and violaxanthin, having a correlation coefficient of 0.76 and 0.71 respectively. Heatmaps (Additional file 2: Fig. S7b) showed that the 'dark-red' module-specific genes were over-represented in samples (yellow, orange, white) rich in antheraxanthin and violaxanthin. Cla004704 encodes a photosystem II oxygen evolving complex protein PsbP, Cla005429 encode an oxygen-evolving enhancer protein 2, chloroplastic, PsbP. Cla004746 encodes a chlorophyll a-b binding protein 6A (Additional file 2: Fig. S2g, Additional file 1: Table S7 and S9). Cla021635 encodes a photosystem I reaction center subunit II, rank as the top hub gene in this module (Fig. 7c, Additional file 2: Fig. S8, Additional file 1: Table S10). Many genes in this module are also related to chloroplast or photosystem I, II (Additional file 1: Table S7).
The 'mediumpurple 3' module, with 32 identified genes, was highly correlated to the contents of α-carotene, violaxanthin, neoxanthin, lutein, and zeaxanthin with the correlation coefficient of 0.55, 0.70, 0.72, 0.66, and 0.73 respectively (Fig. 7a). The heatmap was shown in (Additional file 2: Fig. S7c). Cla005637, Cla017046, and Cla011297 were identified as candidate hub genes for this module (Additional file 2: Fig. S8, Additional file 1: Table S10). The 'black' module was specific to the contents of lycopene (r=0.57) and lutein (r=0.58). The 'steelblue' module was specific to the contents of γ-Carotene (r=0.55) and phytofluene (r=0.55), respectively (Fig. 7a). The transcription factor Cla019630 (MADS) gene mentioned above was also in the 'saddlebrown' module, which was correlated to the contents of zeaxanthin (r=0.55), neoxanthin(r=0.56), and antheraxanthin (r=0.51) (Fig. 7a). Their hub genes were listed in Additional file 1: Table S10 and expression levels are shown in Additional file 2: Fig. S8.
By WGCNA, we found that most carotenoid pathway genes are included in the yellow module, which was positively correlated with the yellow and orange color samples (Additional file 2: Fig. S7a). Through network analysis, we identified hub genes in eight main modules for watermelon flesh carotenoid content (Additional file 1: Table S10).
Validation of the expression of key DEGs by qRT-PCR
Twenty-one DEGs were used for qRT-PCR analysis to verify the quality of RNA-Seq data. We found a strong correlation between the RNA-Seq and qRT-PCR data (r = 0.90 ~ 0.99, the correlation was calculated separately for each gene), indicating the reliability of our transcriptome data (Additional file 2: Fig. S9) and the candidate gene from this study can be used for further functional verification.