The RNA-seq data for COAD and LIHC were downloaded from TCGA portal (18). The data type derived from TCGA was used only for STAR-Counts. We obtained 437 COAD and 424 LIHC RNA-seq datasets. To identify metabolic alterations during the early stages, stage I cancer data were selected by comparison with the metadata derived from TCGA. Finally, we obtained 39 normal and 62 tumor samples from COAD, and 50 normal and 171 tumor samples from LIHC.

To ensure data quality, we filtered the STAR counts by removing those with average counts of less than one in all patients. We then applied DESeq2 in Bioconductor (19) to normalize the filtered count data and extract differentially expressed genes (DEGs) from normal and tumor tissues with an adjusted p-value cutoff of 0.01. To visualize the DEGs, we used a cutoff of |log2foldchange (log₂FC)| > 0.58 and converted any genes with p-adjust value (p_adj) or Log₂FC as NA to “1” to prevent undetectable error. The DEGs were displayed using Enhanced Volcano in Bioconductor (20), where the gray dots represented “non-DEGs,” red dots represented “log_¬2FC > 0.58 and p_adj <0.01,” and blue dots represented “log₂FC < -0.58 and p_adj <0.01”.

Each gene in the normal and tumor tissues in COAD and LIHC contained numerous dimensions. To visualize the genes, dimensionality reduction was performed using principal component analysis (PCA) and the results were visualized using ggplot2 in R (21). The PCA plot visualizes PC1 on the x-axis and PC2 on the y-axis, and the normal and tumor groups are represented by ellipses.

To comprehensively understand the functions of the DEGs, we conducted a Gene Ontology (GO) enrichment analysis using ClusterProfiler in Bioconductor (22). Specifically, we used a p-adjusted value cutoff of 0.01 for genes with a log₂FC > 0.58 and log₂FC < -0.58 to indicate upregulated and downregulated genes, respectively. To confirm the metabolic process alterations in the early stages of tumorigenesis, we focused only on biological process (BP) terms that indicate cellular or physiological effects. The results of the GO enrichment analysis are displayed as a heatmap with -log10 p-values, where the upregulated gene set is depicted in red, and the downregulated gene set is depicted in blue. After conducting the GO analysis, we visualized the results using a heatmap. A heat map was generated using the pheatmap function in Bioconductor, which showed the expression levels of the identified genes (23).

To further investigate the metabolic processes involved in COAD and LIHC, we utilized the Gene Set Enrichment Analysis (GSEA) tool provided by ClusterProfiler in Bioconductor (24). The analysis was conducted using a p-value cutoff of 0.05, and only BP (biological process) gene set terms were considered to compare metabolic processes in both cancers. The GSEA results are presented using an enrichment plot in Bioconductor (25) and include the normalized enrichment score (NES) and corresponding p-value.

In this study, we performed constraint-based simulations using two genome-scale metabolic models (GSMs) to elucidate the functional role of HMGCS2 in cancer metabolism. Specifically, we utilized the colon cancer model (26) and the iHepatocytes2322 curated liver model (27) and conducted simulations using the COBRA Toolbox v.3.0[28] and the method of minimization of metabolic adjustment (28). We generated HMGCS2 knock-out colon models by limiting the lower bounds of the HMGCS2-related reactions (HMR1437, HMR4604, and HMR1573) to nine, while the HMGCS2-overexpressed colon models had upper bounds of 4000 for these three reactions. Similarly, HMGCS2 knock-out liver models were derived from iHepatocytes2322 by limiting the lower bounds of HMGCS2-related five reactions (HMR1437, HMR4604, HMR1573, HMR0027, and HMR0030) to 0, while HMGCS2-overexpressed liver models had a lower bound of 2000 and an upper bound of 4000 for these five reactions.

To investigate the functional role of HMGCS2 in cancer metabolism, we observed changes in reaction flux by genetically altering HMGCS2. Specifically, we defined reactions whose flux decreased in HMGCS2 knock-out models and increased in HMGCS2 overexpression models as “flux decreasing” reactions, while reactions whose flux increased in HMGCS2 knock-out models and decreased in HMGCS2 overexpression models were defined as “flux increasing” reactions. We then counted the number of flux-increasing and decreasing reactions per subsystem and categorized these numbers by the total number of reactions in each subsystem to summarize flux changes.

Next, we analyzed the effects of gene perturbation of HMGCS2 in glycolysis and lipid metabolism in colon and liver models. Specifically, we calculated flux changes by subtracting the fluxes of the original models from those of the perturbation models and considered flux changes higher than 10% of the original flux with positive and negative signs as “up-regulated” and “down-regulated,” respectively. Reactions whose changes were neither up- nor down-regulated were assigned as “no change,” while reactions that were unidentified in the model were indicated as “unidentified”.

Colon cancer (Caco-2) cells, derived from human colorectal adenocarcinoma, were procured from ATCC and maintained in Dulbecco’s Modified Eagle’s Medium supplemented with 10% fetal bovine serum and 1% penicillin-streptomycin-amphotericin B at 37°C with 5% CO₂. To target HMGCS2 [NM_001166107.1 and NM_005518.3], siRNA sequences were purchased from Bioneer (Korea), and Lipofectamine RNAiMAX (Thermo Fisher Scientific, Inc., MA, USA) was used to transfect the siRNA according to the manufacturer’s instructions.

To measure the Oxygen Consumption Rate (OCR) and Extracellular Acidification Rate (ECAR) of Caco-2 monolayers, we employed a Seahorse XFp Extracellular Flux Analyzer (Agilent Technologies, Santa Clara, CA, USA). The Seahorse XFp Sensor Cartridge was pre-hydrated with XFp Callibrant solution one day prior to the test and incubated overnight at 37°C in a CO₂-free incubator to eliminate CO₂, which could interfere with pH-sensitive measurements. Subsequently, Caco-2 cells were seeded onto XFp Miniplates at a density of 2×10⁴ cells/well and allowed to settle overnight. On the day of the assay, the complete growth medium was replaced with 180 ul/well of XF assay medium, which was maintained at 37°C in a non-CO₂ incubator for 1 h to allow pre-equilibration with the XF assay medium. We then analyzed the mitochondrial function of the cells by sequentially injecting oligomycin (1 µM), carbonyl cyanide-4 (trifluoromethoxy) phenylhydrazone (FCCP, 0.5 µM), and a mix of rotenone and antimycin A. Finally, OCR and ECAR values were normalized using cellular protein content.

The present study aimed to identify common genetic foundations and related signaling pathways in GI malignancies. We extensively analyzed the publicly accessible TCGA database, focusing on the COAD and LIHC datasets comprising 437 and 424 samples, respectively. To investigate the metabolic changes in early tumorigenesis, we used only Stage I cancer data for further analysis, resulting in 39 normal samples and 62 tumor samples for COAD, and 50 normal samples and 171 tumor samples for LIHC.

As is demonstrated in Supplementary Figure 1A , the PCA plot clearly displays distinct elliptical clusters that effectively separated the normal and tumor samples. This supports the notion that the expression profiles of GI systems change substantially due to tumorigenesis. Using a list of DEGs, we generated volcano plots ( Figure 1A ) to identify significant differences in gene expression profiles between normal and cancer tissues. We found 7837 and 8767 up-regulated genes and 7232 and 3642 down-regulated genes in the colon and liver tissues, respectively. Tables 1 , 2 show the top ten upregulated and downregulated DEGs in both COAD and LIHC tissues based on p-values. In COAD, ETS variant transcription factor 4 (ETV4), keratin 80 (KRT80), and forkhead box Q1 (FOXQ1) were the top three upregulated genes, whereas estrophin 4 (BEST4), glycolipid transfer protein (GLTP), and carbonic anhydrase 7 (CA7) were the top three downregulated genes. Similarly, in LIHC, plasmalemma vesicle-associated protein (PLVAP), collagen type XV alpha 1 chain (COL15A1), and gamma-aminobutyric acid type A receptor subunit delta (GABRD) were the top three upregulated genes, whereas ADAM metallopeptidase with thrombospondin type 1 motif 13 (ADAMTS13), oncoprotein induced transcript 3 (OIT3), and stabilin 2 (STAB2) were the top three downregulated genes.

To gain further insight into the metabolic pathways that were enriched during the early stages of tumorigenesis, we conducted a pathway enrichment analysis using GSEA. As shown in Figure 1B and Supplemental Figure 1B , the heatmap displays the enriched pathways in cancer and normal tissues. The analysis revealed that the genes differentially expressed during early tumorigenesis are associated with various aspects of cell cycle regulation. Notably, genes involved in “DNA replication,” “mitotic nuclear division,” and “cell cycle G2/M phase transition” were found to be positively enriched in both COAD and LIHC. Furthermore, the results indicate that genes related to lipid metabolism were significantly enriched in COAD and LIHC. Specifically, “fatty acid beta-oxidation (FAO)” and “cellular lipid catabolic process” were found to be negatively associated with early tumorigenesis in both cancer types ( Figures 1B–D ). These findings suggested that dysregulated lipid metabolism is crucial in the development and progression of COAD and LIHC.

To assess glucose metabolism in both COAD and LIHC groups, we compared the gene expression of key irreversible enzymes involved in regulating glycolysis and gluconeogenesis ( Figure 1E ). The major rate-limiting enzymes in glycolysis, including phosphofructokinase homologs (PFKP and PFKM) and pyruvate kinase (PKM), which were significantly increased in both COAD and LIHC. Conversely, the levels of key enzymes related to gluconeogenesis, such as pyruvate kinase (PC), phosphoenolpyruvate carboxykinase (PCK1 and PCK2), and glucose-6-phosphatase (G6PC1 and G6PC2), were significantly decreased. These findings were consistent with the expected alterations in glucose metabolism in COAD and LIHC, commonly known as the Warburg effect (29), suggesting a shift towards increased glucose uptake and utilization through glycolysis in these malignancies.

Mitochondria are key organelles in cellular energy metabolism, as they serve as the primary sites for oxidative phosphorylation (OXPHOS) and FAO, and for ATP production (30). When analyzing the DEGs in COAD and LIHC, we identified a specific subset of 426 and 325 genes, respectively, that were significantly linked to mitochondrial function and structure (31) ( Figure 2A ). Notably, among the mitochondrial genes identified, 164 were common DEGs between the two cancers ( Figure 2B ).

As shown in Figure 2C and Supplementary Figure 2 , our results demonstrate the enrichment of mitochondrial genes based on the DEGs identified between cancer and normal tissue samples. Interestingly, we observed an upregulation in genes involved in “mitochondrial protein import” and “mitochondrial respiratory complex assembly,” which are critical components of mitochondrial biogenesis and energy generation (32, 33), in both COAD and LIHC. Conversely, we noted a downregulation of genes related to “FAO” and “lipid catabolic process.” Our findings suggest a potential shift in the metabolic profile of GI cancers towards an increased reliance on mitochondrial biogenesis and a decreased dependence on lipid metabolism.

The heat map displayed in Figure 2D shows the common DEGs involved in mitochondria-related metabolism in both COAD and LIHC. Our analysis revealed a significant increase in the expression of genes associated with fatty acid synthesis, whereas most genes related to FAO were downregulated. Furthermore, we observed a decrease in several genes involved in tryptophan metabolism, including kynurenine 3-monooxygenase (KMO) (34) and monoamine oxidase A (MAOA) (35). Additionally, we noted a decrease in the expression of the succinate dehydrogenase complex subunit D (SDHD) gene, which encodes a subunit of the mitochondrial enzyme responsible for succinate oxidation and is a well-known tumor suppressor (36). These results provide important insights into the altered metabolic pathways in GI cancers, which may contribute to their development and progression.

To identify crucial candidates that regulate energy metabolism in GI malignancies, we conducted a correlation network analysis using the GeneBridge toolkit (37). This newly developed bioinformatics tool allows the imputation of gene functions and biological connectivity using large-scale multispecies expression datasets (37). The analysis revealed that 285 genes in COAD and 2399 genes, including 3-Hydroxy-3-Methylglutaryl-CoA Synthase 2 (HMGCS2), were associated with “fatty acid oxidation” (GO:0006635) ( Figure 3A ). Among these genes, 25 genes including 3-Hydroxy-3-Methylglutaryl-CoA Synthase 2 (HMGCS2) are common mitochondrial genes between COAD and LIHC. To identify crucial mitochondrial genes associated with GI malignancies, we calculated the hazard ratio (HR) for each gene’s related all-cause mortality in COAD and LIHC. Figure 3B displays the HR of common mitochondrial genes, with HMGCS2 being one of the most highly expressed HR genes in both cancers. Patients with low HMGCS2 expression had higher HR than those with high HMGCS2 expression in both malignancies.

Then, we performed survival analyses of cancer patients based on the expressions of DEGs that are commonly observed in COAD and LIHC using the GEPIA tool (38). By analyzing common DEGs, we identified a set of 25 genes that were particularly linked to FAO. Moreover, our investigation revealed 6 genes that have a noteworthy impact on the survival of patients with cancer. Of these 6 genes, HMGCS2 was the only gene that displayed a statistically significant difference in the overall survival rates of patients with both COAD and LIHC ( Figure 3C ; Supplementary Figure 3 ). Notably, the expression of HMGCS2 was found to be considerably reduced in lung cancer and rectosigmoid junction cancer, and in COAD and LIHC, compared to normal tissues ( Figure 3D ). In addition, HMGCS2 expression was also found to be significantly lower in colon and liver cancer, as shown in Supplementary Figure 4 , where we analyzed public cancer datasets for colon and liver cancer. These results indicated that HMGCS2 may play a critical role in the pathogenesis of GI malignancies.

To gain further insight into the metabolic functions of HMGCS2 in GI malignancies, we conducted genome-scale metabolic simulations using the COAD and LIHC models. In the COAD model, the suppression of HMGCS2 led to a significant increase in the fluxes of over half of the reactions in the fatty acid synthesis subsystems (i.e., fatty acid biosynthesis and elongation), whereas the fluxes in the fatty acid degradation subsystems (i.e., fatty acid destruction, beta-oxidation, and mitochondrial carnitine shuttle) were significantly reduced ( Figures 4A, B ; Supplementary Figure 2 ). Furthermore, HMGCS2 inhibition resulted in a notable upregulation in the flux of glycolysis subsystems and downregulation in the flux of oxidative phosphorylation. Remarkably, the suppression of HMGCS2 resulted in similar changes in metabolic flux in a normal liver tissue model ( Supplementary Figure 5 ). However, in the LIHC model, suppression of HMGCS2 did not cause significant changes in metabolic flux prediction.

To further investigate the role of HMGCS2 in energy metabolism in cancer cells, we measured oxidative phosphorylation and glycolysis using a Seahorse extracellular flux analyzer (39). Our investigation focused on human Caco-2 cells and aimed to explore the effects of HMGCS2 inhibition on these metabolic pathways. HMGCS2 knockdown resulted in a discernible decrease in the OCR of Caco-2 cells, suggesting decreased oxidative phosphorylation ( Figure 4C ). We also noticed a corresponding increase in the ECAR in these cells, indicating enhanced glycolysis ( Figure 4D ; Supplementary Figure 6 ). These results support the notion that the inhibitory effects of HGMCS2 alter the metabolic flux, which is in line with the predictions made by our model.