A Novel Prognostic Signature Integrating Immune and Glycolytic Pathways for Enhanced Prognosis and Immunotherapy Prediction in Hepatocellular Carcinoma

Introduction

Hepatocellular carcinoma (HCC) ranks as the fourth leading cause of cancer-related mortality, contributing to approximately 8.3% of all cancer deaths globally.¹ In recent years, advancements in alternative therapies and immunotherapy have brought hope for improving the prognosis of HCC patients. For instance, metabolic reprogramming and targeted inhibition of key signaling pathways have emerged as potential strategies to overcome treatment resistance and enhance survival rates.² Immunotherapy, particularly immune checkpoint inhibitors (ICIs) such as anti-PD-1 and anti-CTLA-4, has transformed the treatment landscape for HCC by enhancing anti-tumor immune responses.³ However, the prognosis for HCC remains bleak, particularly for patients in advanced stages, where the median survival is less than 1 year.⁴ This grim outlook underscores the urgent need for reliable prognostic biomarkers to more accurately predict clinical outcomes, enable personalized treatment approaches, and ultimately improve patient survival rates.

One of the defining characteristics of tumors, including HCC, is an increase in glycolysis, heightened lactate production, and disruption of the tumor immune microenvironment (TIME). Tumor cells often exhibit elevated glycolytic activity, relying heavily on glycolysis for energy production even in the presence of oxygen, a phenomenon known as the Warburg effect.⁵ This metabolic shift results in the generation of substantial amounts of lactic acid, which accumulates and leads to the acidification of the TIME. Such acidification creates a conducive environment for tumor cells to evade immune detection and attack.⁶ In recent years, research has increasingly focused on regulating TIME by manipulating glycolysis and enhancing lactate production, positioning this as a promising area for developing novel therapeutic strategies. Studies have shown that the overexpression of key glycolytic enzymes (HK2, PKM2, and LDHA) in HCC is closely associated with tumor malignancy and poor prognosis.⁷ Immunosuppressive cells, including tumor-associated macrophages (TAMs), regulatory T cells (Tregs), and myeloid-derived suppressor cells (MDSCs), extensively infiltrate the HCC microenvironment. They suppress anti-tumor immune responses by secreting immunosuppressive cytokines (TGF-β, IL-10) and expressing immune checkpoint molecules (PD-L1).⁸ Additionally, lactate produced by HCC cells through glycolysis can drive the M2 polarization of TAMs by activating the MCT-HIF1-α pathway, or directly inhibit the function of cytotoxic T cells (CTLs) and natural killer (NK) cells, further weakening the immune system’s anti-tumor capabilities.⁹ A diagram was plotted to illustrate the recognized pathways related to immunity and glycolysis in HCC (Supplementary Figure 1). Nonetheless, the complex biochemical and molecular mechanisms governing the Warburg effect within the TIME of HCC demand further investigation into immune- and glycolysis-related genes.

In this study, the characteristic gene set of HCC was identified using weighted gene co-expression network analysis (WGCNA). By intersecting the immune-related gene set, glycolytic gene set, and HCC characteristic gene set, a total of 196 immune and glycolytic-related genes (IGRGs) ere obtained. These 196 IGRGs underwent further screening, and ultimately, 13 key survival-related immune and glycolytic-related genes (SRIGRGs) were identified through a combination of F-test, univariate Cox regression, and Lasso regression analysis. An immune-glycolysis-related prognostic signature (IGRPS) was then constructed and validated in both the training and validation cohorts. The Kaplan-Meier (K-M) survival curves, time-dependent receiver operating characteristic (ROC) curves, and the concordance index (C-index) of the prognostic nomogram demonstrated excellent performance in both cohorts, highlighting the superior predictive capability of the IGRPS. Gene set enrichment analysis (GSEA) revealed significant differences in multiple biological processes and pathways between the high- and low-risk groups. TIME analysis showed that the infiltration of immune cells, such as activated dendritic cells (DCs), mast cells, and natural killer (NK) cells, along with various immune-related functions, differed significantly between these groups. Further analyses examined the association between the risk score and clinical features, immune subtypes, glycolysis profiles, immune escape potential, and response to immunotherapy. Finally, the expression levels and functional roles of the genes within the IGRPS were validated through in vitro experiments. In summary, we constructed a valuable prognostic model based on IGRPS and explored the relationship of the model with TIME which might provide insights for physicians to develop personalized treatment strategies for HCC patients.

Materials and Methods

Data Acquisition

RNA-seq data, clinical information, and gene mutation profiles for 424 samples, including 374 hCC specimens and 50 adjacent non-tumor tissue samples, were acquired from The Cancer Genome Atlas (TCGA) database (https://www.cancer.gov/tcga). The mRNA transcriptomic data (GSE76427) and corresponding clinical data for 167 matched pairs of HCC and adjacent tissues were sourced from the Gene Expression Omnibus (GEO) database (https://www.ncbi.nlm.nih.gov/geo/).¹⁰ Additionally, single-cell sequencing datasets GSE175793, comprising three HCC specimens and three non-tumor tissue samples, were also obtained. All data followed strict screening criteria: a) samples with a missing data rate exceeding 20% were excluded, b) HCC samples without prior chemotherapy or radiotherapy were selected. The TCGA cohort served as the training group, while the GEO cohort was used as the validation group to construct the prognostic model. Lists of immune-related genes were retrieved from the ImmPort (https://www.immport.org/resources), InnateDB (https://www.innatedb.ca/), and GSEA Molecular Signatures databases (https://www.gsea-msigdb.org/gsea/msigdb). Concurrently, a gene set associated with glycolysis was identified from the GSEA Molecular Signatures database.

Identification of HCC IGRGs

Characteristic genes of HCC were identified using WGCNA in the TCGA cohort. First, we calculated the network properties, fit indices (R²), and average connectivity for each soft threshold (β=1 to 20). The β value with R² above 0.8 and moderate average connectivity density was selected as the optimal soft-thresholding power value. A scale-free co-expression network was constructed based on the optimal soft-thresholding power value, which was then converted into an adjacency matrix and a topological overlap matrix (TOM). Pearson correlation analysis was used to integrate disease traits with gene modules, and it was found that the disease trait-gene module correlation exhibited minimal variation under the optimal soft threshold, indicating high stability. HCC characteristic modules were defined based on a correlation coefficient greater than ∣0.09∣ and a P-value less than 0.05. Genes from these characteristic modules were extracted and classified as HCC characteristic genes. Finally, IGRGs were identified by taking the intersection of the HCC characteristic gene set, the immune-related gene set, and the glycolytic-related gene set.

F-Test, Univariate Cox Regression, and Lasso Regression Analyses

IGRGs with scores exceeding 60 were initially selected using an F-test and subsequently analyzed with Univariate Cox regression to identify SRIGRGs. The results were visualized in a forest plot. To further refine the selection and prevent overfitting, Lasso regression was employed. Using a 10-fold cross-validation approach to minimize the risk of overfitting, the genes associated with the optimal λ value were selected. This stepwise screening process, incorporating F-test, Univariate Cox regression, and Lasso regression, enabled the identification of crucial SRIGRGs critical for constructing the predictive model.

Construction of IGRPS

Multivariate Cox regression was performed to construct IGRPS in the training group. The risk score of each patient was calculated as follows: risk score = gene A expression × gene A regression coefficient + gene B expression × gene B regression coefficient +. + Gene N expression × gene N regression coefficient. The accuracy of the model was verified in the validation group. The “survminer” and “survival” R packages were then applied to obtain the median value of the risk score, which was used to divide HCC patients into high-risk and low-risk groups.

The Predictive Value and Clinical Efficacy of IGRPS

The prognostic value of the IGRPS was evaluated and validated in both the training and validation groups. K-M survival curves, time-dependent ROC curves, and a prognostic nomogram were generated using the “maxstat”, “survival”, and “pROC” R packages. Additionally, the nomogram was constructed with the “rms” and “survivalROC” packages to further assess the model’s predictive accuracy. Finally, both univariate and multivariate Cox regression analyses were conducted to determine whether the risk scores could serve as independent prognostic and predictive factors for HCC.

Correlation Analysis of Risk Score with Clinical Features, Immune Subtypes, and Glycolysis Profiles

The “ComplexHeatmap” package was utilized to analyze the correlation between the risk score and various clinical features, including age, gender, tumor grade, clinical stage, and TNM stage. The “RColorBrewer” package was employed to investigate the association between the risk score and different immune subtypes. Additionally, glucose transporters and essential kinases (glycolysis-related genes) that encode and directly regulate cellular glycolysis were identified from previous studies.¹¹ The differential expression of these glycolysis-related genes between high-risk and low-risk groups was then analyzed to further understand their roles in the prognosis of HCC.

GSEA and Genomic Stability Analysis

To explore the potential molecular mechanisms of IGRGs in HCC, GSEA was conducted to identify the relevant functions and pathways. The top five pathways most closely associated with the risk score were selected and visualized to provide insights into their biological significance. Additionally, genetic mutation data of HCC samples were obtained from the TCGA database. The “maftools” R package was used to analyze and visualize the genetic mutation profiles in the high-risk and low-risk groups, allowing for a comprehensive comparison of mutation patterns between these groups.

TIME Analysis

The single-sample gene set enrichment analysis (ssGSEA) algorithm was employed to compute the relative infiltration abundance of 16 types of immune cells in each sample. The differences in immune cell infiltration between the high-risk and low-risk groups were visualized using a boxplot. Immune cells associated with overall survival (OS) were identified using the “limma”, “survival”, and “survminer” R packages, and their significance was visualized through K-M survival curves. To analyze differences in immune-related functions between the high-risk and low-risk groups, the “CIBERSORT”, “limma”, “reshape2”, and “ggpubr” packages were applied. Additionally, the correlation between the risk score and Stromal scores, Immune scores, and ESTIMATE scores was assessed using the “ESTIMATE” package.

Furthermore, the glycolytic metabolic pathway reshapes the local TIME. Immune cells such as γδ T cells undergo metabolic reprogramming to trigger metabolic adaptation, meeting their own needs and exerting immunosuppressive effects.¹² Immune cell infiltration related to the IGRPS and the expression levels of IGRPS genes in γδ T cells was analyzed using single-cell sequencing data. The GSE175793 dataset was normalized, and high-variance genes were selected with the “Seurat” R package. Dimensionality reduction was performed through principal component analysis (PCA) and visualized using the uniform manifold approximation and projection (UMAP) method. Cell clustering was then conducted using the “FindNeighbors” and “FindClusters” functions. The “FindAllMarkers” function was utilized to identify marker genes for each cell cluster. Finally, cell types were automatically annotated using the “SingleR” package, based on reference transcriptome data.

Validation of Immunotherapy Response

Immunophenotype scores (IPS) for 378 hCC patients were obtained from The Cancer Immunome Atlas (TCIA) (https://www.tcia.at/home). The “ggpubr” R package was used to analyze the differential efficacy of anti-CTLA-4, anti-PD-1 monotherapies, and their combination therapy between high-risk and low-risk groups. Additionally, clinical data, molecular characteristics, and survival data of patients treated with different immunotherapy regimens were sourced from the IMvigor 210 database (https://github.com/ixxmu/mp_duty/issues/4275). The “limma” and “ggpubr” packages were utilized to evaluate the effectiveness of these immunotherapies.

Moreover, Tumor Immune Dysfunction and Exclusion (TIDE) data for HCC were downloaded from the TIDE database (http://tide.dfci.harvard.edu/). The “limma” and “ggpubr” packages were applied to investigate the association between the risk score and various TIDE components, including dysfunction, exclusion, microsatellite instability, and overall TIDE score. Finally, the tumor mutation burden (TMB) for each sample was calculated, and differences between high-risk and low-risk groups were compared. The correlation between TMB and risk scores was also analyzed to provide further insights into the impact of mutational load on prognostic outcomes and immunotherapy response in HCC.

Validation of Expression Levels of IGRGs

The differential expression of IGRGs was assessed between HCC and adjacent non-tumor tissues, as well as between high-risk and low-risk groups. The “ggpubr” and “ggExtra” packages were employed to analyze the correlation between the expression of IGRGs and the risk scores. The “pROC” package was utilized to construct ROC curves to determine the diagnostic value of each IGRG.

A total of 15 paired HCC and adjacent non-tumor tissues were collected from patients in the Department of Hepatobiliary Surgery at Shandong University Qilu Hospital. This study adhered strictly to the principles outlined in the Declaration of Helsinki, with ethical approval granted by the Ethics Committee of Shandong University Qilu Hospital. Informed consent was obtained from all patients prior to surgery. The expression levels of IGRPS genes at both the mRNA and protein levels were evaluated using reverse transcription-quantitative polymerase chain reaction (RT-qPCR) and immunohistochemistry (IHC), respectively.

Cell Culture and Transfection

Human HCC cell lines, Bel-7402 (No.C6109) and Huh-7 (No.SCSP-526), were obtained from the Cell Bank of the Chinese Academy of Sciences. Bel-7402 cells were cultured in RPMI-1640 (Solarbio, No.11875) medium supplemented with 10% fetal bovine serum (FBS) (Solarbio, No.S9030) and 1% penicillin-streptomycin (Solarbio, No.P1410), while Huh-7 cells were maintained in DMEM (Solarbio, No.11965) with the same supplements (10% FBS and 1% penicillin-streptomycin). All cell lines were maintained at 37°C in a humidified atmosphere containing 5% CO₂. Small interfering RNA (siRNA) (Beijing SyngenTech) transfections were conducted using Lipofectamine 3000 (Solarbio, No.L3000008).

Determination of Glucose Consumption, ATP, Lactic Dehydrogenase Activity, Lactate, and Extracellular Acidification Rate (ECAR) in HCC Cells

Bel-7402 and Huh-7 cell lines (5×10⁵ cells/well) were cultured in 96-well plates (Servicebio, No.KB-96U20-L) and transfected with either negative control siRNA (NC-siRNA) or specific siRNAs targeting PRKAG1 and B3GAT3. Once the cells reached 80% confluence, the cell supernatant was collected for further analysis. The levels of glucose consumption, ATP, lactate dehydrogenase activity, and lactate were measured using a glucose consumption detection kit (Abbkine, No.KTB1300), an ATP detection kit (Abbkine, No.KTB1016), a lactate dehydrogenase activity detection kit (Abbkine, No.KTB1110), and a lactate detection kit (Abbkine, No.KTB1100), respectively. Additionally, the ECAR was assessed using the Seahorse XFe96 Flux Analyzer to evaluate glycolytic activity.

Cell Proliferation and Migration Ability Assays

Cell proliferation was assessed using the Cell Counting Kit-8 (CCK-8) assay (Abbkine, No.KTA1020). Bel-7402 and Huh-7 cell lines were seeded into 96-well plates at a density of 2000 cells per well and transfected with either NC-siRNA or specific siRNAs targeting PRKAG1 (siRNA1, 2) and B3GAT3 (siRNA1, 2). The experiment included nine groups for each cell line: NC-siRNA, PRKAG1 siRNA1, PRKAG1 siRNA2, B3GAT3 siRNA1, B3GAT3 siRNA2, lactate added to NC-siRNA, lactate added to PRKAG1 siRNA1, lactate (Servicebio, No.GM3029) added to PRKAG1 siRNA2, and 2-deoxy-D-glucose (2-DG) (Solarbio, No.D8930) added to NC-siRNA or siRNA1, 2. After culturing for various time points, CCK-8 reagent was added to each well. Absorbance was measured at a wavelength of 450 nm using a microplate reader (BioTek), and cell viability was determined based on the optical density (OD) values.

For the migration assay, cells were seeded in the upper chamber of a 24-well Transwell plate with an 8-μm porous polycarbonate membrane (Servicebio, No.WG3422) at a density of 1×10⁵ cells per well. The upper chamber was filled with medium without FBS, while the lower chamber contained medium supplemented with 10% FBS. The experiment was designed with the same nine groups for each cell line as in the proliferation assay. After incubation, cells that had migrated to the lower chamber were fixed and stained. The number of migrating cells was then counted under an inverted microscope (Thermo Fisher Scientific) to assess the migration capability of each group.

In the scratch assay, cells were seeded onto a 6-well plate (Servicebio, No.CCP-6H) and cultured for 24 hours to allow them to adhere to the surface. A sterile 200-µL pipette tip was then used to create a straight scratch through the middle of the cell monolayer, ensuring that the width of the scratch was as uniform as possible. Images of the scratched area were captured using an inverted microscope at the start of the experiment (0 hours) and after 24 hours of incubation to assess cell migration into the scratch area.

Immunofluorescence (IF) was performed to evaluate the proliferative capacity of HCC cells at the protein level. The cells were fixed, permeabilized, and blocked to prepare for staining. An anti-Ki-67 antibody (Proteintech Group, No.27309-1-AP) was added to the cells, and they were incubated overnight to detect the proliferation marker. The following day, a fluorescently labeled secondary antibody was applied and incubated for 1 hour. Finally, the cells were mounted with a mounting medium and observed under a fluorescence microscope to capture images, which were used to assess the expression of Ki-67 as an indicator of cell proliferation.

Statistical Analysis

All statistical analyses were conducted using R version 4.3.1 and GraphPad Prism version 8.0. Differences between groups were evaluated using either a one-way analysis of variance (ANOVA) or a T-test, where appropriate. A p-value of less than 0.05 (P<0.05) was considered statistically significant.

Results

Screening of IGRGs

Eleven distinct HCC characteristic modules were identified using WGCNA, from which 4016 hCC characteristic genes were extracted (Figure 1A-C). A total of 11,429 immune-related genes were downloaded from the ImmPort, InnateDB, and GSEA Molecular Signatures databases, while 930 glycolytic-related genes were obtained from the GSEA Molecular Signatures database. The intersection of these three gene sets yielded 196 IGRGs (Figure 1D).

Figure 1 Identification of 13 key SRIGRGs. (A-C) HCC characteristic genes were screened by WGCNA. The optimal power value was 7, the color depth in Figure C represented the correlation, and the numbers in brackets were p-values. (D) 196 IGRGs were obtained from the intersection of the HCC characteristic gene set, immune-related gene set, and glycolytic gene set. (E) 24 SRIGRGs were obtained by univariate Cox regression. (F and G) Lasso regression selected the optimal λ value 13 and obtained 13 key SRIGRGs.

Additionally, 28 genes with scores greater than 60 were identified using the F-test (Supplementary Table 1). From these, 24 SRIGRGs were found to be associated with prognosis through univariate Cox regression analysis (Figure 1E). The K-M survival curves of these SRIGRGs are presented in Supplementary Figure 2. Lasso regression analysis, utilizing an optimal λ value of 13, further narrowed down the critical genes to ALDH2, PCK1, FBP1, AAAS, GMPPA, PRKAG1, PSMC4, HAX1, ADH1A, EIF6, B3GAT3, BPGM, and HSPA5 (Figure 1F and G).

Establishment of IGRPS

Multivariate Cox regression analysis was conducted to develop a prognostic signature based on three key SRIGRGs (ALDH2, PRKAG1, and B3GAT3) within the training group. The risk score for each patient was calculated using the following formula:

The positive regression coefficients of PRKAG1 and B3GAT3, as well as their hazard ratios (HR) greater than 1 in the univariate Cox analysis, suggested that these genes might function as oncogenic drivers in HCC. In contrast, ALDH2, with a negative coefficient, appeared to act as a protective factor (Figure 1E).

Validation of Predictive Value and Clinical Efficacy of IGRPS

Patients in the training group were categorized into high-risk and low-risk groups based on the median risk score. K-M survival analysis revealed that OS was significantly longer for patients in the low-risk group compared to those in the high-risk group, in both the training and validation cohorts (Figure 2A and B). The area under the curve (AUC) values for both the training and validation groups exceeded 0.7, indicating good prognostic performance (Figure 2C and D).

Figure 2 Validation of IGRPS. (A and B) Survival curves of the training group and the verification group. (C and D) 1, 3, and 5-year time-ROC curves for the training and validation groups. (E and F) Univariate and multivariate Cox analysis confirmed that IGRPS was an independent predictor and prognostic factor of HCC. (G and H) Prognostic nomograms and calibration curves for the training group.

Additionally, univariate and multivariate Cox regression analyses were conducted using the training group data. The univariate Cox analysis indicated that both clinical stage (HR = 1.680, P < 0.001) and risk score (HR = 1.893, P < 0.001) were independently associated with OS. In multivariate Cox analysis, the risk score remained an independent predictor of survival for HCC (HR = 1.675, P < 0.001) (Figure 2E and F). These findings confirmed that the IGRPS was an independent prognostic factor for HCC.

The IGRPS prognostic nomogram further demonstrated strong predictive performance, with a C-index of 0.73 for the training group (95% CI: 0.67–0.79, P = 7.58e-13). Calibration curves indicated that the nomogram had moderate accuracy in predicting outcomes (Figure 2G and H). Overall, these results validated the excellent predictive value of the IGRPS.

Finally, we developed a prognostic prediction website for HCC (https://xever.shinyapps.io/hcc_igrps/), named the HCC IGRPS Model, to enhance the clinical applicability of this model. By inputting the expression data of ALDH2, PRKAG1, and B3GAT3, the website can predict the prognostic risk for the person.

Association of Risk Score with Clinical Features, Immune Subtypes, and Glycolysis Profiles

The results demonstrated that the risk score was significantly associated with tumor grade, clinical stage, and T stage. Patients in the low-risk group were predominantly found in Grade 1, Grade 2, Stage I, and T1, whereas those in the high-risk group were mainly concentrated in Grade 3, Grade 4, Stage II, Stage III, T2, T3, and T4 (Supplementary Figure 3A-D). Furthermore, a significant correlation was observed between the risk score and immune subtypes; the low-risk group was primarily composed of C3 and C4 subtypes, while the high-risk group was mainly composed of C1 and C2 subtypes (Supplementary Figure 3E). Analysis of the glycolysis profile revealed that the expression levels of glycolysis-related genes were higher in the high-risk group compared to the low-risk group, indicating that glycolytic activity increased with the risk score (Supplementary Figure 3F and G).

GSEA and Genomic Stability Analysis

The GSEA results showed distinct pathway enrichment patterns between the high-risk and low-risk groups in both the GO and KEGG datasets. In the GO dataset, the high-risk group was significantly associated with pathways related to B cell-mediated immunity, complement activation, digestion, humoral immune response mediated by circulating immunoglobulins, and immunoglobulin complexes. Conversely, the low-risk group was linked to processes such as drug catabolic process, drug metabolic process, aromatase activity, oxidoreductase activity acting on paired donors with incorporation, and steroid hydroxylase activity.

In the KEGG dataset, the high-risk group showed significant enrichment in pathways related to cell cycle, cytokine-cytokine receptor interaction, ECM-receptor interaction, hematopoietic cell lineage, and neuroactive ligand-receptor interaction. In contrast, the low-risk group was significantly enriched in pathways such as drug metabolism-cytochrome P450, fatty acid metabolism, metabolism of xenobiotics by cytochrome P450, primary bile acid biosynthesis, and retinol metabolism (Supplementary Figure 4A-D). The gene mutation analysis revealed that TP53 had the highest mutation rate in the high-risk group, whereas CTNNB1 exhibited the highest mutation rate in the low-risk group (Supplementary Figure 4E and F).

The Relationship of Risk Score with TIME

Immune cell infiltration analysis showed that the levels of activated DCs, immature DCs, macrophages, T follicular helper cells, and Tregs were significantly higher in the high-risk group compared to the low-risk group. Additionally, antigen-presenting cell co-stimulation and MHC class I activities were more prominent in the high-risk group. Conversely, mast cells, NK cells, and Type II IFN response were more abundant in the low-risk group (Figure 3A and B). Furthermore, the infiltration levels of naive B cells, M0 macrophages, M2 macrophages, CD8⁺ T cells, and T follicular helper cells were significantly correlated with OS (Figure 3C-G). The risk score was also found to be negatively correlated with Stromal scores, Immune scores, and ESTIMATE scores (Figure 3H-J).

Figure 3 TIME analysis. (A) The difference in immune cell infiltration between high-risk and low-risk groups. (B) Differences in immune-related functions between high-risk and low-risk groups. (C-G) Immune cells associated with HCC survival. (H-J) Differences in Stromal scores, Immune scores, and ESTIMATE scores between high-risk and low-risk groups. *: P<0.05; **: P<0.01; ***: P<0.001.

Single-cell analysis revealed that cells from three normal and three HCC tissues in the GSE175793 cohorts were grouped into 17 distinct clusters. These clusters mainly consisted of six cell types: γδ T cells, B cells, myeloid cells, endothelial cells, epithelial cells, and hepatocytes (Supplementary Figure 5A and B). Notably, among the six cell types, the differential expression of PRKAG1 and B3GAT3 between normal and HCC tissues was most pronounced in γδ T cells (Supplementary Figure 5C and D).

Analysis and Validation of IGRPS Immunotherapy

When comparing the efficacy of two immune checkpoint inhibitors (ICIs: anti-CTLA4 and anti-PD-1) between the high-risk and low-risk groups, the results demonstrated that the low-risk group experienced better outcomes with combined anti-CTLA4 and anti-PD-1 treatment compared to the high-risk group. However, there were no significant differences in efficacy between the two groups when treated with anti-CTLA-4 alone or anti-PD-1 alone (Figure 4A-C). Validation using the IMvigor210 cohort further indicated that patients who responded to immunotherapy had a lower risk score than those who did not, suggesting that patients in the low-risk group were more likely to benefit from immunotherapy (Figure 4D).

Figure 4 Analysis of the effect of IGRPS on HCC immunotherapy response. (A-C) Differences in treating the ICIs anti-PD-1 and anti-CTLA-4 in high-risk and low-risk groups. (D) The IMvigor210 cohort was applied to verify the difference in immunotherapy response between high-risk and low-risk groups. (E and F) Relationship between risk score and TMB. (G-J) Relationship between risk score and immune escape characteristics. ***: P<0.001.

Additionally, the results revealed a significant negative correlation between TMB and risk score, with patients in the high-risk group exhibiting lower TMB (Figure 4E and F). Analysis of the relationship between risk score and immune escape showed that T-cell dysfunction, T-cell exclusion, microsatellite instability, and TIDE scores were all higher in the high-risk group than in the low-risk group. This indicated a higher likelihood of immune escape and poorer response to immunotherapy in the high-risk group (Figure 4G-J). Collectively, these findings confirmed the sensitivity and utility of the IGRPS model in evaluating the efficacy of immunotherapy for HCC patients.

Validation of IGRPS Gene Expression

Differential expression analysis revealed that the expression levels of IGRPS genes PRKAG1 and B3GAT3 were higher in the high-risk and tumor groups compared to the low-risk and normal groups, while the expression of ALDH2 showed the opposite trend (Figure 5A and B). Correlation analysis further demonstrated that the expression levels of PRKAG1 and B3GAT3 were significantly positively correlated with the risk score, whereas ALDH2 expression was negatively correlated (Figure 5C).

Figure 5 Validation of IGRPS genes ALDH2, PRKAG1, and B3GAT3 expression. (A and B) Gene expression levels of IGRPS genes in high-low risk groups as well as HCC and its adjacent tissues. (C) Relationship between risk score and expression of IGRPS genes. (D) The diagnostic-ROC curve evaluated the diagnostic value of IGRPS genes. (E) The differential expression of the IGRPS gene at the mRNA level in HCC and adjacent tissues was confirmed by RT-qPCR. (F) The differential expression of the IGRPS gene at the protein level in HCC and adjacent tissues was verified by IHC (Scale bar: 50 μm; 20 X objective). *: P<0.05; ***: P<0.001; ****: P<0.0001.

Moreover, the diagnostic ROC curve analysis showed that the AUC values for PRKAG1, B3GAT3, and ALDH2 were 0.86 (95% CI: 0.89–0.82), 0.80 (95% CI: 0.85–0.76), and 0.62 (95% CI: 0.67–0.57), respectively, indicating their potential diagnostic value (Figure 5D). Finally, the results from RT-qPCR and IHC confirmed that the expressions of PRKAG1 and B3GAT3 at the mRNA and protein levels were upregulated in HCC tissues, while ALDH2 was downregulated, aligning with the bioinformatics findings (Figure 5E and F).

PRKAG1 and B3GAT3 Induce HCC Cell Proliferation and Migration by Enhancing the Glycolytic Capacity of HCC Cells

Depletion of PRKAG1 or B3GAT3 led to a significant reduction in glucose consumption, ATP levels, lactate dehydrogenase activity, lactate levels, and ECAR in HCC cells (Figure 6A-E). Additionally, the proliferation and migration abilities of HCC cells were inhibited in the PRKAG1 or B3GAT3 siRNA interference groups. Interestingly, the addition of lactate to the transfected cells restored their proliferation and migration capacities, which were subsequently reversed by adding 2-DG, an inhibitor of glycolysis (Figure 7A-C).

Figure 6 PRKAG1 and B3GAT3 regulate the glycolytic capacity of HCC cells. The levels of glucose consumption (A), ATP (B), lactate dehydrogenase activity (C), lactate (D), and ECAR (E) were measured and analyzed in HCC cells transfected with NC-siRNA, PRKAG1-siRNA, or B3GAT3-siRNA. ***: P<0.001.

Figure 7 The effect of PRKAG1 and B3GAT3 on glycolytic capacity, proliferation, and migration. (A) CCK-8 assay. (B) Transwell migration assay (Scale bar: 200 μm; 10 X objective). (C) Wound-healing assay (Scale bar: 200 μm; 10 X objective). *: P<0.05; **: P<0.01; ***: P<0.001.

Moreover, IF analysis showed that depletion of PRKAG1 or B3GAT3 resulted in a downregulation of the proliferation marker Ki-67 in the nucleus, further confirming at the protein level that PRKAG1 and B3GAT3 promoted the proliferation of HCC cells (Figure 8A). These findings suggested that PRKAG1 and B3GAT3 might play oncogenic roles in HCC by enhancing the glycolytic capacity of the cancer cells, thereby promoting tumor growth.

Figure 8 The effect of PRKAG1 and B3GAT3 on proliferation at the protein level. (A) IF assay. Blue represented the cell nucleus, and red indicated the localization and expression of the Ki-67 protein. The deeper the red, the higher the expression.

Discussion

HCC is the sixth most common cancer globally and is characterized by a complex mechanism of occurrence and progression, leading to a poor prognosis.¹³ Although numerous prognostic biomarkers have been identified, they often focus on the impact of a single biological process, neglecting the integration with other biological phenomena. For instance, Kong et al have identified five glycolytic genes associated with HCC prognosis and used two glycolytic gene features to predict patient outcomes.¹⁴ Similarly, Chen et al have developed a model based on five prognostic genes related to glycolysis/gluconeogenesis.¹⁵ Other previous studies enrolled mitochondrial-related genes (MRGs), RNA-binding proteins, or cytotoxic T lymphocyte-evasion genes to construct HCC prognostic models, have also considered only one kind of biological marker.^16–18

In this study, glycolysis and TIME were considered two crucial biological markers of tumors. Therefore, we constructed an IGRPS using IGRGs. The IGRPS demonstrated excellent predictive performance and was an independent prognostic factor for HCC. Additionally, it provided insights into the state of the TIME and the changes in tumor glycolysis. Notably, the IGRPS comprised only three genes, ALDH2, PRKAG1, and B3GAT3, making it practical for clinical application. Most importantly, the roles of these model genes in the development of HCC were experimentally validated, underscoring their relevance in tumor biology.

PRKAG1 and B3GAT3 were found to be upregulated in HCC, suggesting that they might serve as risk factors and possess oncogenic effects in the progression of HCC. Conversely, the downregulation of ALDH2 in HCC might help inhibit tumor progression. Notably, the diagnostic ROC curves confirmed that these model genes had excellent diagnostic value, indicating that they were valuable diagnostic and prognostic markers for HCC.

While there has been extensive research on the phenotypic and mechanistic roles of ALDH2 in HCC, studies focusing on PRKAG1 and B3GAT3 in HCC are relatively limited. Consequently, this study centered on PRKAG1 and B3GAT3. PRKAG1 (protein kinase AMP-activated non-catalytic subunit gamma 1) is one of the γ subunits of the AMPK (AMP-activated protein kinase) complex and is involved in AMPK activation. AMPK serves as a critical cellular energy sensor that regulates metabolic pathways to maintain energy homeostasis.¹⁹ Upon activation, AMPK enhances catabolic processes that generate ATP, such as glycolysis.^20,21

In breast cancer, PRKAG1 regulates AMPK activity by sensing the intracellular AMP/ATP ratio. When activated, AMPK can directly phosphorylate and activate 6-phosphofructo-2-kinase (PFK-2), thereby increasing levels of fructose-2,6-bisphosphate and promoting glycolysis.^22,23 Additionally, PRKAG1 plays a significant role in immune regulation. Studies have shown that PRKAG1, as part of the AMPK complex, assists immune cells in switching metabolic pathways to meet energy demands during activation and function.²⁴ The PRKAG1-mediated AMPK signaling pathway is crucial for T cell activation and functional maintenance, influencing the anti-tumor immune response.²⁵ However, no studies have yet examined the role of PRKAG1 in HCC, highlighting the novelty and importance of this research focus.

B3GAT3 (β-1,3-glucuronyltransferase 3) is a glycosyltransferase primarily involved in synthesizing glycosaminoglycans.²⁶ B3GAT3 modulates cellular metabolic states by regulating the interactions of glycosaminoglycans with cell surface receptors, thereby influencing signaling pathways such as PI3K/AKT and MAPK. Additionally, the products and energy levels derived from glycolysis feedback regulate these signaling pathways, creating a complex regulatory network.^27,28

Moreover, studies have identified B3GAT3 as a risk factor in the progression of various tumors. For instance, B3GAT3 expression is upregulated in osteosarcoma stem cells (OSC) and is associated with poor survival rates in osteosarcoma patients.²⁹ Zhao et al have screened nine genes related to amino acid metabolism and identified B3GAT3 as an independent risk factor for HCC.³⁰ Similarly, Wu et al have identified four glycolysis-related genes linked to the prognosis of clear cell renal cell carcinoma (ccRCC), with upregulated expression of B3GAT3 in ccRCC correlating with a poor prognosis for patients.³¹ These findings suggest that B3GAT3 plays a significant role in tumor progression by influencing metabolic and signaling pathways, making it an essential focus for understanding its impact on cancer biology, including HCC.

This study demonstrated, for the first time, that PRKAG1 and B3GAT3 could enhance the glycolytic capacity of HCC cells, leading to the production of a substantial amount of lactate, which in turn induced the proliferation and migration of HCC cells. This finding aligned with the glycolysis profile analysis conducted in this study. Both bioinformatics analyses and experimental validations confirmed that PRKAG1 and B3GAT3 were upregulated in HCC, and they served as strong prognostic indicators for the disease. These results suggested that PRKAG1 and B3GAT3 could act as risk factors for HCC and play oncogenic roles.

TMB quantifies the total number of genetic mutations within a tumor. A high TMB indicates a greater frequency of mutations, which may lead to the production of more neoantigens recognizable by the immune system, thereby potentially enhancing the immune response.³² Consequently, tumors with high TMB are more likely to respond to ICIs, such as anti-PD-1/PD-L1 therapies.³³ In this study, TMB was found to be negatively correlated with the risk score; the lower TMB observed in the high-risk group suggested a poorer response to immunotherapy.

The TIDE score represents a measure of immune escape potential in tumors. It predicts the efficacy of tumor immunotherapy and patient survival based on two primary mechanisms: immune antigen presentation and T-cell exhaustion in the TIME. A lower TIDE score indicates that patients are more likely to respond effectively to immunotherapy.³⁴ In this study, the expression levels of PRKAG1 and B3GAT3, as well as the TIDE score were significantly lower in the low-risk group compared to the high-risk group, indicating a reduced risk of immune escape potential and a better response to immunotherapy. This may be due to PRKAG1 promotes glycolysis by activating AMPK, leading to increased lactate production. The accumulation of lactate lowers the extracellular pH, which inhibits the cytotoxicity of T cells and NK cells, thereby enhancing immune suppression and contributing to immune escape.³⁵ Similarly, B3GAT3 promotes the synthesis of glycosaminoglycans (GAGs) which can interact with immune checkpoint molecules such as PD-L1 to suppress T cell function, thereby facilitating tumor cell immune escape.³⁶

Furthermore, the immunotherapy analysis results indicated that the combination of anti-CTLA-4 and anti-PD-1 treatments was more effective in the low-risk group compared to the high-risk group. However, there was no significant difference in efficacy between the two groups when anti-CTLA-4 or anti-PD-1 was administered as monotherapy. Combination therapy offers complementary anti-tumor effects by targeting different immune escape mechanisms.³⁷ For instance, anti-CTLA-4 therapy enhances T cell cytotoxicity by activating the initial T cell response, while anti-PD-1 therapy improves T cell effectiveness by releasing the inhibition on effector T cells.³⁸ Meanwhile, patients in the low-risk group exhibited higher TMB and lower immune escape characteristics, which allowed combination therapy to more effectively stimulate a multifaceted immune response. Furthermore, TIME analysis revealed that the infiltration scores of mast cells and NK cells were higher in the low-risk group, while the scores of macrophages and Tregs were lower. Studies have shown that mast cells can release pro-inflammatory mediators to activate other immune cells (such as DCs and T cells) or interact with NK cells, thereby enhancing anti-tumor immune responses.³⁹ In contrast, M2-type macrophages and Tregs promote tumor growth and immune escape by secreting immunosuppressive factors such as IL-10 and TGF-β.^40,41 These findings suggested that patients with a low-risk profile might benefit more from combination immunotherapy due to their more favorable TIME.

This study systematically explored the application value of the IGRPS in HCC, offering new perspectives for the diagnosis, prognosis, and immunotherapy of HCC. However, there are some limitations to this study. First, the data were derived from RNA-seq and clinical information available in public databases such as TCGA and GEO. These samples come from diverse sources and may include data generated from different laboratories and platforms, leading to heterogeneity and batch effects that could impact the consistency and generalizability of the model. Therefore, the predictive value of IGRPS needs further validation in large, prospective clinical cohorts to ensure its robustness and reliability.

Conclusions

This study introduced novel prognostic features for HCC by integrating immune status and glycolytic processes, significantly enhancing the ability to predict patient outcomes and response to immunotherapy. These features offered a promising framework for personalized therapeutic strategies in HCC, allowing for more tailored interventions that align with individual patient profiles. Moreover, these insights provided valuable therapeutic targets for advancing HCC immunotherapy, paving the way for more effective and precise treatment approaches in clinical practice.