A hypoxia-linked gene signature for prognosis prediction and evaluating the immune microenvironment in patients with hepatocellular carcinoma.

BACKGROUND: Previous research indicates that hypoxia critically affects the initiation and progression of hepatocellular carcinoma (HCC). Nevertheless, the molecular mechanisms responsible for HCC development are poorly understood. Herein, we purposed to build a prognostic model using hypoxia-linked genes to predict patient prognosis and investigate the relationship of hypoxia with immune status in the tumor microenvironment (TME).
METHODS: The training cohort included transcriptome along with clinical data abstracted from The Cancer Genome Atlas (TCGA). The validation cohort was abstracted from Gene Expression Omnibus (GEO). Univariate along with multivariate Cox regression were adopted to create the prediction model. We divided all patients into low- and high-risk groups using median risk scores. The estimation power of the prediction model was determined with bioinformatic tools.
RESULTS: Six hypoxia-linked genes, HMOX1, TKTL1, TPI1, ENO2, LDHA, and SLC2A1, were employed to create an estimation model. Kaplan-Meier, ROC curve, and risk plot analyses demonstrated that the estimation potential of the risk model was satisfactory. Univariate along with multivariate regression data illustrated that the risk model could independently predict the overall survival (OS). A nomogram integrating the risk signature and clinicopathological characteristics showed a good potential to estimate HCC prognosis. Gene set enrichment analysis (GSEA) revealed that genes associated with cell proliferation and metabolism cascades were abundant in high-risk group. Furthermore, the signature showed a strong ability to distinguish the two groups in terms of immune status.
CONCLUSIONS: A prediction model for predicting HCC prognosis using six hypoxia-linked genes was designed in this study, facilitating the diagnosis and treatment of HCC. 2021 Translational Cancer Research. All rights reserved.

Introduction

Liver cancer is a global leading cause of mortality ranked as the 4th cause of cancer-associated deaths (1,2). Hepatocellular carcinoma (HCC) is the most frequent primary liver malignancy in adults (3). As curative therapeutic HCC options consisting of surgical resection, ablation treatments along with liver transplantation are effective for early-stage HCC, the timely diagnosis coupled with early intervention have remarkable impacts on the prognosis of individuals with HCC (4). Because of the lack of symptoms in the early HCC stages coupled with the rapid progress of the disease, most of the individuals with HCC are diagnosed at an advanced stage and are unable to undergo surgical resection (5-7). To enhance the survival of individuals with HCC, it is necessary to find reliable tumor markers for early diagnosis and molecular targets for effective treatment.

Hypoxia is a common phenomenon in the majority of human tumors, and its occurrence results from an imbalance between the supply and consumption of oxygen (8,9). Due to rapid proliferation, HCC cells quickly exhaust nutrients and the oxygen supply, causing hypoxia in the tumor microenvironment (TME) (10,11). Furthermore, hypoxia has been reported to be a critical microenvironmental factor that enhances cancer proliferation, invasion, metastasis, and chemotherapy resistance in HCC (12-14). However, the detailed mechanism underlying hypoxia-induced proliferation and metastasis in HCC has not yet been fully elucidated.

To understand the complex mechanisms, attention has been increasingly paid to the crosstalk between HCC cells and tumor-infiltrating immune cells in the TME (15-17). Immune repressor cells consisting of regulatory T cells (Tregs) (18), tumor-associated macrophages (TAMs) (19), as well as myeloid-derived suppressive cells (MDSCs) (20), in the HCC tumor immune microenvironment (TIME) can facilitate tumor immune suppression or immune escape via establishment of an immune repressive TME. Moreover, immune checkpoint inhibitors (ICIs), including inhibitors of cytotoxic T lymphocyte-associated-4 (CTLA-4), programmed cell death protein 1 (PD-1), and programmed death-ligand 1 (PD-L1) (21), can help HCC cells escape cytotoxic T cell-mediated death (22). Understanding the HCC immune microenvironment will improve the development of immunotherapies. Interestingly, increasing evidence suggests that hypoxia affects components of the HCC immune microenvironment (23-25). Therefore, investigating the gene expression changes associated with hypoxia pathways might contribute to survival benefits for patients with HCC.

Herein, we first abstracted mRNA expression data along with clinical data of individuals with HCC from public data resources. Then, we established a prediction model based on hypoxia-linked genes using a The Cancer Genome Atlas (TCGA) cohort and verified its estimation ability in a Gene Expression Omnibus (GEO) data set. Furthermore, we explored the relationship of the hypoxia-linked signature with immune status in HCC. We hope that this study will make a positive contribution to the development of treatment approaches for individuals with HCC. We present the following article in accordance with the TRIPOD reporting checklist (available at https://dx.doi.org/10.21037/tcr-21-741).

Methods

Data collection

The transcriptome along with the matching clinicopathological information of individuals with HCC were abstracted from TCGA (https://portal.gdc.cancer.gov/) (including 374 HCC and 50 normal tissue samples) as a training cohort. Similarly, data from 81 HCC patients from a GEO (https://www.ncbi.nlm.nih.gov/geo/) dataset were downloaded as a verification cohort. Then, a set of hypoxia-linked genes (n=200) was abstracted from the Molecular Signature Database (MSigDB, V.6.0, USA) (http://www.gsea-msigdb.org/gsea/index.jsp). The study was conducted in accordance with the Declaration of Helsinki (as revised in 2013).

Construction of a protein-protein interaction (PPI) network

We constructed a PPI network through the Search Tool for Recurring Instances of Neighboring Genes (STRING) data resource (http://string-db.org) to analyze the interactions of the hypoxia-linked genes. Then, we screened core network nodes based on the number of interconnections, and those with the most interconnections were identified as the key genes for downstream analysis.

Construction and verification of a prediction model with hypoxia-linked genes

We first conducted univariate Cox regression to determine the genes with potential for predicting overall survival (OS) (P<0.05). Further, multivariate regression model was implemented to uncover prognostic genes to include in the prediction model. Computation of the risk score was determined using the formula:

where β designates the regression coefficient, and Exp designates the expression levels of hypoxia-linked genes. As per the median risk score, we divided the individuals with HCC in the TCGA cohort into high- and low-risk groups. Furthermore, 81 HCC samples extracted from the GSE54236 dataset were employed to confirm the estimation potential of the prediction model.

Assessment of risk model’ independence from other clinicopathological characteristics

Univariate along with multivariate regression model were implemented to explore if the risk signature was an independent OS predictive factor compared to conventional clinicopathological characteristics.

Establishment and verification of a prognostic nomogram

The factors that were remarkable in the multivariate assessment were employed to generate a nomogram to assess the one-, three-, and five-year survival probabilities of individuals with HCC. In addition, we conducted time-dependent ROC analyses to explore the accuracy of our predictive model. Furthermore, we generated calibration curves to estimate the predictive potential of the nomogram.

Gene set enrichment analysis (GSEA)

GSEA was adopted to explore remarkable differences in the expression of gene sets between the low- and high-risk groups based on the MSigDB collection (h.all.v7.2.symbols.gmt, USA) using GSEA software (http://www.gsea-msigdb.org/gsea/index.jsp). A thousand-fold gene set permutations were set for each analysis.

Single-sample gene set enrichment analysis (ssGSEA)

On the basis of the gene expression profiles, the ssGSEA approach was employed to quantify the proportions of 16 kinds of invading immune cells and the activity of 13 immune-linked cascades or roles (26).

Statistics analysis

Data analyses were implemented in the R software (V.4.0.2, MathSoft, USA). The six-gene estimation model was created by univariate along with multivariate regression models. Kaplan-Meier approach was adopted to compare the prognosis between the high- and low-risk groups. Besides, the log-rank test was adopted to determine the differences. ROC curves coupled with area under the ROC curve (AUC) values were adopted to explore the estimation ability of the hypoxia-linked risk signature. Univariate along with multivariate Cox regression were conducted to uncover independent OS predictors. A nomogram predicting the prognosis of individuals with HCC was generated based on the statistically significant predictors. ROC curves were created to estimate the accuracy of the model. We plotted a calibration curve to determine the estimation capacity of the nomogram. P<0.05 signified statistical significance, with all the P values being two-tailed.

Results

Construction of the prediction model with the TCGA cohort

The flow chart of this study is given in . A gene set associated with the hypoxia signaling cascade was abstracted from the GSEA web resource (n=200). To identify hub genes among the hypoxia-linked genes, we first constructed a PPI network using the STRING database. Then, we identified the top 50 genes with the highest cross talk degrees as key genes ().

Figure 1

Flow chart illustrating data collection along with analysis. HCC, hepatocellular carcinoma; TCGA, The Cancer Genome Atlas.

Figure 2

Construction of an estimation model on the basis of hypoxia-linked genes. (A) The top 50 genes according to their interaction degrees; (B) construction of the prediction model following univariate coupled with multivariate regression. Relationship of six hypoxia-linked genes in the TCGA (C) and GEO (D) cohorts. TCGA, The Cancer Genome Atlas; GEO, Gene Expression Omnibus.

Furthermore, univariate along with multivariate Cox regression were adopted to create a prediction model. According to the univariate analysis, 18 genes were remarkably related with the OS of individuals with HCC. According to the results of the multivariate analysis, six genes were chosen to build the prediction model, as illustrated in . The risk score of each sample was computed using the following formula:

These six genes were remarkably related with one another in the TCGA () and GEO () cohorts.

Evaluation of the prediction model with the TCGA cohort

To explore the estimation potential of the model, we first divided samples in the TCGA data set into high- (n=185), as well as low-risk (n=185) groups, as per the median risk score (). Following a series of bioinformatics analyses, a risk plot was generated, and it illustrated that high-risk score patients harboured a shorter survival in contrast with those with a low-risk score (). Furthermore, Kaplan-Meier analysis illustrated that the high-risk group harboured dismal OS in contrast with the low-risk group (). ROC curve assessment was adopted to determine the AUC values, which were 0.748 at one year, 0.714 at two years, and 0.719 at three years (). Additionally, the heatmap analysis illustrated that the levels of expressions of the six genes were remarkably different between different groups, implying that the six-gene predictive model can estimate the prognosis of individuals with HCC ().

Figure 3

Predictive ability of the hypoxia risk signature in HCC. The risk scores distribution and median value in the TCGA (A) data set and the GEO (F) cohort. Distribution of the OS status, survival time along with risk scores in TCGA (B) data set and the GEO (G) data set. Kaplan-Meier curves between the high- and low-risk groups for the OS of patients in the TCGA (C) data set and the GEO (H) data set. Time-dependent ROC curves verifying the predictive ability of the risk signature in the TCGA (D) data set and the GEO (I) cohort. Heatmap of the differences in the expression of the six hypoxia-linked genes between the low- and high-risk groups in the TCGA (E) data set and the GEO (J) data set. AUC, area under the curve; HCC, hepatocellular carcinoma; TCGA, The Cancer Genome Atlas; GEO, Gene Expression Omnibus; OS, overall survival.

Verification of the prediction model with the GEO data set

To explore the performance of the prediction model, the samples from the GEO data set were adopted as a validation dataset. According to the same formula generated from the TCGA data set, the samples in the GEO data set were similarly stratified into high- (n=39), as well as low-risk (n=42) groups (). Similar to the results obtained from the TCGA data set, the risk plot illustrated that high-risk score patients exhibited shorter survival in contrast with those harbouring low-risk scores (). In addition, Kaplan-Meier data illustrated that high-risk score patients harboured a worse prognosis in contrast with those with low-risk scores (). Likewise, the risk signature AUC values were 0.736 at one year, 0.682 at two years, and 0.633 at three years (). In addition, heatmap analysis also revealed that the expression levels of the six genes were different between the high- and low-risk groups (). The results of the validation analysis demonstrated that our prediction model had satisfactory predictive ability.

Assessment of the independent prognostic role of the hypoxia risk signature score with the TCGA cohort

Univariate along with multivariate regression model were adopted to explore the independent predictive significance of the hypoxia risk model score for OS. In the univariate Cox regression analysis, we identified two independent factors consisting of TNM stage and the risk score (). In the multivariate regression, TNM stage along with the risk score remained independent predictors of OS (). Next, we built a nomogram that incorporated the independent predictors from the multivariate analysis (). This allowed us to estimate the one-, three-, and five-year survival rates of individuals with HCC by plotting a vertical line between the total score axis and each prognosis axis. The ROC curve data illustrated that the nomogram AUC values at one-, three-, and five-years were 0.744, 0.772, and 0.738, respectively (), revealing high accuracy. Moreover, we drew calibration curves to explore the predictive efficiency of the nomogram. As shown in , satisfactory agreement was observed between the predicted and observed outcomes for one-, three-, and five-year OS, indicating that the nomogram had good efficacy in estimating the survival of individuals with HCC.

Figure 4

Independent predictive significance of the hypoxia risk signature score in the TCGA data set. Univariate (A) along with multivariate (B) regression of OS for individuals with HCC in the TCGA data set. (C) The nomogram was created on the basis of two independent predictive factors to predict OS in individuals with HCC at one-, three-, and five years. (D) Time-dependent ROC curves of the model for one-, three-, and five-year OS in HCC to explore the estimation efficiency of the nomogram. (E-G) Calibration plot of the nomogram for one-, three-, and five-year OS in HCC to explore the accuracy of the nomogram. AUC, area under the curve; TCGA, The Cancer Genome Atlas; OS, overall survival; HCC, hepatocellular carcinoma.

GSEA between the high- and low-risk groups

To further investigate associated signaling cascades activated in the high-risk group in contrast with the low-risk group, GSEA was performed. The data of the analysis showed that signaling cascades associated with HCC cell proliferation and metabolism, such as the glycolysis, PI3K/AKT/mTOR, and G2/M checkpoint pathways, were activated in the high-risk group as illustrated in . Besides, these findings were verified by the data from the GEO dataset ().

Figure 5

GSEA between the low- and high-risk groups. (A-C) The genes enriched in the high-risk group were closely related to glycolysis, PI3K-AKT-mTOR signaling, and G2/M checkpoint signaling cascades in the TCGA cohort. (D-F) The data were further verified by the data from the GEO cohort. GSEA, gene set enrichment analysis; TCGA, The Cancer Genome Atlas; GEO, Gene Expression Omnibus.

Differences of immune status between the high- and low-risk groups

To uncover the relationship of the risk signature score with immune status, we calculated the enrichment score of immune cell subclasses and immune-linked roles and cascades using the ssGSEA strategy. As illustrated in , the scores of activated dendritic cells (aDCs), plasmacytoid dendritic cells (pDCs), antigen-presenting cell (APC) coinhibition, APC costimulation, HLA and MHC class I were remarkably higher in the high-risk group in contrast with that in the low-risk group, and these components primarily participate in the antigen presentation process. Moreover, the scores of Th2 cells, macrophages, Tfh cells, Th1 cells, Treg cells and the activation of pathways related to cytokine and cytokine receptors (CCRs), checkpoints, and parainflammation were higher in the high-risk group in contrast with that in the low-risk group, while the scores of mast cells and pathways related to the type II IFN response were elevated in the low-risk group than in the high-risk group. Interestingly, in addition to components of antigen presentation, such as aDCs, APC costimulation, HLA, as well as MHC class I, the macrophage scores were remarkably different between the high- and low-risk groups, as validated with the GEO cohort ().

Figure 6

Comparison of the immune status enrichment scores between the high- and low-risk groups using ssGSEA. The scores of 16 immune cells (A,C) and 13 immune-linked roles or cascades (B,D) in the TCGA along with GEO data sets are exhibited in boxplots. Asterisks represent levels of significance *P<0.05, **P<0.01, ***P<0.001. aDCs, activated dendritic cells; pDCs, plasmacytoid dendritic cells; Treg, regulatory T cells; APC, antigen-presenting cell; CCR, cytokine-cytokine receptor. ssGSEA, single-sample gene set enrichment analysis; TCGA, The Cancer Genome Atlas; GEO, Gene Expression Omnibus.

Analysis of the differences between high- and low-risk groups leading to the development of an immunosuppressive microenvironment

The “cancer-immunity cycle” constitutes a cycle of processes consisting of dynamic antitumor immune responses in 7 steps (27). Dendritic cells capture antigens released during tumour cell apoptosis (Step 1). Captured antigens are presented to T cells by the dendritic cells (Step 2), leading to the priming coupled with triggering of responses of effector T cells against cancer-distinct antigens (Step 3). The stimulated effector T cells transit to (Step 4) and invade the tumor tissue (Step 5), particularly recognize and dock to cancer cells (Step 6), where they kill their target cancer cells (Step 7). Dying cancer cells release antigens, and the cancer-immunity cycle returns to Step 1 and continues (28).

We first downloaded the gene set associated with the cancer-immunity cycle from the Tracking Tumor Immunophenotype web resource (29). Then, we determined the differential expression of genes inversely modulating this cycle between the high- and low-risk groups. As illustrated in , genes participating in the negative modulation of the tumor-immunity cycle were remarkably upregulated in the high-risk group, illustrating that hypoxia can inhibit the activity of this cycle.

Figure 7

Correlation between the hypoxia risk signature score and the immunosuppressive microenvironment. Heatmap of the expression of the genes participating in the inverse modulation of the tumour-immunity cycle between the low- and high-risk groups in the TCGA (A) and GEO (B) cohorts. Relationship of PD-1 expression with the hypoxia risk score in the TCGA (C) and GEO (F) cohorts. Relationship of PD-L1 expression with the hypoxia risk score in the TCGA (D) and GEO (G) cohorts. Correlation of CTLA-4 expression with the hypoxia risk score in the TCGA (E), as well as GEO (H) data sets. Asterisks represent levels of significance *P<0.05, **P<0.01, ***P<0.001. TCGA, The Cancer Genome Atlas; GEO, Gene Expression Omnibus; PD-1, programmed cell death protein 1; PD-L1, programmed death-ligand 1; CTLA-4, cytotoxic T lymphocyte-associated-4.

Furthermore, we assessed the differential expression of genes linked to immune checkpoints between the high- and low-risk groups. illustrates that the expression of PD-1, CTLA-4, along with PD-L1, were positively related with the hypoxia risk score, revealing that upregulated immune checkpoint molecules are involved in the poor prognosis for patients with high-risk score by inducing immune escape. Therefore, these findings indicated that hypoxia facilitates the formation of an immunosuppressive microenvironment in HCC, leading to HCC progression and metastasis.

Discussion

HCC is a frequent malignant tumor with a complicated pathogenesis, as well as high mortality that severely threatens public health worldwide. Although therapeutic strategies have been continuously improved, the prognosis of HCC remains unfavorable, which mainly results from late diagnosis and poor response to current treatment. Therefore, searching for effective diagnostic signatures and treatment targets to improve the prognosis of HCC is necessary.

Hypoxia is a common characteristic of solid tumors, including HCC. Ample evidence has indicated that hypoxia participates in the onset and progress of HCC. Nevertheless, further studies will be required to explore the complex molecular mechanism. Recently, with the advancements in high-throughput sequencing technology combined with bioinformatics analysis, multiple molecular markers have been identified for the construction of prediction models and the development of therapeutic targets, greatly improving the long-term survival of patients.

Herein, we constructed a prediction model including six hypoxia-linked genes, HMOX1, TKTL1, TPI1, ENO2, LDHA, and SLC2A1. HMOX1, as a stress-inducible rate-limiting enzyme, was upregulated in HCC tissues, and high expression promotes HCC cell proliferation via the STAT3 pathway (30). Similarly, TKTL1, as a member of the transketolase gene family, was reported to enhance HCC proliferation because transketolase metabolism controls the pentose phosphate pathways (PPP), which is involved in nucleic acid ribose synthesis in tumor cells (31). TPI1, as a glycolysis-related gene, can be activated by HIF-1α in HCC tissues, and patients exhibiting upregulated expression had shorter OS than those without increased expression (32). Yuan et al. investigated the relationship between circadian genes and glycolysis in HCC cells under hypoxic conditions and found that NPAS2 (one of the core circadian genes) can positively promote the expression of ENO2 by transcriptionally regulating HIF-1α, resulting in the proliferation along with metastasis of HCC cells (33). Sheng et al. documented that the overexpression of LDHA correlated with high metastatic potential in HCC cells, whereas knockdown of the gene strongly inhibited proliferation and metastasis in HCC (34). SLC2A1 (also known as GLUT1), encoding a glucose transporter, was illustrated to be positively linked with the activity of HCC cell glycolysis, and repressing SLC2A1 expression reduced the proliferation along with invasion of HCC cells (35). The findings related to genes other than ENO2 were consistent with the results of our study. However, because there are few relevant studies on the relationship between ENO2 and HCC cells, further research is necessary to identify the role of ENO2.

In routine clinical practice, the TNM staging system is a vital prognostic determinant in HCC. In this study, we also identified the staging system as an independent predictor of OS following univariate along with multivariate Cox regression. However, patients who are at the same stage may harbour a different prognosis, indicating that the current staging system cannot fully reflect intratumor heterogeneity. Hence, it is critical to determine novel molecular signatures to develop a prediction model for the diagnosis and treatment of HCC. For instance, a 21-gene expression assay (termed as Oncotype DX, Genomics Health) is recommended by the National Comprehensive Cancer Network (NCCN) to assess the prognosis of individuals with hormone receptor-positive breast cancer, and this panel has been meaningful for evaluating prognosis and guiding treatment. In our study, a prediction model integrating six hypoxia-linked genes showed satisfactory predictive ability overall. Moreover, a nomogram combining the TNM staging system with the risk signature could efficiently estimate individual survival for individuals with HCC, and this nomogram may enhance the accuracy of determining the prognosis of HCC patients.

According to GSEA, we found that genes associated with the glycolysis pathway, including cascades involved in the growth and metabolization of HCC cells under hypoxic conditions, were remarkably abundant in the high-risk group. It has been well demonstrated that elevated glycolysis provides adequate energy to meet the demands for rapid proliferation of HCC cells in a hypoxic environment (36,37). Moreover, Li et al. found that the circRNA MAT2B promoted HCC progression by enhancing glycolysis, indirectly indicating the key role of glycolysis in the proliferation of HCC cells (38).

The TIME is crucial for the initiation and progression of HCC. According to ssGSEA, we found that the enrichment score of macrophages was remarkably higher in the high-risk group in contrast with that in the low-risk group. This finding is congruent with previous studies showing that hypoxia can induce polarization of macrophages or their migration into the hypoxic TME and that an increased macrophage numbers in the TIME positively correlates with the progress of HCC along with poor prognosis (39-42). In addition, in the context of the cancer-immunity cycle, the results of our study illustrated that genes related with negative modulation of this cycle were remarkably upregulated in high-risk group, illustrating that hypoxia can repress the antitumor immune response. Furthermore, we found that PD-1, CTLA-4 along with PD-L1 were upregulated in the high-risk group, illustrating that hypoxia can induce immune escape. Chen et al. demonstrated that anti-PD1 treatment can boost the antitumor immune response in HCC (43). Similarly, Agdashian et al. reported that anti-CTLA4 treatment can improve the treatment response of patients with advanced-stage HCC (44). Therefore, understanding the key mechanisms by which hypoxia affects immunosuppression and immune escape may facilitate advances in immunotherapy.

However, several limitations to this study need to be acknowledged. First, the construction along with the validation of our prediction model were based on retrospective data from public databases, and further prospective collection of clinical data for validation is necessary. In addition, the sample size of the validation cohort was small and a larger data set is required to validate the findings. Finally, the detailed mechanism of how these identified genes affect HCC cells proliferation and the anti-tumor immune response under hypoxic conditions requires further experimental explanation.

Conclusions

In summary, we used six hypoxia-linked genes to create and validate a prediction model to estimate the OS of individuals with HCC. This hypoxia risk signature has the potential to serve as an independent predictive factor of OS in individuals with HCC that could enable clinicians to estimate the prognosis of individuals with HCC. Besides, the risk signature also reflected the relationship between hypoxic conditions and immune status in the HCC microenvironment, which may help researchers to develop novel immunotherapeutic targets for HCC.