Stemness-hypoxia genes PPARGC1A, PFKFB3, and SLC2A5 are associated with prognosis and tumor immune microenvironment in hepatocellular carcinoma

Introduction

As the sixth most frequent cancer, liver cancer claims over 800,000 lives annually worldwide (1). Particularly, hepatocellular carcinoma (HCC) is a major type of liver malignant tumors and one of the deadliest cancers worldwide. Considerable success has been achieved for HCC treatment in the past decade. However, only 40% patients are diagnosed during the early stages, and the survival time of HCC patients in the advanced stage is usually <1 year (2,3). Although potentially curative treatments are available for patients with early-stage HCC, the recurrence rate remains high (4,5).

Cancer stem cells (CSCs) are one of the major contributors to the recurrence of cancers (6,7). The CSCs show unlimited growth, invasive capacity, and called cancer stemness (8). The phenotype, appearance and function of CSCs are changed compared with those of normal tumors. Such changes are generally induced by aging tumor cells, radiation, and chemotherapy. These triggers can directly contribute to cancer stemness by increasing the plastic phenotype of CSCs, activating the stemness pathway in non-CSCs and promoting aging escape and subsequent activation of the stemness pathway (9). HCC is a highly heterogeneous cancer and recent evidence reveals the existence of CSCs in HCC (10,11). The stemness of HCC is closely associated with resistance to chemotherapy, distant metastasis, and cancer recurrence (11). Therefore, the evaluation of cancer stemness holds great potential in facilitating the clinical management of HCC. Several markers of HCC stemness are currently available, and new markers are emerging to advance this field (11). Recently, a machine learning-based bioinformatics algorithm was proposed for the evaluation of cancer stemness using transcriptomic and epigenetic data from The Cancer Genome Atlas (TCGA) or other databases (12).

Hypoxia is prevalent in solid tumors (13,14), and it is a critical regulator of stemness in multiple cancers (15-18). Several studies suggest that hypoxia drives cancer stemness through a hypoxia-inducible factor 1α (HIF-1α)-dependent manner in HCC (19,20). Hypoxia is a potent predictor of cancer prognosis in human tumors, including HCC (21-25). Additionally, some hypoxia-related genes are involved in cancer development. For instance, peroxisome proliferator-activated receptor-gamma coactivator 1 alpha (PPARGC1A) has been reported to be hypoxia-responsive and contribute to the progression and prognosis of various cancers, including HCC (26). The glycolytic activator 6-phosphofructo-2-kinase/fructose 2,6-bisphosphatase 3 (PFKFB3) can be O-GlcNAcylated in response to hypoxia and plays a key role in promoting cell cycle progression and tumor cell proliferation under hypoxia (27). The expression of SLC2A5 can be regulated by tumor hypoxia, and overexpression of SLC2A5 promotes cell growth and metastasis in lung adenocarcinoma by enhancing fructose utilization (28). Despite the importance of stemness and hypoxia in cancers and the close association between cancer stemness and hypoxia, the prognostic capability of gene related to stemness-hypoxia in HCC has not been fully reported.

Herein, we calculated the mRNA expression-based stemness index (mRNAsi) and downloaded the hypoxia-related genes. Based on the above work, we aimed to construct a stemness-hypoxia gene signature to evaluate the prognosis of HCC patients. The workflow is provided in Figure S1. After the validation of our stemness-hypoxia gene signature, we further investigated the association between this model with clinical features, tumor immune microenvironment (TIME), and expression of immune checkpoints, which revealed a close relationship between the stemness-hypoxia gene signature and intra-tumoral immune landscape. Our study provides novel insights into the prognosis evaluation and clinical management of HCC. We present this article in accordance with the TRIPOD reporting checklist (available at https://tcr.amegroups.com/article/view/10.21037/tcr-24-2030/rc).

Methods

Data collection and processing

The study was conducted in accordance with the Declaration of Helsinki and its subsequent amendments. Clinical information and the RNA-seq data of 365 samples were obtained from TCGA-HCC database (29). Combining ID/Gene annotation file with Gencode database, the Ensembl_ID was converted into Symbol_ID. For those Symbol_ID corresponding with multiple Ensembl_IDs, we designated the mean as the expression value of the genes. Genes containing 0 in more than 80 % of the samples were excluded in this study.

Characterization of mRNAsi

mRNAsi of HCC samples was calculated by using one-class logistic regression (OCLR) machine-learning from gelnet package (version 1.2.1 in R.3.6.1) (12).

Acquisition of stemness-hypoxia genes

Hypoxia-related genes were downloaded from hallmark dataset in Molecular Signatures Database (MSigDB) (30), followed by the calculation of Pearson’s correlation coefficient (PCC) of the expression between mRNAsi and hypoxia-related genes using cor function in R. Genes with absolute PCC value >0.15 and with P<0.05 were selected as stemness-hypoxia genes.

Construction of stemness-hypoxia gene interaction network and enrichment analysis

The protein-protein interaction (PPI) network was constructed by searching the Search Tool for the Retrieval of Interacting Genes/Proteins (STRING) database (version 11.0) (31) with threshold of interaction score >0.7 and the degree ≥1. The PPI network was then visualized using Cytoscape (version 3.6.1) (32). Kyoto Encyclopedia of Genes and Genomes (KEGG) and Gene Ontology Biological Processes (GO-BP) annotation were performed using Database for Annotation, Visualization, and Integrated Discovery (DAVID) Bioinformatics Resources. False discovery rate (FDR) <0.05 was used as a cut-off value.

Screening of prognostic stemness-hypoxia genes

A total of 344 samples were obtained after exclusion of 21 samples from patients with overall survival (OS) <30 days. OS and OS status were aligned for each sample. The samples were divided into the training set (229 samples) and the validation set (115 samples) randomly by using survival (version 3.2-7) (33). Using the survival package and the clinical information, the genes related to stemness-hypoxia, which had significant association with prognosis, were chosen from the training set by using univariable Cox hazard analysis. The selected genes were further analyzed using multivariate Cox hazard analysis to select genes with independent prognostic value. Log-rank P<0.05 was used as a cut-off value.

Construction of the prognostic model

To eliminate the mutual influence of the independent prognostic genes screened above, the optimal gene combination for constructing the prognostic model was further screened using least absolute shrinkage and selection operator (LASSO) Cox regression analysis the glmnet package (version 2.0-18) in R3.6.1 (34). The prognostic model was established based on the regression coefficient (β) and the expression level of genes related to stemness-hypoxia obtained from multivariate Cox hazard analysis. The formula is:

Riskscore=∑βstemness-hypoxiagene×Expstemness-hypoxiagene[1]

where β_{stemness-hypoxia gene} represents regression coefficient and Exp_{stemness-hypoxia gene} represents the expression of stemness-hypoxia genes in TCGA dataset.

Validation of the prognostic model

The regression coefficient (β) obtained from the training set was used for validation in the validation set and whole dataset. The independent prognosis value in the training set, validation set, and whole dataset was evaluated using Kaplan-Meier survival analysis. Patients were divided into High- and low-risk groups based on the median value of the risk score, and the significance between the high-risk group and the low-risk group was examined.

Analysis of the prognostic characteristics of genes in the model

The genes in our prognostic model were further divided into high-expression (equal to or higher than the median) and low-expression (lower than the median) group. The gene expression-based Kaplan-Meier survival curve was drawn to evaluate the association of prognosis and gene expression.

Association between the risk groups and clinical features

The differences in clinical features [pathologic Tumor-Node-Metastasis (TNM) stage, tumor stage, age, gender, and stemness index] between the high- and low-risk groups were examined by using ggstatsplot (version 0.5.0) in R package. The percentage of each clinical feature was calculated. The P values were calculated by using the Chi-squared test.

Analysis of the independence of the prognostic model and construction of the nomogram

The independence of the prognostic model was determined by using univariate and multivariate Cox analyses with a cut-off P<0.05. Nomogram was built using rms (version 6.1-0) (https://hbiostat.org/R/rms/) in R (35), and C-index was calculated to assess the prediction ability of the nomogram.

Association between the high-risk/low-risk group and TIME

The CIBERSORT algorithm was used to analyze 22 types of tumor infiltrating immune cells [T cell CD8, T cell CD4 naïve, T cell CD4 memory resting, T cell CD4 memory activated, T cell follicular helper, T cell regulatory (Treg), T cell gamma delta, B cell naïve, B cell memory, plasma cells, resting natural killer (NK) cells, activated NK cells, resting mast cells, activated mast cells, eosinophils, neutrophils, monocytes, macrophages M0, macrophages M1, macrophages M2, resting dendritic cells, and activated dendritic cells] in the high-risk and low-risk groups (36).

Association between the high-risk/low-risk group and the expression of immune checkpoint genes

Immune checkpoint genes, such as programmed cell death 1 (PDCD1) and PD-1/PD-L1, which are involved in immune co-stimulatory and inhibitory pathways, should be regarded as potential targets for immune checkpoint blockade therapy (37). Hence, the expression levels of immune checkpoint genes, including PDCD1, CD274 (also known as PDL1), cytotoxic T-lymphocyte associated protein 4 (CTLA4), hepatitis A virus cellular receptor 2 (HAVCR2, also known as TIM3), lymphocyte activating gene 3 (LAG3), and V-Set domain containing T cell activation inhibitor 1 (VTCN1, also known as B7H4), were determined in the high-risk and low-risk groups.

KEGG pathway enrichment analysis in the high-risk and low-risk groups

The KEGG pathway enrichment was performed using the gene set variation analysis (version 1.36.2) (38) based on c2.cp.kegg.v7.2.symbols.gmt in MSigDB (version 7.2). The enrichment scores were further analyzed using the limma package in R, and P<0.05 were used as a cut-off value. Significant KEGG pathways were selected to build a heatmap using pheatmap (version 1.0.8).

Cell culture

Human HCC cell lines (HepG2 and Hep3B) and human normal liver cell line (LO-2) were obtained from Procell (Wuhan, China) and Shanghai MEIXUAN Biological science and technology LTD (Shanghai, China), respectively. All cells were maintained in Dulbecco’s modified eagle medium (DMEM) (Invitrogen, Carlsbad, CA, USA) supplemented with 10% fetal bovine serum (Gibco) and 1% penicillin (Gibco) at 37 ℃ and 5% CO₂.

Cell treatment

To investigate the effect of hypoxia on the expression levels of model genes, HepG2 and Hep3B cells were cultured in a hypoxia incubator: the oxygen concentration was adjusted to 1% to simulate the actual hypoxic environment inside tumors (1% O₂, 5% CO₂, and 94% N₂; 1% O₂ group). Moreover, HepG2 and Hep3B were treated with 150 µM cobalt chloride (CoCl₂) in 21% O₂ for 24 h to simulate the effects of hypoxia, which served as a positive control group for comparison with the actual hypoxic environment (CoCl₂ group). The normal conditions (21% O₂ and 5% CO₂; 21% O₂ group) were used as a negative control.

Quantitative reverse transcription polymerase chain reaction (qRT-PCR)

Total RNA was extracted using TRIzol reagent (TaKaRa, Otsu, Japan) and synthesized into cDNA using primeScript RT Master MIX (TaKaRa, Otsu, Japan) following the manufacturer’s instructions. qRT-PCR was performed using SYBR Green PCR Master Mix (Thermo, Waltham, MA, USA). Glyceraldehyde 3-phosphate dehydrogenase (GAPDH) was utilized as an endogenous reference. The primers sequences are listed in Table S1.

Statistical analysis

The experimental data were expressed as mean ± standard deviation. Statistical analyses were completed using GraphPad Prism 9.5.0 (GraphPad Software, San Diego, CA, USA). One-way analysis of variance (ANOVA) was utilized to analyze the differences between groups. The value of P<0.05 indicated statistically significant.

Results

Evaluation of mRNAsi and identification of genes related to stemness-hypoxia

mRNAsi was calculated by using the OCLR machine learning algorithm (12). Among the 200 hypoxia genes, 191 aligned to the expression dataset. A total of 109 genes were designated as genes related to stemness-hypoxia, with a cut-off of corValue >0.15 and P<0.05 (Table S2). The representative stemness-hypoxia genes are shown in Figure S2. Specially, SLC2A5 expression was positively correlated with mRNAsi, whereas PPARGC1A and PFKFB3 expression was negatively related to mRNAsi.

PPI network and KEGG pathway and GO enrichment analyses of stemness-hypoxia signature

After obtaining the stemness-hypoxia genes, we investigated the relationship among the genes. The PPI network of the selected 109 genes was built using the STRING database. A total of 194 edges with 75 gene nodes were identified, and EGFR was localized at the center of the PPI network (Figure S3A).

The biological functions of the genes were also evaluated by conducting KEGG pathway and GO enrichment analyses. The analyses revealed 36 significantly enriched BP and 16 KEGG pathways (Figure S3B,S3C). Both the BP and KEGG pathways exhibited metabolism- and hypoxia-related features, which are closely associated with responses to hypoxia. These functional analyses demonstrated that we successfully identified hypoxia-associated genes.

Construction of stemness-hypoxia prognostic model

We evaluated the prognostic ability of the 109 stemness-hypoxia genes in the training set. Using univariate Cox analysis, we found that 37 genes correlated significantly with HCC prognosis (Table S3). Further, in the multivariate Cox analysis, five genes of the 37 genes, namely cyclin G2 (CCNG2), peroxisome proliferator-activated receptor gamma coactivator 1-alpha (PPARGC1A), 6-phosphofructo-2-kinase/fructose-2,6-biphosphatase 3 (PFKFB3), solute carrier family 2 member 5 (SLC2A5), and RAR-related orphan receptor A (RORA), could independently predict HCC prognosis and used to construct the risk model (Table 1). On the basis of the five independent prognostic genes, three optimal genes (PFKFB3, SLC2A5, and PPARGC1A) were further screened using LASSO COX regression analysis, and their β-coefficient was determined (Figure S4 and Table 2). The risk model was constructed as follows: risk score = PFKFB3 × 0.277590888 + SLC2A5 × 0.274758201 + PPARGC1A × (−0.39065252). PFKFB3 and SLC2A5 were unfavorable factors with hazard ratio (HR) >1, whereas PPARGC1A was favorable factors for HCC patients (HR <1).

We constructed Core Risks model in the training set, validation set, and total set and divided the patients into high-risk (risk score > median) and low-risk (risk score ≤ median) groups. The Kaplan-Meier survival curves were plotted accordingly. In all the three sets, patients in the high-risk group showed a significantly lower survival probability than that seen for the low-risk group (Figure 1). These findings demonstrate that our prognostic model can effectively distinguish high-risk patients from low-risk patients.

Figure 1 KM survival analysis. (A) KM survival analysis of risk model in the training set. (B) KM survival analysis of risk model in the validation set. (C) KM survival analysis of risk model in the total set. KM, Kaplan-Meier.

Survival outcome according to the expression level of the genes in the prognostic model

We examined the association between the expression of three genes in the prognostic model with OS in the HCC transcriptomic dataset. The expression of PFKFB3 had a significant negative correlation with OS (P=0.003) (Figure S5A). SLC2A5 showed a similar but milder trend as PFKFB3 (Figure S5B). Elevated expression levels of PPARGC1A (P=5e−4) were favorable to the OS of HCC patients (Figure S5C).

Comparison of the clinical features of the high-risk and low-risk groups

We compared the clinical features in the high-risk group with those in the low-risk group. A significant difference in tumor stages was observed (P=1.38e−04). In the high-risk group, 62% of the patients (stage I and II) were in the early stage, but 78% of the patients were in the early stage in the low-risk group (Figure 2A). Additionally, significant difference was also exhibited in pathologic T stages (P=1.22e−04, Figure 2B). However, there were no differences in the other clinical features including age, gender, stemness (mRNAsi), pathologic M stages, and pathologic N stages (Figure S6).

Figure 2 Clinical features analysis. (A) The ratio of patients with different tumor stages in the high-risk and low-risk groups. (B) The ratio of patients with different pathologic T stages in the high-risk and low-risk groups. CI, confidence interval; T, tumor.

Independent prognostic clinical features

To investigate the clinical feature that independently predicted the prognostic outcome of HCC, we employed the univariate and multivariate Cox analyses to assess the potential of the clinical features, including pathologic TNM stage, risk group, mRNAsi, gender, age, distant metastatic status, and regional lymph node invasion. We observed a significant association of pathologic T stage, risk score, and mRNAsi with OS (Figure 3A,3B). Subsequently, we constructed a nomogram with the three independent prognostic clinical features to visualize the prognostic outcome (Figure 3C). Furthermore, we evaluated the independency of our prognostic model, and the results demonstrated that the prognostic model could independently predict the prognostic outcome of HCC from TCGA dataset (Figure 3D).

Figure 3 Nomogram of independent prognostic factors. (A) Univariate Cox regression analysis to screen independent prognostic factors. (B) Multivariate Cox regression analysis to screen independent prognostic factors. (C) The nomogram was constructed based on the independent prognostic factors (risk score, pathologic T stage, and mRNAsi). (D) The correlation of nomogram-predicted survival probability with the actual OS proportion in HCC. CI, confidence interval; HCC, hepatocellular carcinoma; M, metastasis; mRNAsi, mRNA expression-based stemness index; N, lymph node; OS, overall survival; T, tumor.

Evaluation of tumoral immune cell infiltration in high-risk and low-risk groups

Immune responses to anti-tumor therapy are closely associated with clinical outcome. Thus, we assessed immune cell infiltration into tumor tissues by using CIBERSORT (36). A total of 22 infiltrating immune subpopulations were evaluated (Figure 4). Of the subpopulations, B cells memory, T cells CD4 memory resting, and macrophages M0 showed significant difference between the high-risk and low-risk groups. These differences in tumor infiltrating immune cells, at least partially, explain the ability of our prognostic model to predict the clinical outcome of HCC patients in the TCGA dataset.

Figure 4 Analysis of tumor infiltration of immune cells. The infiltration of 22 immune subpopulations was evaluated in TCGA-HCC dataset. The red figure indicates high-risk group, whereas the blue figure represents low-risk group. HCC, hepatocellular carcinoma; NK, natural killer; TCGA, The Cancer Genome Atlas.

Association between the expression of immune checkpoint genes and risk groups

The correlation of immune cell infiltration with risk scores prompted us to check the expression of immune checkpoint genes, which are strongly associated with the responses to and the outcome of immunotherapies. Except LAG3, the immune checkpoint genes (PDCD1, HAVCR2, CTLA4, VTCN1, and CD274) showed significantly higher expression in the high-risk group (P<0.05), implying a more suppressive TIME (Figure 5).

Figure 5 Comparison of the expression of immune checkpoints in the high-risk and low-risk groups. The expression level of immune checkpoints PDCD1 (A), LAG3 (B), HAVCR2 (C), CTLA4 (D), VTCN1 (E) and CD274 (F) in the high-risk and low-risk groups. *, P<0.05; **, P<0.01; ***, P<0.001. CTLA4, cytotoxic T-lymphocyte associated protein 4; HAVCR2, hepatitis A virus cellular receptor 2; LAG3, lymphocyte activating gene 3; ns, non-significance; PDCD1, programmed cell death 1; VTCN1, V-Set domain containing T cell activation inhibitor 1.

KEGG pathway enrichment analysis of high-risk and low-risk groups

We investigated the pathways that were correlated with the risk score. The KEGG pathway enrichment analysis identified 37 pathways, dominated by cellular metabolism and signal transduction pathways that were significantly associated with the risk score. The top 5 up- and down-regulated pathways were visualized in Figure 6. These data suggested that the pathways might be involved in the regulation of tumor progression and clinical outcome.

Figure 6 KEGG pathways enrichment analysis. The differences in the enriched pathways between the high-risk and low-risk groups were analyzed. KEGG, Kyoto Encyclopedia of Genes and Genomes; M, metastasis; N, lymph node; T, tumor.

Validation of the expression of independent prognostic genes in HCC

LO-2 and HepG2 cells were used to further explore the expression of the independent prognostic genes (PPARGC1A, PFKFB3, and SLC2A5). The results showed that the expression of SLC2A5 was significantly highly expressed in HepG2 cells relative to LO-2 cells (Figure 7A). However, PPARGC1A expression did not exhibited significant difference between HepG2 and LO-2 cells (Figure 7B). However, there was no result of PFKFB3 because its expression was too low.

Figure 7 The relative expression of model genes in the human HCC cell line (HepG2) and human normal liver cell line (LO-2). (A) The relative expression of SLC2A5 in HepG2 and LO-2 cell lines. (B) The relative expression of PPARGC1A in HepG2 and LO-2 cell lines. *, P<0.05. HCC, hepatocellular carcinoma; PPARGC1A, peroxisome proliferator-activated receptor gamma coactivator 1-alpha; SLC2A5, solute carrier family 2 member 5.

Analysis of the expression of independent prognostic genes under hypoxia condition

We further investigated the effect of hypoxia on the expression levels of the independent prognostic genes (PPARGC1A, PFKFB3, and SLC2A5) in Hep3B and HepG2 cells. As results, in both Hep3B and HepG2 cells, SLC2A5 expression was significantly up-regulated in 1% O₂ and CoCl₂ groups compared to that in 21% O₂ group, whereas PPARGC1A and PFKFB3 expression was remarkably down-regulated in 1% O₂ and CoCl₂ groups (Figure 8). These data indicated that these genes may be key regulators in HCC under hypoxia condition.

Figure 8 The relative expression of model genes (SLC2A5, PPARGC1A, and PFKFB3) in the human HCC cell lines under hypoxic conditions. (A) The relative expression of model genes in Hep3B cells under hypoxic conditions. (B) The relative expression of model genes in HepG2 cells under hypoxic conditions. *, P<0.05; **, P<0.01; ***, P<0.001. HCC, hepatocellular carcinoma; ns, non-significance; PFKFB3, 6-phosphofructo-2-kinase/fructose-2,6-biphosphatase 3; PPARGC1A, peroxisome proliferator-activated receptor gamma coactivator 1-Alpha; SLC2A5, solute carrier family 2 member 5.

Discussion

HCC is one of the deadliest cancers, and the prognosis is always unfavorable (39). The evaluation and prediction of HCC outcomes are critical for the clinical management of the disease. Prognostic biomarkers show promising application value, but extensive studies are necessary for clinical translation. Our present study established a novel stemness-hypoxia prognostic signature for HCC that did not only distinguish high-risk patients from low-risk patients but also showed close associations with TIME, providing insights into HCC prognosis and anti-tumor therapy, such as immunotherapy.

Both cancer stemness and hypoxia condition are critical for HCC progression (10,11,40); therefore, few studies have employed the stemness and hypoxia gene signature to evaluate the prognosis of HCC (21-25,41,42). We established a novel stemness-hypoxia gene signature (PPARGC1A, PFKFB3, and SLC2A5) for predicting HCC prognosis using data from TCGA database. PPARGC1A is the master regulator of mitochondrial biogenesis (43), which plays essential roles in HCC development (44,45). PPAGC1A is lowly expressed in HCC, and PPARGC1A/BAMBI/ACSL5 axis is hypoxia-responsive in the progression of HCC (26). PFKFB3 regulates glycolysis and promotes HCC progression (46). PFKFB3 can be hypoxia inducible and contribute to hypoxia-induced glycolysis in cancer cells (47). SLC2A5, also known as GLUT5, is responsible for fructose transportation, and it is also associated with HCC (48). Under the condition of cancer, SLC2A5 (GLUT5) expression can be induced by hypoxia, and there is a positive correlation between HIF-1α and GLUT5 in breast cancer (49). In this study, a prognostic signature constructed by these genes could distinguish low-risk patients from high-risk patients and had a good ability for predicting patients’ OS. Therefore, we conclude that this three-gene signature is valuable in predicting HCC prognosis and these genes may be novel prognostic biomarkers for HCC. Additionally, we found that SLC2A5 expression was positively correlated with mRNAsi, whereas PPARGC1A and PFKFB3 expression was negatively related to mRNAsi. qRT-PCR showed elevated SLC2A5 expression and reduced PPARGC1A and PFKFB3 expression in both Hep3B and HepG2 cells under hypoxia condition. These data indicate that these genes are related to tumor stemness and are hypoxia-responsive in the progression of HCC. However, future in vitro and in vivo experiments on model genes are required to elucidate their roles in HCC.

Previous studies have shown that tumor stemness and hypoxia are closely associated with intra-tumoral immune landscape (50-53). In this study, we examined the relationship between stemness-hypoxia prognostic signature

with immune characteristics. We observed increased infiltration of tumor-promoting immune subpopulations (Tregs, macrophages M0, and neutrophils) and decreased tumoricidal subpopulations (CD4⁺ memory T cells, and NK activated cells) in the high-risk group. Although the function of each immune subpopulation in the tumor microenvironment is context-dependent, the landscape of the immune subpopulations implies a more suppressive TIME in the high-risk group. In addition, significantly elevated expression of immune checkpoint genes (PDCD1, LAG3, HAVCR2, CTLA4, VTCN1, and CD274) was observed in the high-risk group, suggesting profound immune evasion. Thus, this immune landscape is associated with the risk score and can, at least partially, explain the poor prognosis of patients in the high-risk group. Furthermore, the immune checkpoint genes with elevated levels in the high-risk group indicate that they can be targeted by immune checkpoint blockade therapies.

To explore the biological mechanism underlying the stemness-hypoxia prognostic signature, we analyzed KEGG pathways associated with this signature. Changes in cell metabolism in adaptation to hypoxia occur to promote cell survival and proliferation (54). An enrichment of metabolism-related pathways, such as glycan biosynthesis and purine metabolism was observed. These results suggest that metabolic adaptability is a key feature of tumor development. In addition, the p53 and MAPK signaling pathways are important for tumor progression, and these processes are associated with responses to hypoxia and metabolism adaptation. The p53 signaling pathway is a molecular switch for stem cells to enter the cell cycle and is associated with stem cell self-renewal (55). DNA damage triggers a series of harmful consequences, while p53 signaling is activated to induce cell cycle arrest, cellular senescence or apoptosis, which in turn reduces the damage (56). MAPK/ERK signaling nodes can act as both tumor suppressor and pro-oncogene signal. The main reasons for the different results are the intensity and context of MAPK/ERK signaling (57). In this study, the enrichment of the p53 and MAPK pathways in the low-risk group suggests the association of the pathways with cancer stemness and hypoxia, which may contribute to greater prognosis of HCC patients in the low-risk group.

This study had some limitations. First, only TCGA-HCC dataset was used in this study and experimental results may not be broadly representative due to heterogeneity between samples. In future studies, we will further expand the sample size to reduce the bias of the study. Second, the expression of three model genes was not explored using clinical samples, and the clinical applicability of the stemness-hypoxia prognostic model for HCC was not validated using clinical cohorts. Third, we showed a strong association between our model and TIME, but their complex interaction needs to be explored.

Conclusions

Our work established a stemness-hypoxia gene signature showed favorable prognosis performance. Moreover, we identified suppressive immune signature and elevated expression levels of immune checkpoint genes in the high-risk group, which may be significant for immunotherapy. The findings of this study provide some novel insights into the application of stemness-hypoxia signature in the clinical management of HCC.

Acknowledgments

None.

Footnote

Reporting Checklist: The authors have completed the TRIPOD reporting checklist. Available at https://tcr.amegroups.com/article/view/10.21037/tcr-24-2030/rc

Data Sharing Statement: Available at https://tcr.amegroups.com/article/view/10.21037/tcr-24-2030/dss

Peer Review File: Available at https://tcr.amegroups.com/article/view/10.21037/tcr-24-2030/prf

Funding: None.

Conflicts of Interest: All authors have completed the ICMJE uniform disclosure form (available at https://tcr.amegroups.com/article/view/10.21037/tcr-24-2030/coif). The authors have no conflicts of interest to declare.

Ethical Statement: The authors are accountable for all aspects of the work in ensuring that questions related to the accuracy or integrity of any part of the work are appropriately investigated and resolved. The study was conducted in accordance with the Declaration of Helsinki and its subsequent amendments.

Open Access Statement: This is an Open Access article distributed in accordance with the Creative Commons Attribution-NonCommercial-NoDerivs 4.0 International License (CC BY-NC-ND 4.0), which permits the non-commercial replication and distribution of the article with the strict proviso that no changes or edits are made and the original work is properly cited (including links to both the formal publication through the relevant DOI and the license). See: https://creativecommons.org/licenses/by-nc-nd/4.0/.

References

1. Huang M, Chen H, Wang H, et al. Worldwide burden of liver cancer due to metabolic dysfunction-associated steatohepatitis from 1990 to 2019: insights from the Global Burden of Disease study. Front Oncol 2024;14:1424155. [Crossref] [PubMed]
2. Cheu JW, Wong CC. Mechanistic Rationales Guiding Combination Hepatocellular Carcinoma Therapies Involving Immune Checkpoint Inhibitors. Hepatology 2021;74:2264-76. [Crossref] [PubMed]
3. Llovet JM, Montal R, Villanueva A. Randomized trials and endpoints in advanced HCC: Role of PFS as a surrogate of survival. J Hepatol 2019;70:1262-77. [Crossref] [PubMed]
4. Oh JH, Sinn DH. Multidisciplinary approach for hepatocellular carcinoma patients: current evidence and future perspectives. J Liver Cancer 2024;24:47-56. [Crossref] [PubMed]
5. Akabane M, Imaoka Y, Kawashima J, et al. Advancing precision medicine in hepatocellular carcinoma: current challenges and future directions in liquid biopsy, immune microenvironment, single nucleotide polymorphisms, and conversion therapy. Hepat Oncol 2025;12:2493457. [Crossref] [PubMed]
6. Peitzsch C, Tyutyunnykova A, Pantel K, et al. Cancer stem cells: The root of tumor recurrence and metastases. Semin Cancer Biol 2017;44:10-24. [Crossref] [PubMed]
7. Liu Q, Guo Z, Li G, et al. Cancer stem cells and their niche in cancer progression and therapy. Cancer Cell Int 2023;23:305. [Crossref] [PubMed]
8. Kim J, Orkin SH. Embryonic stem cell-specific signatures in cancer: insights into genomic regulatory networks and implications for medicine. Genome Med 2011;3:75. [Crossref] [PubMed]
9. Walcher L, Kistenmacher AK, Suo H, et al. Cancer Stem Cells-Origins and Biomarkers: Perspectives for Targeted Personalized Therapies. Front Immunol 2020;11:1280. [Crossref] [PubMed]
10. Yamashita T, Kaneko S. Liver cancer stem cells: Recent progress in basic and clinical research. Regen Ther 2021;17:34-7. [Crossref] [PubMed]
11. Tsui YM, Chan LK, Ng IO. Cancer stemness in hepatocellular carcinoma: mechanisms and translational potential. Br J Cancer 2020;122:1428-40. [Crossref] [PubMed]
12. Malta TM, Sokolov A, Gentles AJ, et al. Machine Learning Identifies Stemness Features Associated with Oncogenic Dedifferentiation. Cell 2018;173:338-354.e15. [Crossref] [PubMed]
13. Petrova V, Annicchiarico-Petruzzelli M, Melino G, et al. The hypoxic tumour microenvironment. Oncogenesis 2018;7:10. [Crossref] [PubMed]
14. Brahimi-Horn MC, Chiche J, Pouysségur J. Hypoxia and cancer. J Mol Med (Berl) 2007;85:1301-7. [Crossref] [PubMed]
15. Kang N, Choi SY, Kim BN, et al. Hypoxia-induced cancer stemness acquisition is associated with CXCR4 activation by its aberrant promoter demethylation. BMC Cancer 2019;19:148. [Crossref] [PubMed]
16. Mimeault M, Batra SK. Hypoxia-inducing factors as master regulators of stemness properties and altered metabolism of cancer- and metastasis-initiating cells. J Cell Mol Med 2013;17:30-54. [Crossref] [PubMed]
17. Yun Z, Lin Q. Hypoxia and regulation of cancer cell stemness. Adv Exp Med Biol 2014;772:41-53. [Crossref] [PubMed]
18. Zhou Y, Fan W, Xiao Y. The effect of hypoxia on the stemness and differentiation capacity of PDLC and DPC. Biomed Res Int 2014;2014:890675. [Crossref] [PubMed]
19. Cui CP, Wong CC, Kai AK, et al. SENP1 promotes hypoxia-induced cancer stemness by HIF-1α deSUMOylation and SENP1/HIF-1α positive feedback loop. Gut 2017;66:2149-59. [Crossref] [PubMed]
20. Ling S, Shan Q, Zhan Q, et al. USP22 promotes hypoxia-induced hepatocellular carcinoma stemness by a HIF1α/USP22 positive feedback loop upon TP53 inactivation. Gut 2020;69:1322-34. [Crossref] [PubMed]
21. Jubb AM, Buffa FM, Harris AL. Assessment of tumour hypoxia for prediction of response to therapy and cancer prognosis. J Cell Mol Med 2010;14:18-29. [Crossref] [PubMed]
22. Walsh JC, Lebedev A, Aten E, et al. The clinical importance of assessing tumor hypoxia: relationship of tumor hypoxia to prognosis and therapeutic opportunities. Antioxid Redox Signal 2014;21:1516-54. [Crossref] [PubMed]
23. Gomez CR. Editorial: Tumor Hypoxia: Impact in Tumorigenesis, Diagnosis, Prognosis, and Therapeutics. Front Oncol 2016;6:229. [Crossref] [PubMed]
24. Zhang B, Tang B, Gao J, et al. A hypoxia-related signature for clinically predicting diagnosis, prognosis and immune microenvironment of hepatocellular carcinoma patients. J Transl Med 2020;18:342. [Crossref] [PubMed]
25. Zeng F, Zhang Y, Han X, et al. Employing hypoxia characterization to predict tumour immune microenvironment, treatment sensitivity and prognosis in hepatocellular carcinoma. Comput Struct Biotechnol J 2021;19:2775-89. [Crossref] [PubMed]
26. Zhang Q, Xiong L, Wei T, et al. Hypoxia-responsive PPARGC1A/BAMBI/ACSL5 axis promotes progression and resistance to lenvatinib in hepatocellular carcinoma. Oncogene 2023;42:1509-23. [Crossref] [PubMed]
27. Lei Y, Chen T, Li Y, et al. O-GlcNAcylation of PFKFB3 is required for tumor cell proliferation under hypoxia. Oncogenesis 2020;9:21. [Crossref] [PubMed]
28. Weng Y, Fan X, Bai Y, et al. SLC2A5 promotes lung adenocarcinoma cell growth and metastasis by enhancing fructose utilization. Cell Death Discov 2018;4:38. [Crossref] [PubMed]
29. Goldman MJ, Craft B, Hastie M, et al. Visualizing and interpreting cancer genomics data via the Xena platform. Nat Biotechnol 2020;38:675-8. [Crossref] [PubMed]
30. Liu Y, Wu J, Huang W, et al. Development and validation of a hypoxia-immune-based microenvironment gene signature for risk stratification in gastric cancer. J Transl Med 2020;18:201. [Crossref] [PubMed]
31. Szklarczyk D, Gable AL, Lyon D, et al. STRING v11: protein-protein association networks with increased coverage, supporting functional discovery in genome-wide experimental datasets. Nucleic Acids Res 2019;47:D607-13. [Crossref] [PubMed]
32. Otasek D, Morris JH, Bouças J, et al. Cytoscape Automation: empowering workflow-based network analysis. Genome Biol 2019;20:185. [Crossref] [PubMed]
33. Therneau TM, Grambsch PM. Modeling Survival Data: Extending the Cox Model. Springer; 2000. ISBN 0-387-98784-3.
34. Friedman J, Hastie T, Tibshirani R. Glmnet: Lasso and elastic-net regularized generalized linear models. 2009. Available online: https://cran.r-project.org/web/packages/glmnet/index.html
35. Zhang Z, Geskus RB, Kattan MW, et al. Nomogram for survival analysis in the presence of competing risks. Ann Transl Med 2017;5:403. [Crossref] [PubMed]
36. Newman AM, Liu CL, Green MR, et al. Robust enumeration of cell subsets from tissue expression profiles. Nat Methods 2015;12:453-7. [Crossref] [PubMed]
37. Hu FF, Liu CJ, Liu LL, et al. Expression profile of immune checkpoint genes and their roles in predicting immunotherapy response. Brief Bioinform 2021;22:bbaa176. [Crossref] [PubMed]
38. Hänzelmann S, Castelo R, Guinney J. GSVA: gene set variation analysis for microarray and RNA-seq data. BMC Bioinformatics 2013;14:7. [Crossref] [PubMed]
39. Ren Z, Ma X, Duan Z, et al. Diagnosis, Therapy, and Prognosis for Hepatocellular Carcinoma. Anal Cell Pathol (Amst) 2020;2020:8157406. [Crossref] [PubMed]
40. Wu XZ, Xie GR, Chen D. Hypoxia and hepatocellular carcinoma: The therapeutic target for hepatocellular carcinoma. J Gastroenterol Hepatol 2007;22:1178-82. [Crossref] [PubMed]
41. Bai KH, He SY, Shu LL, et al. Identification of cancer stem cell characteristics in liver hepatocellular carcinoma by WGCNA analysis of transcriptome stemness index. Cancer Med 2020;9:4290-8. [Crossref] [PubMed]
42. Kong W, Gao M, Jin Y, et al. Prognostic model of patients with liver cancer based on tumor stem cell content and immune process. Aging (Albany NY) 2020;12:16555-78. [Crossref] [PubMed]
43. Dorn GW 2nd, Vega RB, Kelly DP. Mitochondrial biogenesis and dynamics in the developing and diseased heart. Genes Dev 2015;29:1981-91. [Crossref] [PubMed]
44. Yun SH, Han SH, Park JI. Peroxisome Proliferator-Activated Receptor γ and PGC-1α in Cancer: Dual Actions as Tumor Promoter and Suppressor. PPAR Res 2018;2018:6727421. [Crossref] [PubMed]
45. Liu R, Zhang H, Zhang Y, et al. Peroxisome proliferator-activated receptor gamma coactivator-1 alpha acts as a tumor suppressor in hepatocellular carcinoma. Tumour Biol 2017;39:1010428317695031. [Crossref] [PubMed]
46. Shi WK, Zhu XD, Wang CH, et al. PFKFB3 blockade inhibits hepatocellular carcinoma growth by impairing DNA repair through AKT. Cell Death Dis 2018;9:428. [Crossref] [PubMed]
47. Yi M, Ban Y, Tan Y, et al. 6-Phosphofructo-2-kinase/fructose-2,6-biphosphatase 3 and 4: A pair of valves for fine-tuning of glucose metabolism in human cancer. Mol Metab 2019;20:1-13. [Crossref] [PubMed]
48. Gillet JP, Andersen JB, Madigan JP, et al. A Gene Expression Signature Associated with Overall Survival in Patients with Hepatocellular Carcinoma Suggests a New Treatment Strategy. Mol Pharmacol 2016;89:263-72. [Crossref] [PubMed]
49. Hamann I, Krys D, Glubrecht D, et al. Expression and function of hexose transporters GLUT1, GLUT2, and GLUT5 in breast cancer-effects of hypoxia. FASEB J 2018;32:5104-18. [Crossref] [PubMed]
50. Krzywinska E, Stockmann C. Hypoxia, Metabolism and Immune Cell Function. Biomedicines 2018;6:56. [Crossref] [PubMed]
51. Noman MZ, Hasmim M, Messai Y, et al. Hypoxia: a key player in antitumor immune response. A Review in the Theme: Cellular Responses to Hypoxia. Am J Physiol Cell Physiol 2015;309:C569-79. [Crossref] [PubMed]
52. Chen P, Hsu WH, Han J, et al. Cancer Stemness Meets Immunity: From Mechanism to Therapy. Cell Rep 2021;34:108597. [Crossref] [PubMed]
53. Miranda A, Hamilton PT, Zhang AW, et al. Cancer stemness, intratumoral heterogeneity, and immune response across cancers. Proc Natl Acad Sci U S A 2019;116:9020-9. [Crossref] [PubMed]
54. Eales KL, Hollinshead KE, Tennant DA. Hypoxia and metabolic adaptation of cancer cells. Oncogenesis 2016;5:e190. [Crossref] [PubMed]
55. Ginestier C, Wicinski J, Cervera N, et al. Retinoid signaling regulates breast cancer stem cell differentiation. Cell Cycle 2009;8:3297-302. [Crossref] [PubMed]
56. Shi T, Dansen TB. Reactive Oxygen Species Induced p53 Activation: DNA Damage, Redox Signaling, or Both? Antioxid Redox Signal 2020;33:839-59. [Crossref] [PubMed]
57. Burotto M, Chiou VL, Lee JM, et al. The MAPK pathway across different malignancies: a new perspective. Cancer 2014;120:3446-56. [Crossref] [PubMed]