Identification of key gene networks controlling organic acid and sugar metabolism during watermelon fruit development by integrating metabolic phenotypes and gene expression profiles

The organoleptic qualities of watermelon fruit are defined by the sugar and organic acid contents, which undergo considerable variations during development and maturation. The molecular mechanisms underlying these variations remain unclear. In this study, we used transcriptome profiles to investigate the coexpression patterns of gene networks associated with sugar and organic acid metabolism. We identified 3 gene networks/modules containing 2443 genes highly correlated with sugars and organic acids. Within these modules, based on intramodular significance and Reverse Transcription Quantitative polymerase chain reaction (RT-qPCR), we identified 7 genes involved in the metabolism of sugars and organic acids. Among these genes, Cla97C01G000640, Cla97C05G087120 and Cla97C01G018840 (r2 = 0.83 with glucose content) were identified as sugar transporters (SWEET, EDR6 and STP) and Cla97C03G064990 (r2= 0.92 with sucrose content) was identified as a sucrose synthase from information available for other crops. Similarly, Cla97C07G128420, Cla97C03G068240 and Cla97C01G008870, having strong correlations with malic (r2 = 0.75) and citric acid (r2 = 0.85), were annotated as malate and citrate transporters (ALMT7, CS, and ICDH). The expression profiles of these 7 genes in diverse watermelon genotypes revealed consistent patterns of expression variation in various types of watermelon. These findings add significantly to our existing knowledge of sugar and organic acid metabolism in watermelon.


Introduction
Watermelon is the fifth most consumed fleshy fruit worldwide, with a global yield of 118 million tons in 2017 1 . It occupies 7% of the total cultivated land of fruits and vegetables. Fruits of watermelon are rich in water and nutrients (amino acid sugars, carotenoids, lycopene, organic acids, etc.) 2,3 . These phytochemicals make watermelon one of the most nutritious fruits, providing substantial nutritional supplementation to the human diet [4][5][6][7] . Sugars and organic acids have a strong influence on organoleptic fruit quality and are crucial components involved in the development of fruit flavor 8 . Unlike staple food crops, for which yield is the ultimate breeding goal, taste and aroma are more vital traits for watermelon, which are both determined by the metabolite composition of the fruit. During the development process, watermelon fruits undergo numerous biochemical changes, including changes in sugar metabolism, organic acid and pigment accumulation, fruit softening, flavor, and aromatic volatile contents 9 .
Organic acids are the key players in maintaining pH and changes in the sensorial quality of fruit 10 . Assessment of fruit maturity depends on the sugar to acid ratio along with the quality of a certain cultivar 11 . The major sugars and organic acids in watermelon are fructose and malic acid 12 . Transport of sucrose to the fruit occurs via the phloem as part of the sugar metabolic pathway, and then neutral invertase converts sucrose into glucose, and fructose or sucrose synthase converts sucrose into UDPglucose and fructose 13,14 . Two enzymes, namely, hexokinase and fructokinase, are involved in the phosphorylation of fructose and glucose into fructose 6 phosphate (F6P) and glucose 6 phosphate (G6P), respectively. Phosphoglucoisomerase is the enzyme responsible for catalyzing the interconversion of F6P and G6P in a reversible reaction. Special transporters confined to the vacuolar membrane can facilitate the transport of sucrose to the vacuole, while acid invertase can convert it into fructose and glucose 13 . These changes are mainly controlled by genetic factors, including variations in individual gene expression patterns and metabolic gene networks 15,16 . Changes that genetically regulate variations among the organic acids in fruits during their development have been reported for many fruits, including peach, lemon, pineapple, apple, and strawberry.
Advances in omics approaches along with quantitative biology offer several ways to identify gene networks and their regulatory mechanisms in living systems, and one such promising approach for identifying coexpressed gene networks using mRNA-Seq data is the weighted gene coexpression network analysis (WGCNA) 17 . WGCNA is useful for identifying modules/networks of coexpressed genes, correlating these modules with phenotypic traits and detecting key genes within the networks. This method has been previously applied to discover the coexpressed genes responsible for fruit acidity 18 , taste, and aromatic qualities of apricot 19 and flavanol biosynthesis compounds in apple 20 . However, no such comprehensive studies are available to explain the gene networks and molecular regulatory mechanisms underlying organic acid and sugar regulation in watermelon fruit. Here, we used RNA sequencing data to uncover the regulatory networks and mechanisms controlling sugars and organic acids in sour watermelon (SrW) and sweet watermelon (SwtW) at 10, 18, 26, and 34 days after pollination. The differentially expressed genes (DEGs) and coexpression network analysis revealed gene networks linked to the metabolism of acids and sugars.

Plant material
For this experiment, we selected two kinds of watermelon materials, namely, sweet watermelon (203Z) and a sour watermelon (SrW), to evaluate the variations in sugar and organic acid metabolic content. The details of the plant materials used in this study were explained in our previous report 21 . Briefly, the pure inbred line watermelon cultivar "203Z" accumulates high total soluble sugar contents and low total organic acid contents at maturity. "SrW", developed by crossing the inbred cultivar "203Z" and the wild subspecies "PI271769", accumulates low total soluble sugar contents and high total organic acid contents. The plant materials were cultivated at the Horticultural Research Base in Henan, China, in spring 2018. Plants were sown following a randomized complete block design with 10 plants in each row and a 50 cm distance between rows. Self-pollination was performed manually for the female flowers, and to record the date of pollination, watermelon plants were tagged. Moreover, samples from watermelon flesh were collected at 4 different developmental stages (10,18,26, and 34 DAP), immediately transferred to liquid nitrogen, and finally stored at −80°C for further sugar and organic acid determinations and extraction of RNA for validation of RNA-Seq data.

Measurement of total soluble solids and pH in watermelon flesh
After homogenizing collected watermelon flesh samples, TSS (Brix %) and pH were measured. For TSS measurements, a refractometer was used (HC-112ATC, Shanghai LICHENKEYI, China). pH was recorded by using a pH meter (PHB-4, Shanghai LICHENKEYI, China) at 10, 18, 26, and 34 days after pollination 22 .

Quantification of organic acids and soluble sugars in watermelon fruit flesh
Watermelon samples were cut diagonally into two parts. Flesh samples were collected from both portions by carefully avoiding seeds in the middle portion. After taking flesh samples from both portions, the two parts were pooled together. Pooling was performed for the flesh material of two parts of the same sample. The samples were instantly frozen and then used to determine the organic acid (citric acid, oxalic acid, and malic acid) and soluble sugar (sucrose, fructose, and glucose) contents as described in previous reports 22,23 . Flesh samples were taken at 4 developmental stages (10,18,26, and 34 DAP) from 3 biological replicates.

Construction and sequencing of the RNA-sequencing library
The RNA-Seq library was prepared by a previously described method 24 , and a 1% agarose gel was used to check for contamination and degradation of RNA. The purity of RNA was analyzed by using a Nano Photometer ® spectrophotometer (IMPLEN, CA, USA). Estimation of RNA concentration was performed using a Qubit ® RNA Assay Kit and a Qubit ® 2.0 Fluorometer (Life Technologies, CA, USA). The integrity of RNA was checked by using an Agilent Nano 6000 assay kit (Santa Clara, California, United States).
Mapping to watermelon genome V2 and quantification of gene expression Clean reads were aligned to the watermelon reference genome (http://cucurbitgenomics.org/organism/21) by using TopHat v2.0.12. Read counting in features (genes, in this case) was performed using HTSeq v0.6.1 25 . Gene lengths and read counts that were mapped to genes were used to calculate FPKM values 26 .

Analysis of differential expression
The differentially expressed genes (DEGs) between SrW and SwtW were identified with the DESeq R package 27 , and Benjamini-Hochberg-adjusted P-values < 0.05 were considered statistically significant 28,29 .
Enrichment analysis of DEGs and WGCNA for identifying correlated gene networks KOBAS software was used to test for enrichment of DEGs in KEGG pathways 30 . WGCNA was performed in R using default parameters to simplify genes into coexpressed modules 17,31 . FPKM values were normalized, and an adjacency matrix was constructed. The phenotype data were imported into the WGCNA package, and correlation-based associations between phenotypes and gene modules were calculated using the default settings. The WGCNA package was used to convert the adjacency matrix into a topological overlap matrix (TOM). After constructing a network, the transcripts with identical patterns of expression were grouped into one module, and eigengenes were also calculated for these modules. The genes from each module were exported using the default parameters for Cytoscape export.

Validation of intramodular candidates through RT-qPCR analysis
Gene expression analysis using RT-qPCR was carried out using three independent biological replicates for each sample 32 . Extraction of total RNA was performed by using a TIANGEN kit (Beijing, China). cDNA was synthesized by using a Takara Tokyo, Japan) kit with reverse transcriptase. The target candidate genes were selected from the gene modules based on their intramodular gene significance and the annotation information from the watermelon reference genome database (http:// cucurbitgenomics.org/organism/21). Primer design and RT-qPCR were carried out as described previously 12 . Actin (Cla016178) was used as a reference gene 33 .

Results
Variations among soluble sugars during watermelon fruit flesh development TSS was measured at 10, 18, 26, and 34 DAP. Significant variation was observed at different developmental stages. TSS increased from 4.13 and 5.60 Brix% at 10 DAP to 12 and 10.5 Brix% at 34 DAP in SwtW and SrW, respectively (Fig. 1a).
A significant difference was observed in soluble sugar content, i.e., glucose, sucrose, and fructose, over the course of watermelon fruit development. The glucose content was measured as 11.6, 37.01, 30 (Fig. 1d). These results confirmed that as the fruit moves towards maturity, fructose and glucose are broken down to generate the energy necessary for various metabolic processes and fruit ripening 11 . The sucrose content was very low at 10 DAP but increased significantly at 18 DAP and 26 DAP (9.67, 13.28, and 5.2, 10.60 mg/g FW) and 34 DAP (25.83 and 22.39 mg/g FW) in SwtW and SrW (Fig. 1e).

Variations among organic acids throughout the development of watermelon fruit flesh
The difference in pH was recorded at different developmental stages. For example, SwtW had a pH of 4.92, 4.98, 4.76, and 5.33 at 10, 18, 26, and 34 DAP, respectively, whereas, in SrW, a pH of 4.53, 4.61, 4.38, and 4.48 was observed, respectively, at the 4 stages (Fig. 1b).
In the present study, malic acid was observed as the major acid in both materials. An increase in the malic acid content was observed from 10 DAP to 26 DAP in SwtW and SrW ( Fig. 1f). At maturity, the malic acid content decreased but was still higher in SrW than in SwtW. The concentration of citric acid was extremely high at 26 DAP ( Fig. 1g) but decreased at maturity. In the case of oxalic acid, the concentration was higher at 10 DAP in SwtW and SrW but decreased near maturity, i.e., to 0.88 mg/g FW in SwtW and 1.35 mg/g FW in SrW (Fig. 1h). Transcriptome analysis and identification of DEGs in watermelon fruit For transcriptome profiling, samples were collected at 4 stages of watermelon fruit development (10,18,26, and 34 DAP) from 3 replicates of both SwtW and SrW to construct 24 cDNA libraries. In total, 191.86 Gb of clean data was filtered, and 6 Gb of clean data for each individual sample was obtained with more than a 91% Q30 base percentage. When sequencing data were aligned to the reference genome, more than 90.1% of clean reads were mapped. Less than 3% of the reads were not mapped (Table 1).
We compared the transcriptome profiles of SwtW and SrW at 4 fruit developmental stages to identify the DEGs during watermelon fruit development ( Supplementary  Fig. 1). The DEGs were filtered according to an expression level | log2 (fold-change) | > 1 and an adjusted P-value < 0.05 in each pairwise comparison. The downregulated DEGs were more abundant than the upregulated DEGs at each developmental stage, but the upregulated DEGs were more abundant than the downregulated DEGs when the developmental stages were compared between the two materials ( Table 2).
KEGG analysis of DEGs identified from pairwise comparisons between developmental stages during watermelon fruit development provided additional information about the enriched biological pathways, including secondary metabolite biosynthesis, starch and sucrose regulatory metabolism, pentose and glucoronate interconversions, ubiquitin-mediated proteolysis, terpenoid backbone biosynthesis, glycolipid metabolism, flavonoid biosynthesis, fatty acid biosynthesis, elongation and degradation, and ascorbate and aldarate metabolism ( Supplementary Fig. 2). Based on the statistical significance criterion for multiple testing correction (adjusted P-value), phenylpropanoid biosynthesis, glycerophospholipid metabolism, MAPK-signaling pathway, secondary metabolite biosynthesis, metabolism of sucrose and starch, metabolic pathways, and plantpathogen interaction pathways were significantly enriched (q-value ≤ 0.001) and displayed a constant increase in expression throughout watermelon fruit development ( Supplementary Fig. 3).

Identification of coexpressed gene networks and key candidates
From the transcriptome data entered into the WGCNA module with FPKM values, 11 distinct gene modules were identified based on the coexpression patterns of individual genes. These gene modules are labeled in distinct colors and presented as a clustergram and network heatmap (Fig. 2a). The TSS, pH, and glucose, sucrose, fructose, malic, citric, and oxalic acid contents at each developmental stage were used as phenotypic data for the analysis of module-trait correlations. A sample dendrogram and trait heatmap were constructed to illustrate the expression of each phenotypic parameter at different developmental stages (Fig. 2b).
Of the 11 coexpressed gene networks, 3 showed significant correlations with sugars and organic acids. The blue module contained 1301 genes and showed a significant association with sucrose with a correlation coefficient (r 2 ) of 0.92. The yellow module with 410 genes was positively correlated with fructose (r 2 = 0.79) and glucose (r 2 = 0.83). The brown module consisting of 741 genes showed significantly high correlations with malic acid (r 2 = 0.75) and citric acid (r 2 = 0.85). A detailed description of the gene module and trait correlations is presented in Fig. 3a, b. The genes from each of these 3 modules were selected, and their edges and nodes were calculated with the WGCNA R package for gene-network visualizations. To further search for candidate genes with major contributions within the gene networks, annotation information of all these genes was extracted from the watermelon reference genome annotation database. This integrative approach of combining intramodular hub genes consistent with DEGs and comparisons with annotation information was followed by RT-qPCR validation of 23 selected genes for precise identification of key candidates. These genes comprised 15 genes involved in the metabolism of sugar, 6 genes linked to malic acid metabolism, and 2 genes linked to citrate metabolism. The results obtained from the RT-qPCR and RNA-Seq data were consistent with each other ( Supplementary  Fig. 4). Finally, 3 genes (Cla97C01G000640, Cla97C05G087120 and Cla97C01G018840) from the yellow module, correlated with glucose (r 2 = 0.83), were identified as sugar transporters (SWEET, EDR6 and STP), and 1 gene (Cla97C03G064990) correlated with sucrose (r 2 = 0.92) was identified as a sucrose synthase (SS) from the information available for other crops. Similarly, 3 genes (Cla97C07G128420, Cla97C03G068240, and Cla97C01G008870) from the brown module, having a strong correlation with malate (r 2 = 0.75) and citrate (r 2 = 0.85), were annotated for malate and citrate transporters (ALMT7, CS, and ICDH, respectively). These putative candidates are highlighted in red in the gene networks ( Fig. 3c-e). The genes described as key candidates were selected based on (1) intramodular significance (hub genes from modules correlated with sugars and organic acids) and (2) the gene expression profiles. In addition to the 7 genes predicted as key candidates, other hub genes from coexpression modules are listed with their annotations in Table 3, Supplementary Fig. 5, and Supplementary File 1.

Genes related to soluble sugar metabolism
Some of the important genes that were identified as associated with sugar metabolism had a considerable impact on the traits measured at various stages of fruit development. Five genes linked to SS and 4 genes associated with SPS genes showed differential expression between SrW and SwtW. One of these sucrose synthase genes (Cla97C03G064990) was also present in the blue module that strongly correlated (r 2 = 0.92) with sucrose content. The expression patterns of 2 sucrose phosphate synthase genes (Cla97C06G120770 and Cla97C11G215270) were consistent with that of the sucrose content at 4 developmental stages. Five fructokinase (FK), 4 fructose bisphosphate aldolase (FBA), and 4 fructose 1-6, bisphosphate (FBP) genes were found to be differentially expressed at various stages of watermelon fruit development. The expression patterns of FK (Cla97C05G094630), FBA (Cla97C04G076620) and FBP (Cla97C08G147060) were consistent with the fructose Fig. 3 Gene networks and key candidate genes involved in sugar and organic acid regulation during watermelon fruit development as identified by WGCNA. a Hierarchical clustering presenting eleven modules having coexpressed genes. Each leaflet in the tree corresponds to an individual gene. b Module-trait associations based on Pearson correlations. The color key from green to red represents r 2 values from -1 to 1. c Gene network for the blue module, which positively correlated with the sucrose content (r 2 = 0.92, P = 0.001). d Gene network for the brown module, which positively correlated with the fructose (r 2 = 0.73, P = 0.04), malic acid (r 2 = 0.75, P = 0.03) and citric acid (r 2 = 0.85, P = 0.007) contents. e Gene network from the yellow module, which positively correlated with the fructose (r 2 = 0.79, P = 0.02) and glucose (r 2 = 0.83, P = 0.01) contents. Hub genes (key candidates) within each network are highlighted in red due to the highest weight within the module and coded for gene descriptors based on annotations content at different developmental stages. Neutral invertase (Cla97C05G088310) was downregulated. Among the 5-hexose kinase (HK) genes, Cla97C03G063270, and Cla97C11G208480 were upregulated. Of the 12 SWEET genes screened, only WMSWEET1 (Cla97C01G000640) was highly upregulated throughout the development of watermelon fruit. Expression of the sugar transporters EDR6 (Cla97C05G087120) and STP (Cla97C01G018840) showed decreasing trends at early stages, but during the last stages of fruit development, a significant increase in their expression was observed, suggesting their roles in increasing sugar transport and accumulation (Fig. 4).

Genes linked to the metabolism of organic acids
Malic acid is the key acid in watermelon and is thought to play a critical role in developing the sour taste and reducing the pH 12 . In total, 8 ALMT genes were found in the transcriptome data, but only ALMT7 (Cla97C07G128420) was consistent with the brown module and had a significant association with malate content. The expression of this gene was higher in SrW than in SwtW. One gene, namely, 'malate synthase' (Cla97C03G054690), involved in the synthesis of malate was more highly expressed in SrW than in SwtW. The transcript levels of 2 malate dehydrogenase genes (Cla97C06G125760 and Cla97C05G103110) and 1 PEPC gene (Cla97C11G223580) increased at early developmental stages but remained stable later. One gene for transmembrane malate transport (Cla97C02G049340) was also found to be upregulated. For malate degradation, 2 pyruvate kinase genes (Cla97C04G076530 and Cla97C11G218660) showed a considerable increase in expression. One citrate synthase gene (Cla97C03G068240) was found in the brown module, and 1 isocitrate dehydrogenase gene (Cla97C01G008870) was found in the transcriptome data related to citric acid synthesis. The expression of these 2 genes was consistent with the citric acid content in SrW and SwtW at 4 developmental stages (Fig. 4).

Validation of candidate genes in different watermelon accessions
For further validation of candidate genes linked to sugar and organic acid transport, the expression profiles were investigated in 11 watermelon accessions. These 11 genotypes were selected based on the genetic diversity of sugars and organic acids to investigate the patterns of expression in diverse types of watermelon. Among these genotypes, Egusi and Amarus have low TSS and pH values, the edible-seed watermelon has low TSS and high pH values, and the improved lines have high TSS and pH values. The expression profiles of 7 key candidates (genes) identified from coexpression networks were checked via RT-qPCR at 4 developmental stages, i.e., 10, 18, 26, and 34 DAP. The genes linked to acids showed higher expression in accessions with low TSS and pH values. The genes linked to sugar transport were highly expressed in accessions with high pH and TSS values, which is consistent with the low acid contents and high sugar content. Thus, our results confirm that SWEET, STP, ERD6, SS, ALMT7, CS, and ICDH are potential candidate genes linked to sugar and organic acid transport and regulation during watermelon fruit development (Fig. 5).

Discussion
The development and maturation of fruit are composite biological processes that are regulated by various environmental, hormonal, and gene regulation factors 7 . Sugar and acid regulatory pathways are vital metabolic constituents during fruit development and maturation. In this study, by using WGCNA for the first time in watermelon and integrating transcriptome and metabolome data, we identified complex genetic factors underlying variations in sugar and organic acid pathways in watermelon.

Genetic factors underlying metabolic variations at various developmental stages
Our data showed that the fructose and sucrose contents were incredibly low in SwtW and SrW at the initial fruit development stages. Toward maturity, both of these compounds accumulated in large amounts, and a significant difference appeared between SwtW and SrW, especially for the sucrose content, which showed a substantial increase near maturity. Altogether, the fructose content was higher than the sucrose and glucose contents, which is similar to the findings of some earlier studies 12,34 . Likewise, our data showed that sugar metabolism in fruits mainly depends on developmental processes. The enzymes that regulate these extreme changes in sugar contents at different developmental stages, i.e., sucrose phosphate synthase (SPS), neutral invertase, fructose bisphosphate aldolase, sucrose synthase (SS), fructokinase (FK), and hexokinase (HK), exhibited varied expression at various developmental stages, which is consistent with the findings of the previous studies 35 . This consistency in developmental stage-dependent variation in genetic factors is also caused by variations in the sugar content itself. Sugar molecules, similar to hormones, act as signal regulators to alter gene expression levels. These signals are induced by the interaction of sugar molecules and sensor proteins 36 .
Organic acids usually exhibit high accumulation during the early developmental stages of fruit and tend to decrease near maturity; a similar trend was observed in this study. This decrease in organic acids at later stages results from the enhancement of sugar synthesis and secondary metabolic pathways towards ripening 37 . In this study, early accumulations of malate and citrate were observed, and as the maturity process continued, the malate and citrate contents decreased; these acids showed similar trends, but their points of degradation occurred early during development. These results were consistent with those of studies conducted on different fruits, including apple, peach, plum, and loquat 38,39 . Physiological changes during watermelon fruit development lead to substantial changes in fruit size, weight, flesh firmness, TSS, SSC, and acid content, thus affecting the overall quality of watermelon fruit. Therefore, it is vital to identify the key genetic factors, including gene networks and major contributors, controlling variations in these compounds, especially sugars and organic acids.

Gene networks regulating sugar and organic acid metabolism
Complex biochemical pathways such as sugar synthesis and organic acid synthesis are regulated by multigene responses and cannot be explained by individual genes. Genes in these metabolic pathways are often regulated via coordinated expression, and hence, correlation-based models are often applied to identify gene networks. Previous studies have characterized several genes regulating sugars and organic acids in watermelon and other fruits. Some examples are the tonoplast sugar transporter ClTST2 in watermelon 40 , tonoplast sugar transporter PpTST1 in peach 41 , and tonoplast sugar transporter CsPH8 and citrus transcription factor CitVHA-c4 in citrus fruits 42 . Similarly, ALMT-family genes have often been reported to be involved in organic acid regulation in fruits and other plant parts [43][44][45] . A notable finding in all these reports is that these genes only slightly contribute to the overall variation in sugar and/or organic acid content because these complex metabolic traits are regulated by gene networks. We identified 11 gene modules/networks with WGCNA (a weighted-gene correlation module), among which 3 were highly correlated with sugars and organic acids. Candidate genes for sugar metabolism in watermelon during the fruit maturation process mainly determine sweetness at the harvesting stage. Sugar molecules have to move across membranes during their transport into the chloroplast 46 , the vacuoles 47 , and the Golgi apparatus 48 . Multiple transporter genes are involved in this movement of sugars across membranes; for example, SWEET proteins form a distinct transporter family that consists of 17 members in Arabidopsis and 21 in rice. These members can transport sugars and are involved in the loading and storage of sugar molecules 49,50 . Similarly, proteins from the MST subfamily ERD6 or ESL1 are possibly involved in energy-independent sugar efflux from the vacuole 51,52 . In watermelon, however, there is little information available for the gene networks regulating sugar transport, the candidate genes linked to these networks, and their transcriptional regulation. We identified gene networks based on the coexpression of genes at the transcriptional level and their correlations with sugar content variation at various developmental stages. Within these networks, we found 3 genes (Cla97C01G000640, Cla97C05G087120, and Cla97C01G018840) annotated as sugar transporters, namely, SWEET, EDR6, and STP, and 1 gene (Cla97C03G064990) for sucrose synthesis. These genes were identified as hub genes within the modules, and the expression profiles in diverse watermelon accessions were consistent with the sugar content at various developmental stages. Thus, our confidence in these genes being the true variants for sugar transport is extremely high. None of these genes have been reported or studied in previous watermelon research and, therefore, should be investigated thoroughly in future molecular biology studies to fully understand the genetic regulators of sugar metabolism in watermelon.
Phosphoenolpyruvate carboxylase, aconitase (ACO), citrate synthase (CS), malate dehydrogenase (MDH), isocitrate dehydrogenase (ICDH), and malic enzyme (ME) are enzymes linked to organic acid metabolism and play a vital role in the biosynthesis and degradation of organic acids in fruit. Malic acid, a major organic acid in watermelon, is manufactured during the TCA cycle during different fruit growth and developmental stages. This acid is degraded during gluconeogenesis, and the TCA cycle takes place in the cytosol and mitochondrion. The vacuole is the major storage organ for malate.
We found that the brown module from WGCNA correlated significantly (r 2 = 0.75) with the malate content. From the network built from the top genes in the brown module, Cla97C07G128420 was identified as a hub gene, which encodes an ALMT-family protein. This gene was also identified as a differentially expressed gene between sweet and sour materials, and its relative quantitative expression also showed it to be consistently regulated with the malate content at various developmental stages.
ALMT genes have also been reported in other crops to regulate organic acid contents, for example, AtALMT9 and AtALMT6 in Arabidopsis thaliana and four AtALMT9 homologs in grape berries 53 . Similarly, in apple, genes linked to fruit acidity and/or pH include two ALMT-like genes, and ALMT II has already been reported 54,55 . According to these literature data and our findings, it is highly likely that Cla97C07G128420 is a key candidate in the malate pathway's gene network that contributes to the maximum trait variation. Similarly, Cla97C01G008870 annotated as ICDH and Cla97C03G068240 as citrate synthase from the brown module exhibited relative expression complementary to the citric acid content at all 4 developmental stages under investigation ( Supplementary Fig. 4). These results suggest that these 3 genes are key players in the TCA cycle and pH regulation in watermelon. Fully characterizing these genes using molecular biology techniques would deepen our understanding of organic acid metabolism in watermelon.
Potential genes in sugar and organic acid pathways not detected as key candidates As mentioned earlier, we predicted key candidate genes according to a) hub genes from modules correlated significantly with fructose, glucose, malic acid, and citric acid and b) the expression profiles of selected genes to narrow down the search. However, several other genes encoding proteins for sugar and organic acid pathways also showed differential expression between diverse watermelon types. For example, Cla97C05G094630, which is associated with fructose kinase protein, showed downregulated expression, which would suggest that fructose is less utilized during the metabolic process 56 . Likewise, higher accumulation of fructose was observed at the maturity stage in our material, which suggests that Cla97C05G094630 does play a vital role in fructose regulation; however, this gene was not identified as a potential candidate based on intramodular gene significance. Similarly, we found two hexokinase-related genes (Cla97C03G063270 and Cla97C11G208480) showing upregulated expression at maturity. Hexokinases are commonly known to be involved in the glycolysis process; however, these genes were not among the high-confidence genes based on our selected models. The accumulation of sucrose is primarily linked with the activity of sucrose-metabolizing enzymes, including SPS, and sucrose synthase 57 . Two SPS-related genes (Cla97C11G215270 and Cla97C06G120770) and one sucrose synthase gene (Cla97C03G064990) were found in the blue module, which was positively correlated with sugar content. The relative expression of these genes also complemented the sugar content at various developmental stages, but their intramodular significance was not exceedingly high. Cla97C05G103110, identified in the brown module and correlated with malate content, is associated with the NAD-MDH enzyme based on available data from other crops. This gene also showed expression consistent with the malate content during various developmental stages in our study. Other important genes, such as Cla97C03G066480 (malic enzyme), Cla97C11G223580 (PEPCK), Cla97C04G076530, and Cla97C11G218660 (both annotated for pyruvate kinases), also showed expression trends consistent with changes in relevant traits.
These results indicate that complex metabolic pathways such as those for sugars and organic acids regulate the coordinated expression of gene networks and that certain genes within the networks make key contributions to trait variation. One limitation of our study should also be considered: the statistical power of our analysis was low due to the small sample size. Nevertheless, the expression profiling of key candidates in genetically diverse accessions supported our findings. In future studies, larger and more diverse samples should be investigated to further demonstrate the coexpression patterns of genes/gene networks. Additionally, future studies should focus on how certain changes in the gene networks can affect other pathways or traits, and to do so, high-throughput phenotyping of agronomic and metabolic traits should be incorporated into the analysis.

Conclusions
In the present study, we used transcriptome profiles to investigate gene networks (based on coexpression patterns) controlling the regulation of sugars and organic acids in watermelon. We identified 3 coexpression modules/gene networks that were highly correlated with variations in the sugar and organic acid contents at various developmental stages. Within these networks, we identified 7 key candidate genes that were weighted as module hub genes, and their quantitative expression was correlated with phenotypic variation. The expression profiles of these 7 genes in genetically diverse watermelon accessions further supported our findings. In total, we identified 4 candidate genes for sugar content and 3 for organic acid content. The multiple and stringent screening steps used in our study increased the probability of and confidence in these genes being true candidates in the networks controlling sugar and organic acid pathways. These genes are newly identified, and no previous reports about their role or functions in watermelon are available.

Core ideas
• The integration of metabolic phenotypes and gene expression profiles can reveal gene networks and key candidate genes regulating metabolic variations.
• This study reports 4 new high-confidence candidate genes (Cla97C01G000640, Cla97C05G087120, Cla97C01G018840 and Cla97C03G064990) controlling sugar transport and regulation in watermelon.
• This study also reports 3 new high-confidence candidate genes (Cla97C07G128420, Cla97C03G068240, and Cla97C01G008870) for malic acid and citric acid regulation in watermelon.