Development and Validation of a Novel Nomogram Integrated with Hypoxic and Lactate Metabolic Characteristics for Prognosis Prediction in Hepatocellular Carcinoma

Introduction

Primary liver cancer is still the sixth most common cancer and the third leading cause of cancer-related death worldwide, with more than 1 million new cases expected annually by 2025.^1,2 Hepatocellular carcinoma (HCC) is the most prevalent pathological type of liver cancer, accounting for about 90% of cases. In the recent years, significant advancements have been made in the systematic treatment of HCC, including tyrosine kinase inhibitors (TKIs), immune-checkpoint inhibitors (ICIs) and monoclonal antibodies.³ Lenvatinib, a newly approved drug, was found to be non-inferior to the conventional TKI sorafenib in improving the overall survival (OS) of HCC patients. Additionally, the combined blockade of programmed cell death 1 ligand 1 (PD-L1) and vascular endothelial growth factor (VEGF) has shown remarkable efficacy in improving the prognosis of patients with unresectable HCC.^4,5 However, chemotherapy and immunotherapy resistance still negatively impact the prognosis of many HCC patients.^6–8 Therefore, it is imperative to identify a signature that can predict the survival, immune microenvironment, and treatment responses of HCC patients.

Hypoxia is a prominent hallmark in HCC, and plays a critical role in tumor progression and therapy resistance.^9–11 Studies have revealed that hypoxia can endow HCC cells with sorafenib resistance by upregulating miR-210-5P, thus reducing the expression of ATPase family AAA domain containing 3A and promoting mitophagy.¹² Moreover, hypoxia-activated pleomorphic adenoma gene-like-2/epidermal growth factor receptor (EGFR) signaling can impair the sensitivity of HCC cells to the anti-EGFR drug erlotinib.¹³ In hypoxia-high tumor regions, the interaction between regulatory T cells (Tregs) and type 2 conventional dendritic cells is enhanced, causing an immunosuppressive microenvironment and a more proliferative and invasive tumor phenotype.¹⁴ In addition, hypoxia can reprogram glucose metabolism from oxidation to glycolysis through promoting Yes-associated protein-mediated stabilization of hypoxia-inducible factor 1α (HIF-1α) protein, thus providing sufficient energy for the rapid growth of HCC.¹⁵

Lactate is the main product of anaerobic glycolysis. Increasing evidence suggests that intra-tumoral lactate is a critical regulator that induces the immunosuppressive microenvironment and promotes tumor growth rather than being just a mere metabolic waste.^16–19 Lactate confers HCC cells with ferroptosis resistance by elevating the levels of anti-ferroptosis monounsaturated fatty acids and contributes to enhanced proliferation of tumor cells.²⁰ Intracellular lactate facilitates the proliferation, self-renewal and migration of HCC stem cells by driving the lactylation of H3 histone at K9 and K56.²¹ Furthermore, the high level of MOESIN lactylation in Tregs correlates with a poor response to anti-programmed cell death-1 (PD-1) therapy in HCC. The underlying mechanism is that the lactylation of MOESIN at K92 activated transforming growth factor β (TGF-β)/SMAD3 signaling, promoting the proliferation and function of Tregs.²²

Therefore, this study aims to identify a hypoxia- and lactate metabolism-related prognostic model (HLPM) to evaluate the survival, immune landscape and treatment responses for patients with HCC, and develop a novel nomogram integrated with HLPM and clinical characteristics to predict the prognosis in HCC.

Materials and Methods

Data Collection

RNA expression and clinical data of 350 HCC cases with OS longer than 30 days and 49 normal cases were downloaded from The Cancer Genome Atlas (TCGA) database (https://portal.gdc.cancer.gov/). The TCGA-LIHC cohort was regarded as a training set. In order to determine the stability of the prognostic model, ICGC-LIRI-JP cohort (n = 230) with OS longer than 30 days obtained from HCCDB (https://lifeome.net/database/hccdb/home.html) was served as a validation set. Genotype-tissue expression (GTEx) provided by the University of California Santa Cruz Xena (https://xenabrowser.net/datapages/) was used to obtain gene expression of normal liver tissues (n = 111). The hypoxia-related genes (HRGs) (n = 493) and lactate metabolism-related genes (LRGs) (n = 355) were downloaded from the Molecular Signatures Database (https://www.gsea-msigdb.org/gsea/msigdb/index.jsp).²³

Single-cell RNA sequencing (scRNA-seq) dataset GSE149614 was collected from the GEO database, including 10 primary HCC tissues, 2 portal vein tumor thrombus tissues, 1 metastatic lymph node tissue and 8 normal liver tissues.²⁴ To quality control, following low-quality cells were filtered out, including (1) unique molecular identifiers (UMIs) ≤500, (2) detected gene expression ≤500 or ≥8000, (3) mitochondrial gene expression ≥30%. Therefore, a total of 25,712 genes and 71,910 cells were generated and included in this study. The data were normalized by “LogNormalized” method and used to find 2000 highly variable genes by “FindVariableFeatures”. To reduce the dimensionality of the data, the “RunPCA” function was applied specifically to the set of 2000 highly variable genes. Cell clustering was performed using the “KNN” method with a resolution parameter set to 0.5. A two-dimensional representation of the clustered cells was obtained using the uniform manifold approximation and projection (UMAP) technique. To identify the categories of cells, a combination of information from previous studies and the “SingleR” package was employed.^24,25 Differential expression analysis among the identified clusters was conducted using the Wilcoxon-Mann–Whitney test and the “FindAllMarkers” function from the “Seurat” package. Genes with |log2 fold change (FC)| >2 and an adjusted p-value <0.01 were considered as significant cellular marker genes.

Differentially Expressed HRGs and LRGs

The R package “limma” was used to identify differentially expressed genes (DEGs) with the threshold for |log2 FC| >1.0 and Benjamini–Hochberg adjusted p-value <0.05. Differentially expressed HRGs (DEHRGs) and LRGs (DELRGs) were defined as the intersection of DEGs and HRGs and LRGs. The protein–protein interaction (PPI) network between DEHRGs and DELRGs (HLRGs) in HCC was constructed using the STRING database (https://string-db.org/).

Construction and Validation of HLPM

First, the univariate Cox regression analysis was used for identifying DEHRGs and DELRGs with a significant prognostic value. To avoid overfitting, R package “glmnet” was used for least absolute shrinkage and selection operator (LASSO) regression. Finally, stepwise backward Cox regression analysis of five selected genes was performed to construct a signature of HRGs and LRGs in TCGA-LIHC cohort []. HCC patients were divided into the high-risk and low-risk groups based on the median risk score. Kaplan–Meier (KM) survival analysis, receiver operating characteristic (ROC) curve, distribution of risk score, association between the survival time and survival status of patients, and the expression of the HLRGs were carried out in training set and validation set to assess the performance of the signature in predicting clinical outcomes.

In addition, Kaplan–Meier survival analysis was performed for subgroups of patients with different clinical characteristics to evaluate the applicability of the signature. Immunohistochemical images of prognostic DEGs and HRGs in HCC and normal liver tissues were downloaded from HPA (https://www.proteinatlas.org/).

Nomogram Development for Predicting Survival

A nomogram was constructed based on the risk score and other clinical indicators to calculate the 1-, 3-, and 5-year OS probabilities of HCC patients. ROC curves, C-index, and calibration curves were used to evaluate the performance of the nomogram. In addition, the net benefit of the nomogram and clinical features alone was evaluated by decision curve analysis (DCA).

Functional Enrichment and Somatic Mutations Analysis

To analyze the critical functional features, we analyzed the DEGs between the two groups with Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) using the R package “clusterProfiler”. Gene set enrichment analysis (GSEA) was also used to explore the important functional pathway between the high-risk and low-risk groups (http://www.gsea-msigdb.org/gsea/index.jsp) by R package “clusterProfiler”. (p-value <0.01 and false discovery rate (FDR) <0.20). In GSEA, p-value <0.01 and FDR <0.20 were considered statistically significant. Somatic mutation analysis was conducted by R package “maftools”.

Assessment of the TME, Immune Cell Infiltration, and Immunotherapy Response

The immune cell abundance identifier (ImmuCellAI) was utilized to assess the infiltrated fractions of immune cells.²⁶ The correlation analysis between the five HLRGs and the filtration of immune cells was also performed. The levels of crucial immune checkpoint genes (ICGs) and chemokines and their receptors were compared between two risk groups. The relationship between the risk score and ICGs expression was determined by Spearman correlation analyses. The online tool named TIDE (http://tide.dfci.harvard.edu/) was used to evaluate the response of patients to immunotherapy.²⁷ We further used the GSE78220 cohort to assess the capability of the HLPM in predicting the efficacy of immunotherapy.²⁸ Patients with complete response (CR) or partial response (PR) were classified as responders and patients with progressive disease (PD) were regarded as non-responders.

Assessment of the Sensitivity to Chemotherapy

The R package “oncoPredict” and the Genomics of Drug Sensitivity in Cancer (GDSC) database (https://www.cancerrxgene.org/) were utilized to calculate the half-maximal inhibitory concentrations (IC50) of potential drugs between the high-risk and low-risk groups and the correlation between IC50 and the expression of the five HLRGs.^29,30 Wilcoxon signed-rank tests were used to compare IC50 values between the two risk groups.

Tissue Samples and qRT-PCR

A total of 14 pairs HCC tissue samples and adjacent normal liver tissue samples were collected from the Shulan (Hangzhou) Hospital and Hangzhou First People’s Hospital with the approval by the Medical Ethics Committee of the hospitals. All patients included in the current study have provided written informed consent. Fresh tumor and non-tumor tissues were preserved in liquid nitrogen to maintain their integrity. The extraction of RNA was realized by using TRIzol (Invitrogen, California, USA). The cDNA was synthesized using the Hifair^® V one-step RT-gDNA digestion SuperMix for qPCR (Yeasen, Shanghai, China). Real-time quantitative PCR analysis was performed on a CFX96 Touch real-time PCR system (Bio-Rad, CA, USA) with qPCR SYBR Green Master Mix (Yeasen, Shanghai, China). The amount of target cDNA was analysed through the conversion of the threshold cycle with reference to β-actin. Primers for the genes are displayed in Table S1.

Statistical Analysis

Statistical analyses were performed by the R software (Version 4.0.5). The unpaired Student’s t-test and F-test were used for comparing normal data, and Wilcoxon test was used for comparing nonnormal data. The Log rank test was used for Kaplan–Meier survival analysis. The clinical characteristics among groups were assessed by χ² test or Fisher’s exact test. A p-value of less than 0.05 was considered as statistically significant.

Results

Identification and Functional Enrichment of DEHRGs and DELRGs

A total of 1759 upregulated genes and 2054 downregulated genes were identified and shown in the volcano plot (Figure 1A). The heat map displayed the expression of the top 50 genes with the most significant difference (Figure 1B). DEGs were intersected with HRGs and LRGs, respectively, resulting in the identification of 115 DEHRGs and 42 DELRGs (Figure 1C). The results of PPI analysis suggested the intrinsic association between DEHRGs and DELRGs (Figure 1D). According to GO analysis, the 115 DEHRGs and 42 DELRGs mainly performed a function in response to hypoxia and decreased oxygen levels (Figure 1E). Based on KEGG analysis, DEHRGs and DELRGs were enriched in HIF-1 and PI3K/Akt signaling pathways (Figure 1F). These findings revealed the potential molecular mechanisms of DEHRGs and DELRGs that involved in the progression of HCC.

Figure 1 Identification and functional enrichment of differential expressed hypoxia-related genes (HRGs) and lactate metabolism-related genes (LRGs) in HCC. (A) The volcano plot of the differentially expressed genes (DEGs) between normal and tumor tissues. (B) The heat map showed the detailed expression of the top 50 DEGs. (C) The Veen diagram of the DEGs, HRGs and LRGs. (D) The protein–protein interaction (PPI) network of differential expressed HRGs and LRGs. GO enrichment analysis (E) and KEGG pathway analysis (F) of differential expressed HRGs and LRGs.

Identification and Validation of the Five HLRGs

To identify DEHRGs and DELRGs correlated with prognosis, we performed univariate Cox regression analysis, identifying 38 DEHRGs and DELRGs with significance (Figure 2A). After performing LASSO regression analysis, only 14 DEHRGs and DELRGs were subjected to backward stepwise multivariate Cox regression analysis (Figure 2B and C). Finally, we selected five genes (KIF20A, IRAK1, ADM, PPARGC1A and EPO) to construct the HLPM [risk score = (0.1870 × KIF20A) + (0.3558 × IRAK1) + (0.2233 × ADM) + (−0.0720 × PPARGC1A) + (0.0639 × EPO)] (Figure 2D). In TCGA and GTEx databases, all five biomarkers were significantly differentially expressed between HCC and normal samples, implying the critical role in the tumor progression and patients’ prognosis. As shown in Figure 2E, levels of KIF20A and IRAK1 were upregulated in HCC tissues, while increased levels of ADM, PPARGC1A and EPO were noted in normal tissues. Consistent with the bulk RNA-seq results, qRT-PCR analysis of clinical tissue samples further corroborated the upregulation of KIF20A and IRAK1, along with a downregulation of ADM, PPARGC1A, and EPO in HCC tissues compared to adjacent normal liver tissues (Figure 2F). Furthermore, survival analysis in HCC patients revealed that high expression of EPO, KIF20A and IRAK1 significantly correlated with poor OS, while high expression of PPARGC1A was related to better OS (Figure 2G). HCC is a kind of heterogeneous tumor that has distinct gene expression profiles in different stages.³¹ As shown in Figure 2H, the mRNA levels of EPO, KIF20A, ADM and IRAK1 were significantly differentially expressed among patients at different stages, reflecting that these biomarkers exerted important effects in the tumor progression and could be regarded as the potential therapeutic targets for HCC. Immunohistochemical images from the HPA database further validated the increased levels of KIF20A and IRAK1 and decreased levels of ADM and EPO in tumor tissues compared with normal tissues (Figure 2I). Immunohistochemical staining for PPARGC1A in tumor and normal tissues was not provided in the HPA database. However, several studies have indicated that the protein level of PPARGC1A was downregulated in the human HCC tissues compared with the normal tissues.^32,33 These results suggest the prognostic value of the five DEHRGs and DELRGs in HCC.

Figure 2 Identification and validation of the five differentially expressed hypoxia- and lactate metabolism-related genes (HLRGs). (A) Univariate Cox regression analysis of the differentially expressed hypoxia-related genes (DEHRGs) and differentially expressed lactate metabolism-related genes (DELRGs). (B) Tuning parameter (λ) selection in LASSO model. (C) The LASSO coefficient profile of 14 prognostic DEHRGs and DELRGs. (D) The results of stepwise multivariate Cox regression analysis of HLRGs were shown by forest graph. (E) Differential expression of the five HLRGs between the HCC and normal tissues in TCGA and GTEx databases. (F) The expression levels of the five HLRGs in 14 pairs HCC and matched adjacent normal tissues were examined by qRT-PCR. (G) Kaplan-Meier survival analysis of patients with different expression of five selected HLRGs. (H) The expression levels of five HLRGs in patients at different stages. (I) Immunohistochemical staining for four HLRGs in the normal liver and HCC tissues. *p < 0.05; **p < 0.01; ****p < 0.0001.

To further explore and validate the action of the five prognostic genes, we performed analysis using the scRNA-seq data. We finally obtained 7 cell clusters and annotated the cell type of each cluster, including endothelial cells, fibroblasts, hepatocytes, myeloid cells, T/NK cells, B cells and tumor cells (Figure S1A and B). As shown in Figure S1C, most of fibroblasts and B cells came from the tumor samples, indicating their critical roles in the development of HCC. Consistent with the bulk RNA-seq and qRT-PCR analyses, KIF20A and IRAK1 were primarily expressed in cells from tumor tissues, implying that they exerted vital effects for tumor progression. However, the expression levels of ADM, EPO and PPARGC1A were increased in tumor tissues compared with normal tissues, which is inconsistent with the results from the bulk RNA-seq and qRT-PCR analyses. ADM was mainly expressed in myeloid and endothelial cells, while PPARGC1A was mainly expressed in hepatocytes. The expression of EPO was detected in only a small proportion of cells in both tumor and normal samples (Figure S1C).

Evaluation and Validation of the HLPM

In the training and validation cohorts, patients were divided into high- and low-risk groups based on the median value of risk scores. KM analysis revealed that high-risk patients with HCC had poor OS (Figure 3A and B). In the TCGA-LIHC cohort, HLPM had area under curve (AUC) values of 0.677, 0.726, and 0.716 for predicting 1-, 3-, and 5-year OS, respectively (Figure 3C). In the ICGC-LIRI-JP cohort, the 1-, 3-, and 5-year AUC values were 0.850, 0.747, and 0.606, respectively (Figure 3D). The risk score and survival status plots suggested the negative correlation between the risk score and survival time. Furthermore, the risk heat map revealed the differential expression of the five genes in two risk groups (Figure 3E and F). Overall, these results indicate the accuracy and stability of the HLPM.

Figure 3 Evaluation and validation of the hypoxia- and lactate metabolism-related prognostic model (HLPM). (A and B) Kaplan-Meier survival analysis of two risk groups in the TCGA-LIHC cohort and ICGC-LIRI-JP cohort. (C and D) Receiver operating characteristic (ROC) curve verified the prognostic value of the HLPM in the TCGA-LIHC cohort and ICGC-LIRI-JP cohort. (E and F) Distribution of risk score, survival time and survival status of patients, and expression of the five HLRGs in the TCGA-LIHC cohort and ICGC-LIRI-JP cohort.

To further verify the universal applicability of the HLPM, we performed a stratification analysis. Patients with HCC were divided to different subgroups based on gender (male vs female), age (≥60 vs <60 years), BMI (≥25 vs <25), AFP level (≥100 vs <100 ng/mL), histological grade (grade 1–2 vs grade 3–4), American Joint Committee on Cancer (AJCC) T stage (T1-2 vs T3-4), AJCC N stage (N0 vs Nx), and AJCC M stage (M0 vs Mx). The KM survival analysis showed that patients with high-risk scores were correlated with poor outcomes in many subgroups (Figure S2A–P).

Functional Enrichment and Mutation Analysis of the Two Risk Groups

A total of 255 upregulated genes and 321 downregulated genes were identified in the high-risk group compared with the low-risk group (Figure 4A). The heat map in Figure 4B shows the expression of the top 50 DEGs with the most significant difference. GSEA analysis based on the GO database indicated that cell cycle DNA replication, chromosome segregation and cytosolic ribosome were enriched in the high-risk group (Figure 4C). Conversely, the low-risk group showed enrichment in signatures related to amino acid catabolic process, fatty acid catabolic process and fatty acid oxidation (Figure 4D). KEGG-based GSEA analysis indicated the proliferation-related pathways were enriched in the high-risk group when compared with the low-risk group (Figure 4E and F). As shown in Figure 4G, no strong correlation between the risk score and tumor mutation burden (TMB). The waterfall plot in Figure 4H and I exhibited the top 10 genes with the highest frequency of mutation in both risk groups. TP53 stood out as the most commonly mutated gene in the high-risk group, while CTNNB1 emerged as the most predominant mutated gene in the low-risk group.

Figure 4 Functional enrichment and mutation analysis of two risk groups in TCGA-LIHC cohort. (A) The volcano plot of the differentially expressed genes (DEGs) between high- and low-risk groups. (B) The heat map showed the detailed expression of the top 50 DEGs. GO-based GSEA enrichment analysis in high-risk (C) and low-risk (D) groups. The pathways enriched in the high-risk (E) and low-risk (F) groups according to GSEA analysis in KEGG database. (G) The relationship between the risk score and tumor mutation burden (TMB). Mutation profile in high-risk (H) and low-risk (I) groups.

Different Immune Landscapes of the Two Risk Groups

To investigate the relation between HLPM and tumor immune microenvironment, we analyzed the expression of immune checkpoint genes (ICGs), chemokines and receptors, as well as the infiltration fractions of immune cells in two risk groups. The results showed that the expression of ICGs, including HAVCR2, TIGIT, PDCD1, LAG3 and CTLA4 were upregulated in the high-risk group (Figure 5A), indicating the immunosuppressive tumor microenvironment in patients with high-risk scores. However, their expression levels only had a very weak positive correlation with the risk scores, which could be caused by the heterogeneity of samples. It is worth noting that HAVCR2 was the most significantly elevated in tumor tissues and had the strongest association with risk scores, while TIGIT was the least elevated and had the weakest correlation with risk scores among the five widely known ICGs (Figure 5B). Key chemokines [eg, C-C motif chemokine ligand 20 (CCL20), C-X-C motif chemokine ligand 8 (CXCL8) and CCL3], which are associated with immunosuppression and tumor progression, were also increased in high-risk group (Figure 5C). By performing ImmuCellAI analysis, we found that the infraction of exhausted immune cells and Tregs in the high-risk group was significantly higher than that in the low-risk group (Figure 5D). CCL20 has been reported to recruit the intratumoral infiltration of Tregs,¹⁴ suggesting the reliability of the HLPM and related analysis. Moreover, the results of correlation analysis showed that ADM, EPO, IRAK1, KIF20A and PPARGC1A significantly affected the infiltration of immune cells in HCC, indicating the predictive values of the five biomarkers in the tumor immune microenvironment (Figure 5E). Notably, KIF20A which was the most closely associated with the infiltration of immune cells, had the strongest negative and positive correlation with mucosal-associated invariant T (MAIT) cells infiltration and B cells infiltration respectively. Decreased MAIT cells infiltration has been reported to be associated with poor prognosis of HCC patients,^34,35 which consistent with our results. Infiltrating B cells were reported to have dual roles in the progression of HCC,^36,37 and our results indicated that ADM-associated B cells would primarily exert effects on promoting HCC. These results indicate the good performance of the HLPM in predicting the immune functions in tumor microenvironment.

Figure 5 Different immune landscapes of the two risk groups in TCGA-LIHC cohort. (A) The expression of immune checkpoint genes (ICGs) (HAVCR2, TIGIT, PDCD1, LAG3 and CTLA4) in high- and low-risk groups. (B) The correlation between levels of ICGs and risk score. (C) The expression of chemokines and receptors in two risk groups. (D) The infiltration fractions of immune cells in two risk groups. (E) The correlation between the infiltration of immune cells and the five HLRGs. TPM, transcripts per million; *p < 0.05; **p < 0.01; ***p < 0.001.

Prediction of Response to Immunotherapy Between the Two Risk Groups

Next, we evaluated the response of patients to immunotherapy. As shown in Figure 6A, patients in the high-risk group had higher TIDE scores, suggesting that their tumor have the stronger ability to escape immune surveillance. In addition, risk scores were noted to be higher in non-responders compared to responders (Figure 6B). Consistently, the proportion of patients with response to immunotherapy was higher in the low-risk group than that in the high-risk group (39% vs 30%) (Figure 6C). These results demonstrate that patients with high-risk scores are prone to the immunosuppressive tumor microenvironment and thus less sensitive to immunotherapy.

Figure 6 The assessment of the response to immunotherapy. (A) Differences in TIDE scores between high- and low-risk groups in TCGA-LIHC cohort. (B) Differences in risk scores between the responders and non-responders in TCGA-LIHC cohort. (C) The proportion of patients with response to immunotherapy in two risk groups in TCGA-LIHC cohort. (D) Kaplan-Meier analysis of two risk groups in GSE78220 cohort. (E) Comparison of the risk scores between the responders and non-responders in GSE78220 cohort. (F) The proportion of patients with response to anti-PD-1 immunotherapy in two risk groups in GSE78220 cohort.

Furthermore, we tested the capability of the prognostic model in a common real-world immunotherapy cohort (anti-PD-1 in the GSE78220). The KM survival analysis showed no difference in OS between the two risk groups (Figure 6D). However, the risk scores in non-responders were significantly higher than those in responders (Figure 6E). We also found that the low-risk group had a better CR/PR rate compared to the high-risk group (64% vs 38%), indicating a negative correlation between the efficacy of immunotherapy and risk scores (Figure 6F).

Prediction of Response to Potential Chemotherapeutic Agents

To identify the potential chemotherapeutic agents for HCC, we evaluated IC50 values of the drugs in high- and low-risk groups. Our results showed that patients in the low-risk group were more sensitive to GNE-317, NVP-ADW72, savolitinib, foretinib, BIBR-1532, pyridostatin, BMS-754807, LJI308, AZ6102, AZD6482, AT13148 and JQ1, while patients with a high-risk score were more sensitive to WIIKI4, sepantronium bromide and WEHI-539 (Figure S3A–O). Sepantronium bromide has previously been shown to be more sensitive to TP53-mutant HCC patients with poor prognosis in another study,³⁸ which is consistent with our results.

Additionally, we found that KIF20A was positively correlated with the IC50 values of foretinib, LJI308, AZ6102, AZD6482, AT13148 and BMS-754807, but negatively correlated with the IC50 value of WEHI-539. ADM showed a negative association with the IC50 values for foretinib and pyridostatin. IRAK1 was positively associated with the IC50 values of AZ6102, AT13148, foretinib and BIBR-1532. PPARGC1A exhibited a positive correlation with the IC50 of WEHI-539, whereas it displayed a negative correlation with AZD6482 (Figure S3P). Therefore, the HLPM has the potential to guide patients with HCC in receiving more appropriate drug treatment.

Establishment and Evaluation of a Predictive Nomogram

Univariate and multivariate Cox regression analyses in the TCGA-LIHC cohort were conducted to identify independent factors of the prognosis, including risk scores and other clinical characteristics. The results indicated that age, grade and risk score were significantly correlated with OS and had independent prognostic values (Figure 7A and B). Subsequently, the nomogram was constructed to quantitatively predict the 1-, 3-, and 5-year OS of patients with HCC (Figure 7C). The prediction model had a C-index value of 0.734. The calibration curves indicated good performance of the nomogram in predicting patients’ prognosis (Figure 7D). The AUC of 1-, 3-, and 5-year OS predictions were 0.743, 0.784 and 0.770, respectively (Figure 7E). The DCA curves for 1- and 3-year showed that the nomogram could provide more net benefits than clinical characteristics alone (Figure 7F and G).

Figure 7 Construction and evaluation of nomogram in TCGA-LIHC cohort. The univariate (A) and multivariate (B) Cox regression analysis for risk score and other clinic characteristics. (C) The nomogram for predicting 1-, 3-, and 5-years overall survival (OS) by risk score, age and grade. (D) Calibration curves of the nomogram to predict 1-, 3- and 5-year OS probabilities. (E) Receiver operating characteristic (ROC) curves of the nomogram. Decision curve analysis (DCA) of the nomogram and other clinical characteristics at 1- (F) and 3-year (G) indicated the net clinical benefits.

Discussion

Hypoxia and lactate are crucial regulators in the development and progression of various cancers. Hypoxia has been reported to reshape tumor immune microenvironment through modulating multiple signaling pathways (eg, HIF, STAT and NF-κB), thus promoting tumor growth.³⁹ Accumulating evidence has also indicated that lactate plays an important role in cell metabolism and immune function.^40,41 Lactylation of adenylate kinase 2 at K28 impairs its function and promotes the progression of tumor cells.⁴² Moreover, hypoxia and lactate contribute to the resistance of therapies, resulting poor prognosis of patients with cancer. HIFs such as HIF-1α and HIF-2α have been reported to endow HCC cells with sorafenib resistance.⁴³ Lactate increases the resistance of gastric cancer cells to anlotinib through regulating BDNF-TrkB pathway, enhancing the interaction between fibroblasts and cancer cells.⁴⁴ Establishing a hypoxia- and lactate metabolism-related signature is of great significance for predicting prognosis and evaluating therapeutic response in cancer patients.

HCC is a type of solid tumor that characterized by hypoxia and lactate overproduction. Moreover, it has been reported that hypoxia and lactate were intimately associated with therapy resistance and poor prognosis of HCC patients.^17,45 Prognostic models based on HRGs or LRGs have been constructed for HCC.^46,47 However, no research has been published on a HLPM for HCC, which integrated with the characteristics of both hypoxia and lactate metabolism. Given the heterogeneity of HCC and the strong relationship between hypoxia and lactate metabolism, it remains necessary to construct prognostic models of HLRGs for predicting survival, immune landscape and treatment responses of patients with HCC. To construct the HLPM, we identified five prognostic HLRGs including ADM, EPO, IRAK1, KIF20A and PPARGC1A, whose expression levels and distributions were demonstrated by the bulk RNA-seq, qRT-PCR and scRNA-seq analyses. However, inconsistent with the results from scRNA-seq analysis, the expression levels of ADM, EPO and PPARGC1A were downregulated in tumor tissues compared with normal tissues in the bulk RNA-seq and qRT-PCR analyses, which could be attributed to the heterogeneity among clinical samples and different detection depths. Multidimensional analysis verified the stability and reliability of the HLPM in predicting the prognosis of patients. The GSEA results indicated a significant enrichment of genes related to cell cycle and DNA replication in the high-risk group compared to the low-risk group, potentially resulting in the poor efficacy and survival of the high-risk patients. Although there was no strong correlation between risk scores and TMB, the mutation profiles were different between two risk groups.

Tumor immune microenvironment is significantly correlated with tumor progression and the response to immunotherapy.^48–50 Our study found that the levels of critical ICGs were significantly upregulated in the high-risk group compared to the low-risk group and were positively correlated with risk scores. TIDE score was also noted to be increased in these patients. In the recent years, immunotherapy has made significant advances in the management of HCC, and one of the most promising therapies is the combination of atezolizumab and bevacizumab.^6,51,52 However, the response of each patient to immunotherapy varies, necessitating the need for a signature that can predict efficacy.^53,54 In this study, the HLPM was applied into the prediction of response to immunotherapy, and the results indicated that the negative association between risk scores and immunotherapy response. The HLPM was validated in the immunotherapy cohort GSE78220, showing good performance in predicting immunotherapy efficacy.

Our study also analyzed the sensitivity of potential drugs in treating HCC. We found that IC50 values of major chemical agents were higher in the high-risk group, while some drugs such as WIKI4, sepantronium bromide and WEHI-539 were more sensitive to patients at high risk. Correlation analysis further illustrated the close association between the five HLRGs and the sensitivity of potential chemical agents. Our findings provide new options for personalized treatment of HCC.

Nomograms have been widely used in predicting survival of patients through computing scores of different risk factors. In order to make the prognostic model more targeted and more comprehensive for clinical practice, we constructed a nomogram based on the risk scores of HLPM and clinical characteristics. The novel nomogram integrated with age, histological grade and risk scores of HLPM exhibited a better capacity in prognosis prediction and could provide more net benefits than clinical characteristics or HLPM alone. Our results would offer a novel option for personalized decision-making in HCC management.

However, our study has several limitations. First, the expression profiles of HCC patients were obtained from public databases such as TCGA and ICGC, which may introduce selection bias. Therefore, further prospective studies with large cohorts are needed to validate the reliability of this prognostic signature. In addition, the underlying mechanisms by which the five HLRGs affect tumor progression and therapeutic efficacy need to be demonstrated in cell and animal experiments.

Conclusion

In summary, our study identified five HLRGs that exhibited favorable prediction performance. Using these HLRGs, we constructed a prognostic model for HCC, which can reliably predict the clinical outcome and therapeutic responses of patients with HCC. Finally, we developed a nomogram integrated with age, histological grade, and risk score of HLPM, presenting a novel option for personalized decision-making in HCC management.