Introduction

Primary liver cancer is among the most common malignancies of the digestive system, ranking fifth in global incidence and third in mortality rates.1 By 2025, it is projected that over one million individuals will receive a diagnosis of liver cancer annually.2 HCC accounts for approximately 80% of all primary liver cancers, representing the predominant subtype.3 The primary risk factors implicated in HCC etiology include Hepatitis B Virus (HBV), Hepatitis C Virus (HCV) infection, excessive alcohol consumption, exposure to aflatoxin B1, and nonalcoholic fatty liver disease (NAFLD).4 Due to its hidden and high metastasis, it is rather difficult to accurately diagnose and treat HCC during early occurrence.5 Current treatment modalities for HCC encompass chemotherapy, radiation therapy, liver transplantation, and surgical resection.6 Despite these interventions, the five-year survival rate for advanced HCC remains a disheartening 15%, and a concerning 50% recurrence rate is noted within the first five years post-surgical liver resection for early-stage HCC.7 Thus, there is an urgent need for the development of novel prognostic models to enhance early clinical diagnosis and to improve the survival rates of HCC patients, including their susceptibility to immunotherapy and chemotherapy.

The reprogramming of glucose metabolism has emerged as a hallmark of cancer in recent years. The Warburg effect, observed in aerobic conditions, describes the preference of proliferating tumor cells for aerobic glycolysis as their primary energy source.8 The transformation of glucose metabolism from the oxidative phosphorylation pathway to the glycolysis pathway in HCC is conducive to the rapid proliferation of tumor cells, providing a favorable microenvironment for HCC progression.9 Gluconeogenesis is actually a reverse pathway of glycolysis,10 approximately 80% of endogenous glucose is produced by liver through gluconeogenesis during fasting.11 Gluconeogenesis has a mechanism of action to inhibit glycolysis and block the progression of HCC.12 Genes pivotal to gluconeogenesis, such as phosphoenolpyruvate carboxykinase 1 (PCK1), fructose-1,6-bisphosphatase 1 (FBP1), and glucose-6-phosphatase catalytic subunit (G6PC), typically act as tumor suppressors in normal liver tissue.13 However, in HCC, the expression of these gluconeogenic genes is often downregulated, which suppresses gluconeogenesis and promotes aerobic glycolysis in tumor cells, accelerating tumor growth and leading to poorer patient outcomes.14–16 While research on glucose metabolism reprogramming has historically concentrated on glycolysis, the relationship between gluconeogenesis and tumorigenesis, though less explored, may offer novel therapeutic avenues for HCC. It is crucial to investigate sensitive biomarkers from the perspective of gluconeogenesis to potentially extend the survival rates of HCC patients.

Various components that compose tumors are termed the tumor microenvironment (TME). Within the TME, significant immune cell infiltration occurs, these immune cells are called the TIME. The TIME includes cell surface molecules, external immunological mediators, and immune cells that are innate as well as adaptive that play pivotal roles in tumor progression, recurrence, and metastasis.17–19 Research has established a correlation between glycolysis and TIME and that enhanced glycolysis can attenuate the anticancer efficacy of tumor-active T cells and facilitate tumor immune evasion.20,21 Limited research has explored the association between gluconeogenesis and TIME. Therefore, investigating the role of GRGs within the TIME is essential for a thorough understanding of the potential mechanisms by which gluconeogenesis may influence HCC progression and response to immunotherapy.

In our study, we initially conducted an extensive analysis of mRNA expression profiles and clinical data from HCC patients within the TCGA database. Employing LASSO Cox regression analysis, we successfully identified four GRGs and subsequently developed a prognostic predictive model for HCC. The stability and reliability of this model were corroborated using the ICGC cohort. To further elucidate the role of GRGs in different HCC subtypes, we classified 365 patients into four subgroups based on GRG expression levels. We then performed comprehensive functional analyses, including clinical characteristics, functional enrichment patterns, patterns of immune cell infiltration, and attributes of tumor stemness. Overall, our study not only elucidated the critical role of the gluconeogenesis pathway in HCC progression but also highlighted its impact on immune therapy response and prognosis. Through comprehensive bioinformatics analyses, our study aims to contribute insights for personalized therapeutic strategies in HCC for patients diagnosed with HCC.

Materials and MethodsData Collection

In our study, we obtained gene expression profiles (log-transformed tranpackages per million, TPM) and clinical data using the TCGA-LIHC cohort (374 hCC samples and 50 normal tissue samples) on the TCGA website (https://portal.gdc.cancer.gov/, accessed on 20 November 2023).22 After excluding samples with zero or missing survival times, 365 patients with both genomic expression data and clinical data were included for further analysis. The ICGC website (https://dcc.icgc.org/projects/LIRI-JP, accessed on 20 November 2023) was used to obtain RNA sequencing data and clinical information for an additional 240 hCC samples. The genomic and clinical data of the TCGA cohort is presented in Supplementary Tables S1 and S2. The TCGA and ICGC datasets were both accessible, adhering to their respective data acquisition and release policies.

GRGs Resource

The Gene Set Enrichment Analysis (GSEA) portal (https://www.gsea-msigdb.org/gsea/index.jsp, accessed on 30 November 2023) contains a collection that includes 200 hallmark genes associated with glycolysis and 148 genes related to both glycolysis and gluconeogenesis (Supplementary Materials 1 and 2). We identified and removed 31 overlapping genes from these two gene sets, resulting in a distinct set of 117 genes specifically related to gluconeogenesis.

Construction and Validation of GRGs Prognostic Model

In the TCGA cohort, differentially expressed genes (DEGs) between tumor and non-tumor tissues were identified using the “limma” R package, with fold change (FC) >1.2 and a false discovery rate (FDR) <0.05. Subsequently, univariate Cox analysis was conducted to select genes with prognostic relevance to gluconeogenesis. The “Venn” R package was utilized to ascertain the duplicate genes. LASSO is a regularization method for linear regression problems that can reduce model complexity, mitigating overfitting, and discern pivotal predictors. Its mechanism involves minimizing the sum of squared residuals, leading to certain regression coefficients equal to 0, thereby yielding an optimized model. In this study, we employed LASSO regression analysis facilitated by the “glmnet” package in R to discern a cohort of pivotal genes.23,24 Utilizing the expression profiles of these core genes, a risk score for each patient was calculated as follows:

N denote the number of genes in the model; exp represents the expression level of a gene; coef signifies the corresponding gene’s coefficient. Patients were stratified into high-risk and low-risk groups using the median risk score as the demarcation. To visualize gene expression patterns within the model, t-distributed stochastic neighbor embedding (t-SNE) analysis was implemented to find the distribution patterns across different groups. Furthermore, Kaplan-Meier survival curves were generated using the “survival” and “survminer” R packages to analyze the differences in overall survival (OS) between high-risk and low-risk populations. Time-dependent ROC curve analysis was executed with the “survival”, “survminer”, and “timeROC” R packages to assess the predictive value of prognostic features. Additionally, the independent prognostic ability of the model was investigated in two Cox analyses (monovariate and the multivariate). The ICGC dataset was utilized to validate the prognostic model of GRGs. The risk score for each patient was computed using the aforementioned method. A Kaplan-Meier curve was generated to assess the predictive power of the GRGs prognostic model. The ROC curve was constructed using the aforementioned method.

Building and Verifying the Nomogram

The nomogram is widely used to estimate survival probabilities in cancer patients. Calibration profiles were employed to assess predicted OS values by comparing them to the actual observed outcomes.

Profiling of Functional Enrichment

To investigate the heterogeneity in biological mechanisms related to GRGs, we employed Gene Set Variation Analysis (GSVA) and GSEA utilizing the gene set “c2.cp.kegg.v7.4.symbols.gmt” from the Molecular Signatures Database (MSigDB).25

Analysis of Tumor Microenvironment Infiltration

We utilized the R packages “CIBERSORT” and “ssGSEA” to quantify the levels of immune cell infiltration within TME.26 Estimation of immune cell infiltration in each specimen using CIBERSORT, enabling a comparative analysis between the high and low-risk groups. Additionally, we conducted a Spearman correlation analysis to probe the association between risk and immune cell scores.27 The Spearman correlation was also employed to determine the association between stem cell-like traits and risk scores.

Drug Susceptibility Analysis

The Genomics of Drug Sensitivity in Cancer (GDSC2, https://www.cancerrxgene.org/, accessed on 30 November 2023) obtained drug sensitivity data of cancer cell,28 encompassing 969 cell lines and 297 compounds. The “OncoPredict” R package was used to search susceptibility data in the GDSC2 database.23 Subsequently, a Wilcoxon test was employed to identify drugs exhibiting significant differences in susceptibility between high- and low-risk groups (P < 0.001).

Consensus Clustering of GRGs

We employed the “ConsensusClusterPlus” package to implement unsupervised consensus clustering using the k-means algorithm.29 Unsupervised consensus clustering was performed to divide the patients into distinct molecular subtypes and gene expression patterns.

Cell Culture

Human HCC cell lines HepG2 were sourced from the Cell Bank of Type Culture Collection of the Chinese Academy of Sciences (CBTCCCAS, Shanghai, China). The cells were cultured in Dulbecco’s Modified Eagle Medium (DMEM; Gibco, Cat. No. 11095092), supplemented with 10% fetal bovine serum (FBS; Gibco, Cat. No. 10099141), 100 IU/mL penicillin, and 100 μg/mL streptomycin in 5% CO2 humidified atmosphere at 37°C. The cell culture experiments were performed three times in total.

Lentivirus Transfection

Lentivirus carrying the overexpressed full-length G6PC mRNA was engineered by GenePharma (Shanghai, China). To generate stable cell lines, HepG2 cells were transfected with the designated lentivirus in the presence of 6 μg/mL polybrene. Post-transfection, 4 μg/mL puromycin was added to select positive cells for seven days.

Immunofluorescence

HepG2 Cells were cultured on glass dishes for 24 hours, then fixed with 4% paraformaldehyde for 15 minutes at room temperature before permeabilization with 0.3% Triton X-100 for an additional 15 minutes. Subsequently, the cells were blocked with 10% goat serum in PBS for 60 minutes at room temperature. Following incubation with primary antibodies and subsequent staining with Alexa-conjugated secondary antibodies (Cell Signaling Technology, USA), the cells were stained with DAPI for 5 min. Finally, the cells were photographed under a confocal microscope.

Cell Proliferation Assay

Cell proliferation was evaluated through colony formation assays and the Cell Counting Kit-8 (CCK-8). A density of 1000 transfected cells per well was used to seed 96-well plates, followed by the addition of 10 μL of CCK-8 solution to each well and then incubated at 37°C for 1 hour. The absorbance at 450 nm was then quantified using a microplate reader. To determine clonal proliferation, 500 cells per well of 6-well plates were used for cell seeding. The cells were fixed with methanol and stained with 0.5% crystal violet after growing for two weeks.

Statistical Analysis

To compare the DEGs between the tumor tissue and the surrounding tissue, we ran the Wilcoxon test, and differences between the two groups were further assessed using the chi-square test. The Kaplan-Meier curve was employed to assess the variations in OS between the various groups. Independent predictors significantly associated with OS were identified by univariate and multivariate Cox analysis. Spearman correlation analysis was applied to Examine the connection between tumor stemness, stromal score, and immunological score and the prognostic model hazard score. R software (version 4.3.2) was used for all statistical analyses, which offers various functions such as data processing, statistical analysis, and visualization. For all statistical analyses, A P-value of less than 0.05 was considered to indicate statistical significance.

ResultsIdentification of GRGs

The flow diagram of this study was present in Figure 1. Our study comprised both TCGA (374 patients) and ICGC cohort (240 patients). To identify GRGs, we procured a list of 200 hallmark glycolysis genes and 148 genes involved in both glycolysis and gluconeogenesis from the GSEA website, and subsequently eliminated 31 duplicate genes (Figure 2A). Utilizing the “limma” R package, we identified 88 DEGs in tumors of HCC and nearby non-tumor tissues. Based on univariate Cox regression analysis, 43 of these 88 DEGs showed statistically significant relationships with the prognosis of HCC patients. (Figure 2B and C). 13 genes were upregulated in non-tumor tissues relative to tumor tissues from the heatmap (Figure 2D). Additionally, a network graph was constructed to investigate the correlation between HCC survival rates and the 43 GRGs identified, as well as the interactions among these GRGs, where lines connecting the GRGs signify correlations, with thicker lines indicating stronger correlations. The use of pink and blue colors denotes positive and negative correlations, respectively (Figure 2E). Given that the gluconeogenesis pathway is typically downregulated in HCC progression, we selected nine genes that were underexpressed in HCC tissues and had a strong link to gluconeogenesis (FBP1, G6PC, PPARGC1A, GNMT, ADH1C, ALDH2, ADH1A, SLC2A2, and ADH4), for in-depth investigation. Decreased levels of expression in these genes were substantially related with OS (Supplementary Figure 1). Moreover, these nine genes in cancerous tissues were observed to be lower expression levels than normal tissues (Supplementary Figure 2).

Figure 1 Flowchart of the overall research.

Figure 2 Identification of the candidate gluconeogenesis-related genes in the TCGA cohort. (A and B) Venn diagram to identify GRGs between HCC tissues and normal tissues. (C) 43 GRGs associated with prognosis. (D) The 43 GRGs expression between HCC tissues and normal tissues. (E) Network of 43 GRGs.

Construction and Validation of 4-GRGs Prognostic Model

The expression profiles of the nine genes were analyzed using LASSO-Cox regression, and a prognosis model was established. The optimal values of λ were used to identify markers for four genes (Supplementary Figure 3). Risk score calculation formula: score = (‒0.001 × Exp G6PC) + (‒0.118 × Exp PPARGC1A) + (‒0.066 × Exp ADH4) + (‒0.072 × Exp ALDH2). The patients were divided into two cohorts depending on the median score (Figure 3A). According to a scatter plot, patients at higher risk were more likely to die sooner than patients at lower risk (Figure 3B). Consistent with this, Kaplan-Meier survival curves demonstrated that those with elevated risk had noticeably lower OS (Figure 3C). Survival prediction analysis using the prognosis model generated time-dependent receiver operating characteristic (ROC), with the area under the curve (AUC): 0.680 for one year, 0.666 for two years, 0.680 for three years (Figure 3D). Similar to the TCGA cohort’s findings, the patients were also divided into two cohorts depending on the median score (Figure 3E), high-risk patients in the ICGC cohort were at an increased risk of earlier mortality and had shorter OS (Figure 3F and G). The AUC for ROC curves was 0.779 for one year, 0.672 for two years, 0.653 for three years (Figure 3H).

Figure 3 Prognostic analysis of the four-GRGs prognostic model in the TCGA cohort and ICGC cohort. TCGA cohort (A–D), ICGC cohort (E–H). (A and E) The median value and distribution of the risk scores. (B and F) The distribution of OS status. (C and G) Kaplan-Meier curves. (D and H) ROC curves.

Independent Prognostic Value of the 4-GRGs Model

We investigated whether the risk score and clinical-pathological variables function as separate prognostic indicators for OS using univariate and multivariate Cox analysis. For the TCGA and ICGC cohorts, the univariate Cox analyses revealed a significant connection between the risk score and OS (TCGA cohort: HR = 3.464, 95% CI = 2.113–5.682, P< 0.001; ICGC cohort: HR = 3.663, 95% CI = 1.627–8.250, P = 0.002) (Figure 4A and B). Tumor stage was significantly linked to OS in both cohorts (TCGA cohort: HR = 1.680, 95% CI = 1.369–2.062, P< 0.001; ICGC cohort: HR = 2.203, 95% CI = 1.519–2.195, P< 0.001) (Figure 4A and B). Furthermore, gender in the ICGC cohort was notably associated with OS (HR = 0.502, 95% CI = 0.268–0.940, P = 0.031) (Figure 4B). Upon adjusting for other confounding factors, the multivariate Cox analysis confirmed that tumor stage remained an independent predictive factor for OS (TCGA cohort: HR = 2.767, 95% CI = 1.687–4.541, P< 0.001; ICGC cohort: HR = 2.146, 95% CI = 1.467–3.139, P< 0.001), while gender remained an independent prognostic factor for OS in the ICGC cohort (HR = 0.434, 95% CI = 0.225–0.838, P = 0.013). The risk score was significantly associated with OS solely in the TCGA cohort (HR = 2.767, 95% CI = 1.687–4.541, P<0.001) (Figure 4C and D). An amalgamation of the risk score, risk groups, and clinical stage was created as a nomogram for the TCGA cohort in order to aid in the clinical calculation of survival odds for HCC patients (Figure 4E). Additionally, calibration plots showed that the predictions of the nomograms matched the actual results fairly well (Figure 4F).

Figure 4 OS-related factors were identified, and a nomogram and calibration curve were developed. TCGA cohort (A and C), ICGC cohort (B and D). (A and B) OS-related factors were screened by Univariate Cox regression analyses. (C and D) OS-related factors were screened by Multivariate Cox regression analysis. (E) Nomogram incorporating GRGs and tumor stage. (F) Calibration curve predicting the 1-, 2-, and 3-year performance of the nomogram.

OS, Overall survival; P values were shown as: *P < 0.05; ***P< 0.001.

Risk Score of the Prognostic Model and Clinical Features

After examining the relationship between clinical characteristics and risk scores of participants with HCC, our study found no significant disparity in risk scores between patients aged 65 years or younger and those older than 65 years, and similarly no significant difference was observed in risk scores between male and female patients (Figure 5A and B). However, patients with higher tumor grades and tumor stages exhibited significantly elevated risk scores compared to those with lower tumor grades and stages (Figure 5C and D). Consistent with the findings from the TCGA cohort, there was no significant difference in risk scores between patients aged 65 years or younger and those older than 65 years, and no significant difference between male and female patients in the ICGC dataset (Figure 5E and F). Notably, patients with stage III–IV tumors exhibited a significant increase in risk scores (no information on the tumor’s grade in the ICGC dataset) (Figure 5G). Additionally, differential expression of ALDH2 were detected between patients<=65 years of age and those >65 years, and ALDH2 and ADH4 expression differ in female and male patients (P < 0.05, Supplementary Figure 4A and B). Furthermore, we demonstrated a significant upregulation in the expression of 4 GRGs in individuals with tumors rated 3–4 (P < 0.01, Supplementary Figure 4C). Notably, Patients with tumor stages III–IV had significantly increased expression of G6PC, ALDH2, and ADH4 (P < 0.05, Supplementary Figure 4D).

Figure 5 The risk score in different groups divided by clinical characteristics. TCGA cohort (A–D), ICGC cohort (E–G). (A and E) Age. (B and F) Gender. (C) Tumor grade. (D and G) Tumor stage.

KEGG Enrichment Analysis of Two Risk Groups

Using GSVA and GSEA, we explored the biological processes associated with the two identified risk groups. GSVA revealed that the high-risk group exhibited significantly higher levels of cytoplasmic DNA sensing, heparan sulfate production, and notch signaling (Supplementary Figure 5). Additionally, GSEA indicated that the high-risk group predominantly activated pathways such as neuroactive ligand-receptor interaction, cytokine receptor interaction, extracellular matrix (ECM) receptor interaction, primary immunodeficiency, and lysine degradation (Supplementary Figure 6A). In contrast, the low-risk group primarily activated pathways involved in cytochrome P450 drug metabolism, fatty acid metabolism, glycine/serine/threonine metabolism, and retinol metabolism (Supplementary Figure 6B).

TME and Immune Infiltration Analysis

The link between TME and tumor development is well-documented. We performed an immune infiltration analysis to further investigate the differences between high- and low-risk HCC sufferers. Using the CIBERSORT package of R, we assessed each sample’s relative amounts of various immune cell types. (Figure 6A). Strong relationships were found in our examination of immune cell interactions between active CD4 memory T lymphocytes and CD8+ T cells, functional assistance T cells and CD8-positive T cells, and M0 macrophages and CD8 T cells (Figure 6B). Regarding immune cell infiltration, notable distinctions were found between the high- and low-risk categories, particularly of M1 macrophages, resting CD4 memory T cells, and monocytes (Figure 6C, P< 0.05). Except for ADH4, the correlations between 3GRGs and immune cell subsets are demonstrated, with the risk score showing a pronounced link to the modulation of T cells, resting CD4 memory T cells, M0 macrophages, and M1 macrophages (Figure 6D). Tumor stemness is assessed using RNA stemness score (RNAss) and DNA stemness score (DNAss), which serve as independent indices.30 A negative association was found between the stromal score and DNAss, while correlation studies showed a strong positive relationship between the risk score and the immunological score and RNAss (Figure 6E, P < 0.05).

Figure 6 Relationship between immune microenvironment infiltration and risk scores for HCC. (A) The proportion of immune cells responding in HCC patients with different risk scores. (B) The interrelationship between immune cells. (C) Differences in immune cell levels between different risk groups. (D) Correlation between immune cell populations and four GRGs. (E) The relationship between risk score and RNAss, DNAss, Stromal Score, and Immune Score. P values were shown as: *P < 0.05; **P < 0.01.

Prediction of Anti-Tumor Drug Sensitivity

To further assess the potential therapeutic value of our four-GRGs prognostic model, we employed the GDSC platform to assess the anti-tumor drug responsiveness. Our analysis indicated that relative to the low-risk group, the high-risk group exhibited increased sensitivity to several chemotherapeutic and targeted agents, including 5-fluorouracil, gefitinib, dasatinib, afatinib, lapatinib, pictilisib, navitoclax, temozolomide, osimertinib, sapitinib, staurosporine, and alpelisib (Figure 7, P< 0.001), which suggests a higher level of therapeutic efficacy in high-risk patients, in contrast, six common anti-tumor medications showed increased susceptibility in the low-risk group of individuals, including sorafenib, oxaliplatin, axitinib, irinotecan, gemcitabine, and cytarabine (Supplementary Figure 7, P< 0.001). These findings facilitate more informed predictions for chemotherapy and targeted therapy selection in HCC patients based on their risk scores (Supplementary Table S3).

Figure 7 Application of GRGs for drug sensitivity prediction.

Cell Cycle and Angiogenesis Gene Analysis in Two Risk Groups

The hallmark of malignant tumors is the uncontrolled proliferation of cells stemming from the deregulation of the cell cycle.31,32 Given the cell cycle’s critical role in tumorigenesis, dysregulated expression of cell cycle genes is hypothesized to fuel the progression of HCC. We analyzed the expression patterns of 21 key cell cycle genes across two risk groups. Aside from GSK3B, our analysis revealed a marked upregulation of 20 cell cycle genes in the group with a higher risk, indicating that the four GRGs in our model probably promote the proliferation of tumor cells, thereby promoting tumor growth (Figure 8A). Roy et al highlighted angiogenesis as a crucial factor in HCC development.33 Consequently, we made a comparison of the expression levels of eight angiogenesis-related genes in high- and low-risk HCC groups. High-risk patients showed significantly higher expression of HIF1A, VEGFB, ROBO1, and SLIT1 than low-risk patients (Figure 8B). The ICGC group exhibited significantly higher expression levels of eight angiogenesis in tumors genes and 21 genes about cell cycle in the group with a higher risk (Figure 8C and D).

Figure 8 Expression of cell cycle genes and tumor angiogenesis genes in high- and low-risk groups for the TCGA cohort (A and B) and ICGC cohort (C and D). (A and C) Cell cycle genes in different groups. (B and D) Tumor angiogenesis genes in different groups. P values were shown as: ns, not significant; *P < 0.05; **P < 0.01; ***P< 0.001.

GRGs Subgroups Identification

To delve deeper into the expression patterns of GRGs involved in carcinogenesis, we employed an integrative approach making use of TCGA dataset HCC specimens. On the basis of the expression profiles of the four GRGs, the robust clustering method was applied to divide the HCC samples into different molecular subtypes. K=4 was determined through analysis of the CDF curve (Figure 9A). Subsequently, the TCGA cohort was divided into four discrete subgroups. Survival analysis indicated that subgroup B had a significantly longer OS compared to subgroups A, C, and D (Figure 9B). Furthermore, t-SNE analysis underscored the heterogeneity among the four subgroups (Figure 9C). The expression of GRGs, together with the corresponding clinicopathological characteristics, was shown in a heatmap for each of the four subgroups (Figure 9D).

Figure 9 Subgroups of HCC classified by GRGs. (A) The consensus clustering analysis. (B) Survival probabilities for four HCC subgroups. (C) t-SNE identified four distinct subgroups characterized by variations in the expression levels of GRGs. (D) Clinical and pathological features of four subgroups of GRG expression.

KEGG Enrichment Analysis of Four Subgroups

Given the pronounced differences between subgroups A and B, we conducted a detailed investigation into the differential enrichment of the KEGG pathways using the GSVA and GSEA tools. GSVA indicated that pathways such as Riboflavin metabolism, DNA replication, homologous recombination, and Purine metabolism were significantly enriched in subgroup A (Supplementary Figure 8A). In contrast, GSEA findings revealed that Fatty acid metabolism, Glycine/Serine/Threonine metabolism, Primary bile acid biosynthesis, and Propionate metabolism were particularly active in subgroup B, while these pathways were less active in subgroup A (Supplementary Figure 8B and C).

Immune Cell Infiltration in Four Subgroups

By displaying the expression profiles of each of the four subgroups, the GRG distribution was shown. It is noteworthy that group A had the lowest expression level of GRG among the four groups (Figure 10A). This is in line with the predictive survival patterns observed. Concerning immune cell infiltration, group C stood out from the others, having an extremely reduced proportion of bone marrow-derived inhibitory cells, activated dendritic cells, and regulatory t-cells than the other three distinct categories (Figure 10B). This highlights the close relationship between immune cells and gluconeogenesis in the prognostic subgroups of HCC.

Figure 10 Immunity and gene expression patterns of GRG subgroups. (A) 43 GRGs expression profiles. (B) Patterns of immune infiltration across four subgroups. P values were shown as: ns, not significant; *P < 0.05; **P < 0.01; ***P< 0.001.

GRGs Subgroups and Survival Status

Based on the findings, we identified four subgroups, two risk groups, and two clinical outcomes distributions and relationships. The Alluvial diagram demonstrates that subgroup B is associated with the most favorable prognosis, whereas subgroup A is linked to the poorest prognosis (Figure 11A). Furthermore, a significant disparity in risk scores was noted, with the most pronounced difference observed between subgroups A and B and the least between subgroups C and D. Importantly, high-risk scores were found to be significantly correlated with subgroup A and low-risk scores with subgroup B (Figure 11B).

Figure 11 HCC subgroups and risk, clinical outcomes. (A) Sankey diagrams for four subgroups, two risk groups, and two clinical outcomes. (B) Differences in risk scores between the four subgroups. The x-axis indicates the type of subgroup; the y-axis indicates the risk score.

Cell Cycle and Angiogenesis Gene Analysis in Four Subgroups

Recent research has demonstrated the critical roles that genes linked to the cell cycling and angiogenesis of tumors play in the etiology of HCC.33,34 Building on this, we conducted an in-depth analysis of cell cycling and angiogenesis of tumor genes in the four subgroups of HCC. Among the 21 key cell cycling cells, there were significant differences in expression across the four subgroups, and most cell cycle genes in subgroup A had higher expression levels compared to the other three subgroups, while subgroup B exhibited the lowest expression levels of cell cycle genes (Figure 12A). Given the highest expression of ADH4 in subgroup B, we deduced that ADH4 might have a connection to genes associated to the cell cycling, which could be important for the growth of HCC tumor cells. Regarding tumor angiogenesis genes, only VEGFB, ROBO1, and SLIT1 showed notable differences across the four subgroups, with subgroup A exhibiting the highest levels of gene expression and subgroup B showing the lowest (Figure 12B).

Figure 12 Expression of cell cycle genes and tumor angiogenesis genes in four-GRGs subgroups. (A) Cell cycle genes in different subgroups. (B) Tumor angiogenesis genes in different subgroups. P values were shown as: ns, not significant; **P < 0.01; ***P< 0.001.

G6PC Suppresses Cell Proliferation in HCC

To substantiate the theoretical implications of GRGs in tumor biology with experimental evidence, some experiments was conducted to confirm our findings. Prior research has established G6PC as a critical gene in the gluconeogenesis pathway, which is significantly reduced and functions as a tumor inhibitor in HCC.13,35 Our initial investigation revealed that G6PC expression in HCC tumor tissues is significantly lower than in adjacent normal liver tissues and is strongly associated with the prognosis of HCC patients (Figure 13A and B). Therefore, G6PC was designated as the candidate gene for further research. We successfully established a stable HepG2 cell line with G6PC overexpression, and immunofluorescence confirmed the enhanced expression of G6PC in the lentivirus-transfected HepG2 cells (Figure 13C). The experimental results indicate that G6PC overexpression significantly reduced the proliferative capacity of HepG2 cells (Figure 13D). Furthermore, G6PC overexpression substantially diminished the clonogenicity ability of HepG2 (Figure 13E and F). The findings underscore the role of G6PC as an adverse controller in HCC, impeding the malignant progression of the disease.

Figure 13 Overexpression of G6PC inhibits the proliferation of the HepG2 hepatocellular carcinoma cells. (A) Differential expression of G6PC between HCC tissues and normal tissues in TCGA cohort. (B) Prognosis of G6PC in TCGA cohort. (C) Immunofluorescence assays were used to observe the overexpression of G6PC. (D) Overexpression of G6PC inhibits proliferation of HepG2. (E and F) Overexpression of G6PC inhibits the clonogenicity ability of HepG2. NC, Negative control; P values were shown as: **P < 0.01; ***P< 0.001.

Discussion

HCC is recognized as a highly malignant neoplasm, presenting substantial challenges in cancer management. Research by Bian et al has indicated that the gluconeogenesis pathway may suppress the proliferation of HCC tumor cells.36 Therefore, our study aims to develop a robust prognostic model utilizing GRGs. Additionally, our objective is to examine the correlation between these genes and the prognosis of HCC patients, as well as their influence on tumor microenvironment infiltration and drug sensitivity in HCC.

The prognostic impact of gluconeogenesis in HCC patients has yet to be defined. Our analysis of both TCGA and ICGC cohorts revealed that patients with HCC who scored low on the risk score had better overall survival than those who scored high. Additionally, progressive tumor staging and grading were connected with patients’ higher risk assessment scores, suggesting that gluconeogenesis could serve as a candidate prognostic biomarker for HCC. Chen et al have indicated that gluconeogenesis may be a potential prognostic biomarker for HCC, supporting our findings to some extent.37 We have established a model involving four genes closely linked to gluconeogenesis. Reliability and accuracy were validated using the ICGC cohort. Additionally, we underscored the model’s efficacy in predicting HCC patient outcomes by creating a nomogram and employing calibration curves for model refinement. Prior research has also recognized the role of GRGs in HCC progression. Khan et al reported that upregulation of PCK1 can promote the conversion from glycolysis to the gluconeogenesis pathway,16 thus inhibiting HCC cell proliferation. Downregulation of FBP1 can induce altered glucose metabolism, leading to tumor progression and unfavorable prognosis in HCC, according to Hirata et al.38 To further elucidate the association between GRGs and HCC carcinogenicity, we conducted a consensus clustering analysis. While previous studies by Chen et al and Zheng et al have explored the relationship between glycolysis-related genes and HCC, categorizing HCC patients into two groups based on the expression levels of these genes,37,39 our study differentiated all 365 hCC samples into four subgroups based on GRG expression levels. This stratification may enhance the development of personalized treatment strategies for HCC patients. Survival analysis indicated that patients in subgroup B enjoyed the most favorable survival outcomes, while those in subgroup A fared the worst. A heatmap integrating GRG expression with clinical characteristics across subgroups revealed that subgroup A had significantly lower expression levels of the four GRGs compared to the other groups, suggesting a correlation between GRG expression and HCC carcinogenicity, with lower GRG expression levels implying higher HCC malignancy and a more unfavorable prognosis for patients.

Prior research on the relationship between gluconeogenesis and tumor microenvironment infiltration in HCC was limited. Therefore, the CIBERSORT software was used to perform a thorough analysis. Our results showed that those with low risks had considerably higher invasion rates of M1 macrophages and CD4 memory resting T cells than those with a high risk did. This suggests that an elevation of these immune cells may help to reduce the tumorigenicity of HCC. Studies have shown that increased M1 macrophage presence can postpone HCC recurrence, leading to improved patient prognosis,40 and CD4 T cells play a role in clearing liver cells in precancerous conditions, thereby inhibiting HCC cell proliferation.41 These findings align with our research outcomes. Our findings reveal a significant correlation between G6PC and various immune cells, including CD4+ memory resting T cells, M0 macrophages, and naive B cells. Despite the absence of previous studies exploring the relationship between G6PC and these specific immune cell types, our research suggests that G6PC may modulate the progression of HCC by exerting an influence on immune cells.

Upon analyzing the immune cell infiltration across the four subgroups, it was observed that Subgroup C demonstrated significantly lower proportions of activated B cells, CD8 T cells, dendritic cells (DCs), and myeloid-derived suppressor cells compared to the other subgroups. Previous studies by Liu et al and Luo et al have indicated that the depletion or suppression of CD8 T cells may contribute to the progression and metastasis of HCC.42,43 DCs, being professional antigen-presenting cells (APCs), are crucial for the activation of anti-tumor immune responses.44 Chen et al have shown that DC-based therapy not only enhances anti-tumor immunity but also improves the survival rate and prolongs the survival time of HCC patients.45 These findings corroborate our study’s conclusion that the poorer prognosis of patients in Subgroup C may be attributed to the insufficient infiltration of anti-tumor immune cells. In summary, our research underscores the pivotal role of gluconeogenesis in the TIME of HCC, influencing the infiltration of immune cells and, consequently, the clinical outcomes of HCC patients.

The relationship between gluconeogenesis and tumor drug resistance remains understudied. Xu et al reported that the high-risk group exhibited resistance to Gefitinib, Lapatinib, and 5-Fluorouracil treatments.46 In contrast, our findings indicate that patients within the high-risk group showed increased sensitivity to these three immunotherapeutic and chemotherapeutic agents. Meanwhile, the low-risk group exhibited increased sensitivity to six medications, among which oxaliplatin and gemcitabine are integral to HCC chemotherapy protocols,47 and sorafenib is the frontline treatment for advanced HCC patients.48 Based on our results, our four-GRG prognostic prediction model could reliably predict the responsiveness of HCC to chemotherapy and targeted agents, offering tailored treatment strategies for HCC patients across different risk groups.

The existing literature confirms that all four genes are significantly associated with HCC in our prognostic model for GRGs. G6PC is a recognized key gene in gluconeogenesis, involved in catalyzing the formation of glucose-6-phosphate. Wang et al reported that G6PC deficiency in HCC leads to decreased serum glucose levels and glucose-6-phosphate accumulation, which promotes tumor cell proliferation.49 Li et al found that G6PC serves as a potential prognostic biomarker in HCC, with low expression significantly associated with poor overall survival in HCC patients.50 These findings align with our prognostic model and experimental results. PPARGC1A, an upstream regulator of G6PC, activates the gluconeogenic pathway by modulating G6PC expression.35,51 Zhang et al indicated that PPARGC1A is underexpressed in HCC and associated with a poor prognosis,52 a finding corroborated by our study. The ADH4 gene, part of the alcohol dehydrogenase superfamily, was identified by Liu et al as an independent factor for improved prognosis in HCC, with high expression significantly associated with better survival rates in HCC patients.53 Our results are consistent with these findings. ALDH2, a key enzyme in ethanol metabolism, has been shown in previous studies to be downregulated in HCC cells, causing migration and invasion, and low ALDH2 expression is linked to poor overall survival in HCC patients.54,55 Our findings underscore its essential role in HCC tumor progression.

Despite a clear understanding of the cell cycle and angiogenesis in tumorigenesis, research on the interplay between gluconeogenesis and these processes in HCC is limited. Cyclin B1 (CCNB1), a key promoter of mitosis, exhibits periodic expression throughout the cell cycle.56 Cyclin B2 (CCNB2), a member of the cyclin family, plays a crucial role in the G2/M phase transition.57 Cyclin E (CCNE1) acts as a positive regulator of cell cycle progression, facilitating the G1/S phase transition.58 The cell-division cycle protein 20 (CDC20) is instrumental in chromosome segregation and mitotic exit.59 Studies have shown that elevated expression of cell cycle genes, including CCNB1, CCNB2, CCNE1, and CDC20, is associated with the initiation and progression of HCC, as well as tumor relapse and metastasis.60–63 According to our research, the group with a higher risk had significantly higher expression degrees of CCNB1, CCNB2, CCNE1, and CDC20 than the group with lower risk did. Roundabout Guidance Receptor 1 (ROBO1) is the member of the Roundabout family of receptors and the SLIT/ROBO axis is genetically and epigenetically altered in cancer, promoting proliferation, metastasis, and angiogenesis.64 Vascular endothelial growth factor B (VEGFB) is highly expressed in vascular endothelial cells and other cell types, with numerous studies reporting its anti-angiogenic and anti-tumor effects.65,66 Rossi et al emphasized the critical role of VEGFB in the initiation or progression of HCC,67 while H et al identified overexpression of ROBO1 in HCC.68 Our results corroborate these findings, showing increased expression levels of VEGFB and ROBO1 genes in the high-risk group. In general, our study results indicate a tight link between the four-GRG prognosis model, cell cycling and angiogenesis genes in the progression of HCC.

While our study has yielded significant findings, it is not without limitations. Firstly, the research is constrained by reliance on data from public databases, specifically two datasets, which may restrict the statistical power and generalizability of our results due to a limited sample size. It is imperative to incorporate additional large-sample datasets for the purpose of validating enhance the validity of our prognostic model. Secondly, our conclusions are based on bioinformatic analyses and corroborated by rudimentary cellular experiments; however, further in vitro and in vivo experimental evidence is necessary to bolster the credibility of our findings. Lastly, the study did not explore the molecular mechanisms by which gluconeogenesis influences TIME and drug sensitivity in HCC, and further studies are needed to elucidate these mechanisms. Despite these limitations, our study highlights the considerable utility of the four-GRGs prognostic model, establishing a link between gluconeogenesis, TIME, patient prognosis, and drug treatment strategies in HCC. Our 4-GRG model promises to be a predictive instrument for evaluating the prognosis of patients with HCC, as well as tumor microenvironment infiltration, and drug sensitivity, pointing towards innovative approaches for the clinical, immunological, and pharmacological management for HCC.

Conclusion

In summary, we have developed a prognostic model incorporating 4-GRGs that effectively predicts the prognosis, tumor microenvironment infiltration, and drug sensitivity in HCC patients. This study also identified two distinct risk groups and four molecular subtypes (A, B, C, and D) related to gluconeogenesis in HCC, which exhibit significant differences in gene expression profiles, immune cell infiltration, and clinical outcomes. Preliminary cellular experimental data indicate that G6PC, a component of the 4-GRGs, can suppress the proliferation and migration of HCC tumor cells. These findings may guide the development of innovative strategies and therapeutic interventions for the clinical management of HCC.

Data Sharing Statement

The public databases analyzed in this study can be found in the Cancer Genome Atlas (https://portal.gdc.cancer.gov/) and The International Cancer Genome Consortium (https://dcc.ICGC.org/projects/LIRI-JP).

Ethical Approval and Consent to Participate

All participants in the TCGA and ICGC cohorts had signed informed consent. The study was approved by the ethics committee of Sun Yat-sen Memorial Hospital (Approval number: SYSEC-KY-KS-2021-209) and conducted in compliance with the Declaration of Helsinki.

Acknowledgments

We thank the TCGA and ICGC projects for gene expression and clinical data.

Author Contributions

All authors made a significant contribution to the work reported, whether that is in the conception, study design, execution, acquisition of data, analysis and interpretation, or in all these areas; took part in drafting, revising or critically reviewing the article; gave final approval of the version to be published; have agreed on the journal to which the article has been submitted; and agree to be accountable for all aspects of the work.

Funding

This work was supported by the National Natural Science Foundation of China (82073062).

Disclosure

Authors state no conflict of interest.

References

1. Sung H, Ferlay J, Siegel RL, et al. Global cancer statistics 2020: GLOBOCAN estimates of incidence and mortality worldwide for 36 cancers in 185 countries. CA Cancer J Clin. 2021;71(3):209–249. doi:10.3322/caac.21660

2. Llovet JM, Kelley RK, Villanueva A, et al. Hepatocellular carcinoma. Nat Rev Dis Primers. 2021;7(1):6. doi:10.1038/s41572-020-00240-3

3. Dopazo C, Søreide K, Rangelova E, et al. Hepatocellular carcinoma. Eur J Surg Oncol. 2024;50(1):107313. doi:10.1016/j.ejso.2023.107313

4. Yang JD, Hainaut P, Gores GJ, Amadou A, Plymoth A, Roberts LR. A global view of hepatocellular carcinoma: trends, risk, prevention and management. Nat Rev Gastroenterol Hepatol. 2019;16(10):589–604. doi:10.1038/s41575-019-0186-y

5. Imamura H, Matsuyama Y, Tanaka E, et al. Risk factors contributing to early and late phase intrahepatic recurrence of hepatocellular carcinoma after hepatectomy. J Hepatol. 2003;38(2):200–207. doi:10.1016/s0168-8278(02)00360-4

6. Wu B, Yuan Y, Liu J, et al. Single-cell RNA sequencing reveals the mechanism of sonodynamic therapy combined with a RAS inhibitor in the setting of hepatocellular carcinoma. J Nanobiotechnology. 2021;19(1):177. doi:10.1186/s12951-021-00923-3

7. Krishnan MS, Rajan Kd A, Park J, et al. Genomic analysis of vascular invasion in HCC reveals molecular drivers and predictive biomarkers. Hepatology. 2021;73(6):2342–2360. doi:10.1002/hep.31614

8. Vander Heiden MG, Cantley LC, Thompson CB. Understanding the Warburg effect: the metabolic requirements of cell proliferation. Science. 2009;324(5930):1029–1033. doi:10.1126/science.1160809

9. Shang RZ, Qu SB, Wang DS. Reprogramming of glucose metabolism in hepatocellular carcinoma: progress and prospects. World J Gastroenterol. 2016;22(45):9933–9943. doi:10.3748/wjg.v22.i45.9933

10. Yang Z, Yu M, Li X, et al. Retinoic acid inhibits the angiogenesis of human embryonic stem cell-derived endothelial cells by activating FBP1-mediated gluconeogenesis. Stem Cell Res Ther. 2022;13:239. doi:10.1186/s13287-022-02908-x

11. Sharabi K, Tavares CDJ, Rines AK, Puigserver P. Molecular pathophysiology of hepatic glucose production. Mol Aspects Med. 2015;46:21–33. doi:10.1016/j.mam.2015.09.003

12. Tian H, Zhu X, Lv Y, Jiao Y, Wang G. Glucometabolic reprogramming in the hepatocellular carcinoma microenvironment: cause and effect. Cancer Manag Res. 2020;12:5957–5974. doi:10.2147/CMAR.S258196

13. Wang Z, Dong C. Gluconeogenesis in cancer: function and regulation of PEPCK, FBPase, and G6Pase. Trends Cancer. 2019;5(1):30–45. doi:10.1016/j.trecan.2018.11.003

14. Liu MX, Jin L, Sun SJ, et al. Metabolic reprogramming by PCK1 promotes TCA cataplerosis, oxidative stress and apoptosis in liver cancer cells and suppresses hepatocellular carcinoma. Oncogene. 2018;37(12):1637–1653. doi:10.1038/s41388-017-0070-6

15. Liao K, Deng S, Xu L, et al. A feedback circuitry between polycomb signaling and fructose-1, 6-Bisphosphatase enables hepatic and renal tumorigenesis. Cancer Res. 2020;80(4):675–688. doi:10.1158/0008-5472.CAN-19-2060

16. Khan MW, Biswas D, Ghosh M, Mandloi S, Chakrabarti S, Chakrabarti P. mTORC2 controls cancer cell survival by modulating gluconeogenesis. Cell Death Discov. 2015;1:15016. doi:10.1038/cddiscovery.2015.16

17. Fu T, Dai LJ,