Novel hypoxia features with appealing implications in discriminating the prognosis, immune escape and drug responses of 947 hepatocellular carcinoma patients.

Background: Hypoxia has a profound impact on the development and progression of hepatocellular carcinoma (HCC). This study aimed to explore and elucidate how hypoxia affect prognosis, immune escape and drug responses in HCC.
Methods: HCC-specific hypoxia signatures were identified based on the intersect of differentially expressed genes (DEGs) of GSE41666 and GSE15366. The hypoxia score was calculated using the gene set variation analysis (GSVA) function and validated on GSE18494. We collected five cohorts [The Cancer Genome Atlas (TCGA), GSE14520, GSE39791, GSE36376, GSE57957] for further analysis. First, we analyzed the effect of the hypoxia score on prognosis. Next, we systematically analyzed the potential hypoxia-related immune escape mechanisms and the effect of hypoxia upon immunotherapy. Then, we predicted and screened potential sensitive drugs for HCC patients with high hypoxia levels using machine learning and docking.
Results: We constructed a novel HCC-specific hypoxia score and undertook further analysis in five cohorts (TCGA, GSE14520, GSE39791, GSE36376, GSE57957). We observed that patients with high hypoxia scores exhibited worse overall survival (OS) in TCGA and GSE14520. We also constructed a hypoxia-related nomogram that had good performance in predicting HCC patients' prognosis. Furthermore, patients with lower hypoxia scores had a lower risk of immune escape and thus may benefit from immunotherapy. Finally, we predicted and screened VLX600 as the candidate drug for HCC patients with high hypoxia scores. We further explored and elucidated why VLX600 was more sensitive in HCC patients with high hypoxia than with low hypoxia HCC patients using weighted gene co-expression network analysis (WGCNA). Conclusions: This study provides further evidence of the link between hypoxia and prognosis and immune escape in HCC patients. Moreover, our research screened VLX600 as a potential drug for HCC patients with high hypoxia levels and elucidated the potential mechanism. 2022 Translational Cancer Research. All rights reserved.

Introduction

Hepatocellular carcinoma (HCC) is the fourth leading cause of cancer-related deaths worldwide (1). Unfortunately, patients with advanced HCC have few treatment options, and prognosis is poor. Treatment options for advanced HCC, including surgical resection and non-surgical therapies, are of limited effectiveness (2). Research has shown that metastasis and recurrence rate of HCC are up to 70% after 5 years following curative resection (3). Notwithstanding the merit of immunotherapy in some advanced HCC patients, the survival improvements have been modest (4,5). Therefore, exploring new drug targets or new mechanisms of tumorigenesis is a crucial yet challenging task.

Hypoxia refers to an intrinsic characteristic of solid tumors that is caused by the imbalance between the rate of tumor cell proliferation and the vascular nutrient (6). The association between hypoxia and malignant progression and poor prognosis in HCC has been well documented in the literature (7). Increasing evidence in recent years supports that cancer cells can develop invasive and metastatic features during adaptation to the hypoxic microenvironment (8). Hypoxia can induce HCC cell growth through HMGB1 activation of the TLR9 signaling pathway (9). Dou et al. found that hypoxia can promote the proliferation and migration of HCC cells via upregulating TUFT1 to activate the Ca2+/PI3K/AKT pathway (10). Moreover, recent evidence demonstrates that hypoxia has close relationships with the formation of the immunosuppressive tumor microenvironment. Cui et al. found that hypoxia-inducible lipid droplets promote HCC immune escape from natural killer cells through the interleukin 10/signal transducer and activator of the transcription 3 signaling pathway (11). Ye et al. demonstrated that hypoxia can induce epithelial-mesenchymal transition (EMT) and the establishment of an immunosuppressive tumor microenvironment (TME) via the HIF-1α/CCL20/IDO axis in HCC (12). Although associations between hypoxia and the immunosuppressive microenvironment in HCC have previously been observed in some contexts, the comprehensive interplay between hypoxia and the immunosuppressive microenvironment in HCC is not fully understood. Therefore, further study on the relationship between hypoxia and immune escape in HCC is required in order to develop new therapeutic strategies.

In this study, we analyzed hypoxia-related genes in HCC by using The Cancer Genome Atlas (TCGA) and Gene Expression Omnibus (GEO) databases and constructed a hypoxia score. Then we explored the association of hypoxia with prognosis and the immune microenvironment in HCC. We also predicted and screened potential drugs for patients with high hypoxia levels and discussed the working mechanism of the candidate drugs. These findings may make a meaningful contribution to the development of comprehensive therapeutic strategies for HCC patients. We present the following article in accordance with the MDAR reporting checklist (available at https://tcr.amegroups.com/article/view/10.21037/tcr-22-253/rc).

Methods

HCC data sets downloading and preprocessing

TCGA gene expression data [transcripts per million (TPM) and row counts], copy number data, mutation data, clinical data and sample information of HCC were downloaded from Xena (13). The segment file of copy number was analyzed using GISTIC (14) with default parameters. We retrieved three independent mRNA microarray data sets (GSE41666, GSE15366, and GSE18494) based on hypoxia-treated HCC cells (HepG2) from GEO (https://www.ncbi.nlm.nih.gov/geo/). Data sets and available clinical information on four HCC cohorts including GSE36376 (n=240), GSE14520 (n=225), GSE39791 (n=72), and GSE57957 (n=39) were downloaded from GEO. The raw data of microarray datasets was processed via RMA algorithm background correction, log2 transformation, quantile normalization, and annotation by the Affy R package (15). We also used the IMvigor210 dataset from 623 patients with metastatic urothelial cancer receiving the PD-L1 inhibitor atezolizumab to validate the association between hypoxia and resistance to immunotherapy. The raw transcriptomic and clinical data were retrieved from the IMvigor210 dataset (http://research-pub.gene.com/IMvigor210CoreBiologies) using the IMvigor210CoreBiologies R package (16). We transformed the raw count of IMvigor210 to TPM for further analysis.

Generation of a novel hypoxia gene signature and calculation of hypoxia score

GSE41666 contained the gene expression profles of HepG2 cells exposed to normoxia (21% O2) and hypoxia (0% O2) for 24 h. GSE15366 contained the gene expression profles of HepG2 cells exposed to normoxia (20% O2) and hypoxia (2% O2) for 72 h. GSE18494 collected the gene expression profles of HepG2 cells under normoxic conditions (19% O2, 0 h) and after 4, 8 and 12 h of hypoxia treatment (0.5% O2). We aimed to screen the gene signatures which were stably induced by acute hypoxia and chronic hypoxia. Thus, we calculated the differentially expressed genes (DEGs) in hypoxia and non-hypoxia in GSE41666 (acute hypoxia 24 h) and GSE15366 (chronic hypoxia 72 h). Log-fold changes (logFC) >1 and adjusted P value <0.01 were selected as the threshold. Next, we selected the intersection of up-regulated DEGs in GSE41666 and GSE15366 as the up-regulated hypoxia gene set. The intersection of down-regulated DEGs in GSE41666 and GSE15366 was selected as the down-regulated hypoxia gene set. Gene ontology (GO) enrichment analysis of up- and down-regulated hypoxia gene sets was carried out utilizing the Database for Annotation, Visualization and Integrated Discovery (DAVID) (17). The up hypoxia score was calculated using gene set variation analysis (GSVA) (18) based on the up-regulated hypoxia gene set. The down hypoxia score was calculated using GSVA based on down-regulated hypoxia gene set. The hypoxia score was defined as the up hypoxia score minus the down hypoxia score. To test the reliability of our hypoxia score, we validated our hypoxia score on GSE18494. Finally, we calculated the hypoxia scores for five HCC cohorts (TCGA-LIHC, GSE36376, GSE14520, GSE39791, and GSE57957).

Characteristics of hypoxia-related clinical, molecular, and metabolic features

The clinical data of TCGA and GSE14520 was downloaded from the Xena and GEO websites, respectively. We compared the difference of hypoxia scores in conventional clinical characteristics. To investigate the correlation between hypoxia scores and molecular subclasses of HCC, we used nearest template prediction (NTP) method to make subclass predictions for TCGA and GSE14520 based on subclass-specific signatures. Then we compared the differences in hypoxia scores between different subclasses. Next, we utilized the GSVA method to assess the level of cancer-related hallmarks and metabolism features. The cancer-related hallmarks were downloaded from the Molecular Signatures Database (Msigdb) (19). The metabolism gene sets were obtained from previous research (20). We explored the differences in cancer-related hallmarks and metabolic features between the high hypoxia group and low hypoxia group with respect to TCGA, GSE14520, GSE36376, GSE39791, and GSE57957. Copy number variation (CNV) data on TCGA in the high group and low hypoxia groups were analyzed using the GISTIC with default parameters and visualized using maftools R packages (21). We identified the high hypoxia and low hypoxia peculiar copy number amplification and deletion sites. To further explore how CNV affected hypoxia, we carried out Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analysis for these CNV sites. Mutation data on TCGA in the high hypoxia and low hypoxia groups were analyzed and visualized using maftools R packages.

Survival analysis and building nomograms for hypoxia scores to evaluate the predictive performance of HCC patients’ prognosis

Patients in TCGA and GSE14520 cohorts were stratified into high and low hypoxia groups based on the median value of the hypoxia scores. Kaplan-Meier (K-M) survival analysis was performed to evaluate the overall survival (OS) of patients with high and low hypoxia scores using the survival R package (22). The survival differences were evaluated by the two-sided logrank test. Univariate and multivariate Cox regression analyses were performed to determine whether the hypoxia scores can predict the prognosis independently using survival R package. The nomogram for prognosis prediction was built by integrating the factors that can predict the prognosis independently using the nomogramEx R package (23). The calibration curves of TCGA and GSE14520 were plotted using the rms R package (24).

Association of hypoxia with the immune microenvironment and immunotherapy response

First, the relative levels of distinct immune cell types were estimated using CIBERSORT tools (https://cibersort.stanford.edu) (25). However, CIBERSORT estimates only relative levels of distinct immune cells and does not assess the level of fibroblasts and endothelial cells. Fibroblasts and endothelial cells are the important components of the TME. Thus, we estimated the absolute level of adaptive immune cells and innate immune cells using the GSVA method. The gene signatures of adaptive immune cells and innate immune cells were obtained from The Cancer Immunome Atlas (TCIA) (26). The absolute level of fibroblasts and endothelial cells were assessed using the R package MCP-counter (27). We compared the infiltration of immune cells between the high hypoxia and low hypoxia group using the Wilcoxon test. The list of chemokines was obtained from a previous study (28). The differences in chemokines between the high hypoxia and low hypoxia groups were calculated using the Wilcoxon test and adjusted the P values using the false discovery rate (FDR) method. The potential immunogenomic indicators [mutation load, neoantigen load, homologous recombination defect (HRD) scores, aneuploidy scores, and fraction of loss of heterozygosity (LOH)] were downloaded from the NIH genomic data commons (https://gdc.cancer.gov/about-data/publications#/?groups=&years=&programs=TCGA&order=desc). We collected 76 immunomodulators [10 classical major histocompatibility complex (MHC), 11 non-classical MHC, 25 immune inhibitors, and 47 immune stimulators] from TCIA (26). The correlations between immunomodulator expressions and hypoxia scores were calculated for the TCGA cohort. We calculated T-cell dysfunction scores based on gene sets (TGFB1, CD274, CTLA4, IL10, PDCD1, CD276, HAVCR2, TNFRSF9, LAG3, TIGIT, ICOS) that negatively regulate T-cell function using an ssGSEA algorithm. IFNγ and TGFβ gene sets were obtained from a previous study (29) and the ssGSEA algorithm was used to evaluate the level of IFNγ and TGFβ pathways. The clinical response to immune checkpoint blockade (anti-PD1/CTLA4) was predicted using the Tumor Immune Dysfunction and Exclusion (TIDE, http://tide.dfci.harvard.edu/) tool (30).

Screening sensitive drugs for the high hypoxia group

Expression profile data of human cancer cell lines were obtained from the Broad Institute Cancer Cell Line Encyclopedia (CCLE) project (https://portals.broadinstitute.org/ccle/) (31). Drug sensitivity data on cancer cell lines were obtained from the PRISM repurposing dataset (19Q4, released December 2019, https://depmap.org/portal/prism/) (32). The drug sensitivity data used the area under the dose–response curve [area under the curve (AUC)] values as a measure of drug sensitivity. Lower AUC values indicated higher sensitivity to treatment. We excluded drugs with more than 20% of data missing. Next, we used the K-nearest neighbor algorithm to impute the missing AUC values. Then, we used the pRRophetic R package, which has a built-in ridge regression mode to predict the drug response of patients based on their expression profiles and the AUC value of each drug. The we used the model to predict the drug response for TCGA patients. We further screened the potential drugs for the high hypoxia group. The potential drugs for high hypoxia patients must have a significantly lower AUC than in the low hypoxia group. Thus, we calculated the logFC of predicted differential drug response between the high hypoxia score group (top decile) and low hypoxia score group (bottom decile). We selected logFC <−0.20, and P value <0.05 as the threshold. Moreover, we calculated Spearman correlations between the AUC value and hypoxia to screen drugs that had a negative correlation with hypoxia. We selected correlation <−0.30 and P value <0.05 as the threshold. The selected compounds were studied by docking using Autodock Vina (33) and the binding interaction analyzed using Pymol software. To determine how the sensitive drugs work, we identified the gene co-expression modules using the WGCNA R package (34) and selected the most relevant module. Enrichment analysis using DAVID was carried out using genes in the drug-sensitive relevant module.

Statistical analysis

For comparisons of two groups, we used wilcoxon rank-sum test. For comparisons of more than two groups, we used Kruskal-Wallis tests. The two-sided Fisher’s exact test was performed for the contingency tables. All statistical P values were two-sided. The survival curves of each data set were generated by Kaplan-Meier analysis, and statistically significant differences were determined through the log-rank test. The AUC was calculated using “pROC” package (35). The mutation and copy number of patients were analyzed using the maftools R package (21). All data processing was performed using R 3.6.1.

Ethical statement

The study was conducted in accordance with the Declaration of Helsinki (as revised in 2013).

Results

Identification for a novel HCC-suitable hypoxia signature that includes 74 genes

The overall workflow of this study is shown in . First, we carried out microarray analysis to identify hypoxia-related DEGs with logFC satisfying logFC >1 or logFC <−1 and adjusted P value <0.01 in GSE41666 and GSE15366. A total of 18 negative hypoxia signatures and 56 positive hypoxia signatures were screened (). Functional analysis of these hypoxia signatures strongly indicated that these signatures were hypoxia related (). To further test the hypoxia signatures, we calculated the hypoxia scores on GSE18494. The outcome showed that hypoxic HCC cells had significantly higher hypoxia scores than did the control group and increased with time (). Therefore, the hypoxia scores calculated based on the 74 hypoxia signatures can reflect the hypoxic condition in HCC cells. Thus, we used the 74 hypoxia signatures to assess the hypoxia level for TCGA-LIHC, GSE14520, GSE36376, GSE39791, and GSE57957. Results showed that cancer samples showed a significantly higher hypoxia level than normal samples ().

Figure 1

The workflow of this study.

Figure 2

The 74-gene hypoxia signature indicated hypoxia exposure in HCC cells. (A) Total 18 negative hypoxia signatures and 56 positive hypoxia signatures were screened. (B) The functional enrichment analyses of hypoxia signatures. (C) The hypoxia scores were calculated based on the 74-gene hypoxia signature in the HCC samples with different hypoxia time. (D) The distribution of hypoxia scores in cancer samples and normal samples in 5 cohorts (TCGA-LIHC, GSE14520, GSE36376, GSE39791 and GSE57957). *P<0.05, **P<0.01, ***P<0.001. GO, Gene Ontology; FDR, false discovery rate; HCC, hepatocellular carcinoma; TCGA, The Cancer Genome Atlas.

Associations between hypoxia scores and clinical features

In this section we aimed to analyze the association between hypoxia scores and clinical features. In five cohorts, only TCGA-LIHC and GSE14520 contained clinical information. We ranked the hypoxia scores in TCGA-LIHC and GSE14520 from low to high to explore the associations between hypoxia and clinical features ( and Figure S1). As shown in and Figure S1, patients aged ≤55 years, female and HBV (−) patients had significantly higher hypoxia scores than older, male and HBV (−) patients in the TCGA cohort. However, the age, gender and HBV features did not show significant differences in the GSE14520 cohort. In the GSE14520 cohort, larger tumor size and recurrence indicated higher a hypoxia level, while the hypoxia level did not show significant differences in different recurrence statuses in the TCGA cohort ( and Figure S1). Both the TCGA-LIHC and GSE14520 cohort showed that the later the stage, the higher the hypoxia level. To investigate the correlation between hypoxia and molecular subclasses, we used the NTP method to make subclass predictions based on subclass-specific signatures and compared the difference of hypoxia scores between different subclasses. In the past decade, accumulating robust transcriptome-based subclasses of HCC have been identified, including classifications of Boyault et al. (G1–G6) (36), Chiang et al. (37), Hoshida et al. (S1–S3) (38), Désert et al. classification (39), and Yang et al. (40). Hypoxia scores were found to be higher in Chiang’s unannotated, Désert’s ECM/STEM, Hoshida’s S1, and Yang’s C2 subclasses ( and Figure S2). Notably, these subclasses were demonstrated to have worse prognosis by corresponding studies, which was consistent with their characteristic higher hypoxia scores. Our results indicated that a higher hypoxia level correlated with advanced stage and worse molecular subtypes.

Figure 3

The clinical and biological features associated with hypoxia scores. (A) An overview of the association between hypoxia scores and clinical features of HCC patients in TCGA cohort. (B) An overview of the association between hypoxia scores and clinical features of HCC patients in GSE14520 cohort. (C) An overview of the association between hypoxia scores and HCC molecular subtypes in TCGA. (D) An overview of the association between hypoxia scores and HCC molecular subtypes in GSE14520. (E) The difference of biological features in high hypoxia patients and low hypoxia patients. *P<0.05, **P<0.01, ***P<0.001. HCC, hepatocellular carcinoma; TCGA, The Cancer Genome Atlas; N, normal; AVR-CC, active viral replication chronic carrier.

Associations between the hypoxia scores and molecular characteristics

First, we analyzed the biological hallmarks in the high hypoxia and low hypoxia groups. We found that the hallmark of angiogenesis, EMT, mitotic spindle, and WNT-β catenin signaling, were significantly enhanced in HCC patients with higher hypoxia scores (). The hallmarks of fatty acid metabolism, oxidative phosphorylation, and xenobiotic metabolism were significantly attenuated in HCC patients with higher hypoxia scores (). We then analyzed the distribution of mutations in patients with high hypoxia and low hypoxia scores (). Utilizing the maftools R package, we found that patients with mutations of CTNNB1 and ALB demonstrated significantly lower hypoxia scores (). In order to search for cellular pathways associated with CTNNB1 and ALB, we performed correlation analyses between DEGs of CTNNB1 or ALB mutations and hypoxia-related DEGs. However, no ALB mutation-related DEGs were identified. For CTNNB1, 216 down-regulated and 177 up-regulated genes were identified with logFC >1 or logFC <−1 with an adjusted P value <0.05 (wild type vs. mutation). For hypoxia, 245 up-regulated and 233 down-regulated genes were identified with the same criterion (high hypoxia vs. low hypoxia). We obtained 85 shared up-regulated genes and 90 shared down-regulated genes in both CTNNB1-related DEGs and hypoxia-related DEGs (). KEGG enrichment analysis revealed that shared up-regulated genes were significantly enriched in protein digestion and absorption pathways, the AGE−RAGE signaling pathway in diabetic complications pathway, proteoglycans in cancer pathway, relaxin signaling pathway, and amoebiasis pathway (). Those shared down-regulated genes were significantly enriched in metabolic or biosynthesis-related pathways (). We further analyzed significantly amplified or deleted CNV in patients with high hypoxia () and low hypoxia scores () using GISTIC. We found that the patients with high hypoxia and low hypoxia scores had different CNV inclination (). In patients with high hypoxia scores, 1q21.3, 1q44, 2p24.1, 2p25.1, 2q33.1, 2q33.2, 3q26.2, 3q29, 4p14, 5q35.3, 6p25.2, 7q11.23, 7q21.2, 7q31.2, 8q13.3, 8q22.1, 8q24.21, 11q13.4, 13q34 and 19q13.11 were the characteristic copy number amplification regions (). In patients with high hypoxia scores, 1p36.32, 3p13, 4q22.2, 4q34.3, 8p23.2, 10q23.31, 10q26.2, 12q24.32, 13q14.2, 13q22.2, and 19p13.3 were the characteristic copy number deletion regions (). 1q32.1, In patients with low hypoxia scores,1q32.3, 1q42.3, 3q25.1, 5q35.1, 8q11.1, 8q24.13, 8q24.3, 11q13.2, 13q33.3, 17q23.2, 17q25.1, 19q12 and 20q13.13 were the characteristic copy number amplification regions (). In patients with low hypoxia scores, 1p32.3, 1p36.31, 4q21.23, 4q35.1, 13q14.13, 14q32.33 and 16q23.1 were the characteristic copy number deletion regions (). High hypoxia specific copy number amplification regions significantly enriched in olfactory transduction pathway (). With high hypoxia scores, specific copy number deletion regions were significantly enriched in alcoholic liver disease, pyruvate metabolism, tyrosine metabolism, fatty acid degradation, and glycolysis/gluconeogenesis pathways (). With low hypoxia scores, specific copy number amplification regions did not enrich in any pathway. Low-hypoxia specific copy number deletion regions were significantly enriched in carbohydrate digestion and absorption, pancreatic secretion, complement and coagulation cascades, drug metabolism by cytochrome P450, and metabolism of xenobiotics by cytochrome P450 pathways ().

Figure 4

The mutation features associated with hypoxia scores. (A) Oncoplot of 10 most frequently mutated genes in TCGA HCC patients. (B) The association between mutations and the level of hypoxia. (C) Scatter plot depicting the DEGs between hypoxia (x-axis) and CTNNB1 mutations (y-axis) and each of 58,581 genes. Eighty-five genes were found to be up-regulated in both high hypoxia group and wild type CTNNB1 group (red nodes). Ninety genes were found to be down-regulated in both low hypoxia group and mutation CTNNB1 group (blue nodes). (|logFC| >1 and adjust P≤0.05) (D) Top 5 most significantly enriched gene sets in up-regulated genes (red) and down-regulated genes (blue) (adjusted P<0.05). *P≤0.05, ***P≤0.001. TCGA, The Cancer Genome Atlas; HCC, hepatocellular carcinoma; DEGs, differentially expressed genes.

Figure 5

The copy number variation features associated with hypoxia scores. (A) The significant copy number variation regions in HCC patients with high hypoxia scores. (B) The significant copy number variation regions in HCC patients with low hypoxia scores. (C) The unique copy number amplification and deletion regions of patients with high hypoxia scores (red) and low hypoxia scores (blue). (D) The KEGG enrichment results of high hypoxia-related copy number variations. The copy number amplification regions only enriched in 1 pathway (red bar). The top 5 enriched pathways of copy number deletion regions were displayed (blue bar). (E) The KEGG enrichment results of low hypoxia related copy number variations. The copy number amplification regions did not enrich in any pathways. The top 5 enriched pathways of copy number deletion regions were displayed (blue bar). HCC, hepatocellular carcinoma; KEGG, Kyoto Encyclopedia of Genes and Genomes.

Associations between the hypoxia scores and metabolism features

In the previous section, the outcome of hallmark analysis, mutation analysis, and copy number analysis indicated that metabolism in HCC patients with high hypoxia scores was inhibited. To further explore the metabolism features in a hypoxic microenvironment, we used the gene sets of seven metabolic super-pathways (vitamins and cofactors, lipid, carbohydrate, energy, TCA cycle, amino acids, and nucleotides) published in previous studies for reference (20,41). First, we assessed the level of seven metabolic super-pathways utilizing the GSVA method and compared the differences between high hypoxia and low hypoxia groups. Our outcome showed that all metabolic super-pathways had a lower level in the high hypoxia group than in the low hypoxia group (). Vitamin cofactor, lipid, amino acid, and TCA cycles were significantly down-regulated in the high hypoxia group in most cohorts.

Figure 6

The metabolism features associated with hypoxia scores. (A) The difference of metabolic super-pathways between high hypoxia group and low hypoxia group. (B) Clustering pattern of the seven metabolic super-pathways based on the similarity of the GSVA scores across five cohorts. (C) Univariate Cox regression analysis of metabolic super-pathways in TCGA cohort and GSE14520 cohort. Blue nodes mean HR <1 while red nodes represent HR >1. Only lipid and energy super-pathways were significant in both TCGA cohort and GSE14520 cohort. (D,E) Multivariate Cox regression analysis of hypoxia, lipid and energy super-pathways in TCGA cohort and GSE14520 cohort. Lipid super-pathway did not showed significance in any cohort. Energy showed significance in TCGA cohort but not in GSE14520 cohort. TCA, tricarboxylic acid cycle; GSVA, gene set variation analysis; TCGA, The Cancer Genome Atlas; HR, hazard ratio.

We compared the similarity of different metabolic super-pathways based on sample-level labels and found that vitamins and cofactors, lipid, carbohydrate, and energy pathways formed one tight cluster, whereas integration of the TCA cycle, amino acid, and nucleotide pathways formed another distinct cluster (). Moreover, vitamin and cofactor and lipid pathways formed a tight sub-cluster. There was a high correlation with the level of vitamins and cofactors and lipids. The TCA cycle and amino acids also formed a tight sub-cluster and showed high correlation. To further analyze the effect of metabolic abnormalities on patients’ survival, we carried out univariate Cox regression analysis. We found that lipids and energy were significant protective factors for HCC patients (). However, only energy was a significant protective factor in the TCGA cohort based on multivariate Cox proportional-hazards analysis (). Although energy did not show significance in the GSE14520 cohort, energy exhibited a protective tendency (). Lipids were not significant in both the TCGA cohort and GSE14520 cohort (). These results provide an overall view of the alterations of different metabolic pathways and the effect on survival in the hypoxic microenvironment.

Hypoxia scores can be utilized as an independent prognostic factor in HCC patients

To test the performance of hypoxia scores in predicting prognosis, we drew the time-dependent ROC curves, and the AUCs were more than 0.6 in both TCGA cohort and GSE14520 cohort (). The AUCs indicated the hypoxia scores performed well in predicting prognosis. Next, we explored the relation between hypoxia scores and patients’ survival. According to the median value of the hypoxia scores, 369 patients in the TCGA cohort and 221 patients in the GSE14520 cohort were stratified into high and low hypoxia score groups, respectively. K-M survival analysis suggested that the high hypoxia group had significantly worse OS (). Then, we used univariate and multivariate Cox regression analysis to evaluate whether hypoxia scores can predict the prognosis independently of other traditional clinical characteristics. The results showed that the TNM stage and the hypoxia score were independent prognostic factors for OS both in the TCGA cohort and the GSE14520 cohort (). To accurately predict a certain clinical outcome, we constructed predictive nomograms for the TCGA cohort and GSE14520 cohort (). As the independent factors, hypoxia scores and stage were assigned scores and a total score was calculated by summing up the scores in each individual. The survival probability for the individuals at 1-, 3-, and 5-years was obtained through the function conversion relationship of total scores. We further constructed the calibration plot for internal validation of the nomogram. The calibration plot showed better consistency between the predicted OS outcomes and actual observations both in the TCGA cohort and GSE14520 cohort ().

Figure 7

Evaluation of prognostic value of hypoxia scores. (A) The time dependent ROC curves of hypoxia scores in TCGA cohort. AUCs at 1, 3, and 5 years were used to assess prognostic accuracy. (B) The time dependent ROC curves of hypoxia scores in GSE14520 cohort. AUCs at 1, 3, and 5 years were used to assess prognostic accuracy. (C) The Kaplan-Meier survival analysis of high hypoxia group and low hypoxia group in TCGA cohort. High hypoxia group showed significant worse survival outcome. (D) The Kaplan-Meier survival analysis of high hypoxia group and low hypoxia group in GSE14520 cohort. High hypoxia group showed significant worse survival outcome. (E) Nomogram of the hypoxia score to predict 1-, 3- or 5-year OS in TCGA cohort using hypoxia score and TNM stage. (F) Nomogram of the hypoxia score to predict 1-, 3- or 5-year OS in GSE14520 cohort using hypoxia score and TNM stage. (G) Calibration curves of the nomogram for predicting the probability of OS at 1, 3, and 5 years in TCGA cohort. (H) Calibration curves of the nomogram for predicting the probability of OS at 1, 3, and 5 years in GSE14520 cohort. AUC, area under the curve; TCGA, The Cancer Genome Atlas; LIHC, liver hepatocellular carcinoma; ROC, receiver operating characteristic; OS, overall survival.

Table 1

The univariate and multivariate Cox regression analysis outcome in TCGA cohort

Variables	Univariate			Multivariate
	HR (95% CI)	P value		HR (95% CI)	P value
Age	1.01 (1–1.03)	0.08948		Not include
Gender
Male	Reference			Not include
Female	1.21 (0.85–1.74)	0.28596
BMI	0.97 (0.94–1.01)	0.1179		Not include
Alcohol
No	Reference			Not include
Yes	1.05 (0.72–1.53)	0.80107
Fibrosis ISHAK score
0	Reference			Not include
1 & 2	1.08 (0.7–1.66)	0.63632
3 & 4	0.68 (0.28–1.66)	0.40065
5	0.75 (0.18–3.17)	0.6967
6	0.72 (0.39–1.33)	0.28791
Child
A	Reference			Not include
B	1.63 (0.77–3.44)	0.20077
C	2.14 (0.29–15.55)	0.45356
TNM stage
Stage I	Reference			Reference
Stage II	1.49 (0.91–2.45)	0.11347		1.22 (0.74–2.02)	0.42927
Stage III	2.8 (1.83–4.3)	0.0000023		2.17 (1.39–3.39)	0.00061
Stage IV	5.8 (1.79–18.86)	0.00345		4.02 (1.22–13.25)	0.02214
Recurrence
New primary tumor	Reference			Not include
Locoregional recurrence	1.05 (0.36–3.01)	0.93274
Intrahepatic recurrence	0.48 (0.17–1.38)	0.17464
Extrahepatic recurrence	0.71 (0.24–2.14)	0.54772
Vascular invasion
No	Reference			Not include
Yes	1.34 (0.88–2.04)	0.16919
Recurrence
No	Reference			Reference
Yes	113.68 (15.8–818.06)	2.59E-06		93.61 (12.99–674.74)	1.00E-05
Hypoxia scores	1.9 (1.38–2.62)	0.00008		1.44 (1.02–2.04)	0.03631

TCGA, The Cancer Genome Atlas.

Table 2

The univariate and multivariate Cox regression analysis outcome in GSE14520 cohort

Variables	Univariate			Multivariate
	HR (95% CI)	P value		HR (95% CI)	P value
Age	0.99 (0.97–1.01)	0.40489		Not include
Gender
Male	Reference			Not include
Female	0.59 (0.28–1.22)	0.15329
HBV
Normal	Reference			Not include
CC	1.06 (0.26–4.35)	0.93327
AVR-CC	1.42 (0.34–6)	0.63265
ALT
Low	Reference			Not include
High	1.08 (0.7–1.66)	0.72653
Tumor size
Small	Reference			Reference
Large	1.92 (1.25–2.96)	0.0028		1.08 (0.58–2.02)	0.80389
Cirrhosis
No	Reference			Reference
Yes	4.62 (1.14–18.8)	0.03241		2.82 (0.68–11.69)	0.15257
TNM stage
Stage I	Reference			Reference
Stage II	2.15 (1.24–3.73)	0.00658		1.42 (0.81–2.51)	0.22329
Stage III	5.22 (2.98–9.14)	8.03E-09		2.7 (1.36–5.34)	0.00444
AFP
Low	Reference			Reference
High	1.63 (1.06–2.5)	0.02505		1.15 (0.71–1.86)	0.58331
Recurrence
No	Reference			Reference
Yes	113.68 (15.8–818.06)	2.59E-06		93.61 (12.99–674.74)	1.00E-05
Hypoxia scores	2.29 (1.53–3.44)	6.00E-05		1.77 (1.14–2.74)	0.01037

HBV, hepatitis B virus; CC, chronic carrier; AVR-CC, active viral replication chronic carrier; ALT, alanine aminotransferase; AFP, alpha-fetoprotein..

Associations between hypoxia scores and potential extrinsic immune escape mechanisms

The extrinsic immune escape mechanism mainly composed by four major aspects: lack of immune cells, presence of immunoinhibitory cells [such as type 2 macrophages, regulatory T cells (Tregs) and myeloid derived suppressor cells (MDSCs)], high concentrations of immunoinhibitory cytokines, and fibrosis (42). Based on the results of Cibersort, we found that high hypoxia score groups had a significantly higher proportion of M0 macrophages and lower proportion of M2 macrophages and resting mast cells (). We further used the GSVA method to assess the level of immune cells. The results showed that the high hypoxia group score had a significantly higher level of central memory CD4 T cells, activated CD4 T cells, and MDSC (). We also used the MCPcounter R package to assess the levels of endothelial cells and fibroblasts. We found that high hypoxia score group had a significantly higher level of fibroblasts (). These findings suggest that patients with high hypoxia scores might be in a state with anti-tumor characteristics. However, the high hypoxia group not only had abundant number of active innate and adaptive immune cells but also immunosuppressive cells, such as MDSC. The infiltration of immunosuppressive cells suggests a role in immune escape. In addition, the high hypoxia group also had a higher level of stromal cells than in the low hypoxia group. The abundant stromal cells may also contribute to immune escape.

Figure 8

Association of hypoxia with potential extrinsic immune escape mechanisms of LIHC. (A) The comparison of proportion of immune cells estimated by cibersort between high hypoxia group and low hypoxia group in 5 cohorts. (B) The comparison of absolute level of immune cells estimated by ssGSEA method between high hypoxia group and low hypoxia group in 5 cohorts. (C) The comparison of absolute level of stromal cells using MCPcounter method. (D) The comparison of mRNA expression of chemokines and receptors in high hypoxia group compared with low hypoxia group in 5 cohort. (E) The comparison of mRNA expression of interleukins and receptors in high hypoxia group compared with low hypoxia group in 5 cohort. (F) The comparison of mRNA expression of interferons and receptors in high hypoxia group compared with low hypoxia group. (G) The comparison of mRNA expression of other important cytokines in high hypoxia group compared with low hypoxia group. (A-G) Red block represents up-regulated in high hypoxia group and blue block represents down-regulated in high hypoxia group. The black boxes represent significant difference (adjust P<0.05). Grey block indicates not available in the cohort. LIHC, liver hepatocellular carcinoma; ssGSEA, single sample Gene Set Enrichment Analysis algorithm.

We further analyzed the levels of cytokines in the high hypoxia group and low hypoxia group. We found CCL16, ARG1, IL6R, and FAS were significantly up-regulated in the high hypoxia group (). In contrast, CCL20, CCR6, CXCL6, IL17RD, IL4R, IL8, IFNGR2, EPOR, PDGFA, VEGFA, and VEGFB were significantly down-regulated in the high hypoxia group.

Association of hypoxia scores with potential intrinsic immune escape mechanisms and tumor immunogenicity in HCC

We analyzed the hypoxia-associated potential of intrinsic immune escape mechanisms in HCC patients. Intrinsic immune escape indicates that tumor cells directly develop their ability for immune escape. We analyzed the intrinsic immune escape mechanism from two aspects: tumor immunogenicity and expression of immune checkpoint molecules.

First, we compared potential tumor immunogenicity-related factors: mutation load, neoantigen load, and chromosomal instability level. These factors are the main source of tumor antigens. Patients with high hypoxia levels had significantly higher levels of silent mutation load, although there were no differences in non-silent mutation load and neoantigen load between high hypoxia score patients and low hypoxia score patients (). We observed that high hypoxia score patients had a significantly higher level of chromosomal instability indices (HRD scores, aneuploidy scores, fraction of LOH and CNV load) compared with low hypoxia score patients (). Significantly, the CNV load had a significantly higher level in both focal level and arm level in higher hypoxia score patients (Figure S3). The higher level of chromosomal instability may contribute to the high immunogenicity in the high hypoxia score group. Meanwhile, we observed that high hypoxia patients had significant higher level of cytolytic scores ().

Figure 9

Association of hypoxia with potential intrinsic immune escape mechanisms of LIHC. Comparison of (A) silent mutation load, (B) non-silent mutation load, (C) neoantigen load, (D) HRD scores, (E) aneuploidy score and (F) fraction of LOH altered (G) cytolytic scores between high hypoxia group and low hypoxia group in TCGA cohort. Comparison of (H) MHC class-I molecules, (I) MHC class-II molecules, (J) immunoinhibitors and (K) immunostimulators between high hypoxia group and low hypoxia group in 5 cohorts. (G-J) Red block represents significant up-regulated in high hypoxia score group and blue block represents significant down-regulated in high hypoxia score group. The black boxes represent significant different. (L) The correlation among hypoxia scores, cytolytic scores, absolute level of cytolytic T cells (activated CD4 and CD8 T cells) and mRNA expression of immunoinhibitors. Red dots represent positive correlation and blue dots represent negative correlation. The black boxes represent significant difference (P<0.05). *P<0.05. **P<0.01. ***P<0.001. ns, not significant. LIHC, liver hepatocellular carcinoma; HRD, homologous recombination defect; LOH, loss of heterozygosity; TCGA, The Cancer Genome Atlas; MHC, major histocompatibility complex.

Next, we compared the expression levels of MHC molecules, immune stimulators, and immune inhibitors. Our results showed that TAP1, an MHC class I molecule, were significantly up-regulated in high hypoxia groups in most HCC cohorts (). MHC class II molecules did not significantly alter in most cohorts (). As for immune inhibitors, we found KDR was significantly down-regulated in most cohorts (). IL6R, an immunostimulatory agent, was significantly down-regulated in most cohorts (). Consistently, we observed that hypoxia scores showed significant positive correlation with most immunoinhibitors. Meanwhile, hypoxia scores also positive correlated with activated CD4 T cell and cytolytic scores ().

Association of hypoxia scores with the response to immunotherapy

We explored the relationship between hypoxia scores and immunotherapy. We observed that higher hypoxia scores were significantly correlated with the infiltration of anti-tumor immune cells (such as activated CD4 T cells and central memory CD4 T cells). This implies that patients with high hypoxia scores may have an activated immune state. Meanwhile, high hypoxia score patients had a significantly higher level of fibroblasts (). Research has shown that fibroblasts may have negative effects on anti-tumor immune responses and immunotherapy (43). Patients with higher hypoxia scores showed significantly higher T-cell dysfunction in most cohorts (). This indicated that high hypoxia patients tend to have T-cell exhaustion. These patients’ T-cell function needs to be reinvigorated by immune checkpoint therapy. Furthermore, we compared the level of TGF-β and IFN-γ between high hypoxia group and low hypoxia group. We observed a significantly higher level of TGF-β in high hypoxia patients in all cohorts () but IFN-γ did not exhibit significant differences in most cohorts (). TGF-β has been verified as inducing tumor progression, metastasis, and poor responses to cancer immunotherapy (44). These results suggested high hypoxia patients may not have a good response to immunotherapy. To validate our inference, we used the TIDE tool to predict the response of patients to immune checkpoint blockade (CTLA4 and PD1 therapy). The results showed that high hypoxia patients had significantly higher TIDE scores in almost all cohorts (). This indicated that high hypoxia patients may have a poor response to immunotherapy. In addition, we further validated our conjecture using IMvigor210 in 384 patients with metastatic urothelial cancer receiving PD-L1 inhibitors with atezolizumab. First, we utilized an alluvial diagram to show distribution changes of vital status, overall response, and binary response in individuals in the IMvigor210 cohort (). The K-M survival analysis showed that patients with low hypoxia scores had significant clinical benefits and markedly prolonged survival when compared with patients with high hypoxia scores (). Although the distribution of patients with different overall responses did not show significant differences (Fisher’s exact test, P=0.05251, ), the hypoxia scores of patients with complete response (CR) and partial response (PR) were significantly lower than patients with stable disease (SD) and progressive disease (PD) (Kruskal–Wallis test, P=0.04, ). The proportion of immunotherapy responders was significantly higher in patients with low hypoxia scores compared with patients with high hypoxia scores (Fisher’s exact test, P=0.0185, ). Patients with CR/PR showed significantly lower hypoxia scores than SD/PD patients (Wilcoxon test, P=0.006, ). We also validated the role of hypoxia scores in determining different immune phenotypes in the IMvigor210 cohort, which contained complete information on immune phenotypes. We observed that the patients with desert phenotype remarkably reduced in low hypoxia group (Fisher’s exact test, P=0.02679, ). Meanwhile, the excluded phenotype significantly accumulated in the low hypoxia group. However, the hypoxia scores did not show significant differences among the three phenotypes (Kruskal–Wallis test, P=0.42, ). To further determine whether hypoxia scores can predict the response to immunotherapy, we carried out ROC analysis. The AUC indicated that our hypoxia scores had certain predictive ability independently and can improve the prediction performance of tumor mutation burden (TMB) (). Our results strongly indicated that hypoxia scores were significantly associated with immune escape and could contribute to the prediction of response to immunotherapy.

Figure 10

Association of hypoxia with immunotherapy of LIHC. Comparison of (A) T cell dysfunction scores, (B) TGF-β scores, (C) IFN-γ scores, (D) TIDE scores between high hypoxia group and low hypoxia group in 5 cohorts. (E) Alluvial diagram showing the changes of hypoxia scores, vital status, overall response and binary response in IMvigor210 cohort. (F) Kaplan-Meier survival analysis of high hypoxia group and low hypoxia group in IMvigor210 cohort. Patients with high hypoxia showed significant worse survival than patients with low hypoxia. (G) The proportion of overall response in high and low hypoxia groups in IMvigor210 cohort. The statistical difference was measured with the Fisher’s exact test (P=0.05251). (H) Differences in hypoxia scores among different types of overall response in IMvigor210 cohort. The Kruskal-Wallis test was used to compare the statistical difference between different molecular subtypes (Kruskal-Wallis test, P=0.04). (I) The proportion of binary response in high and low hypoxia groups in IMvigor210 cohort. The statistical difference was measured with the Fisher’s exact test (P=0.0185). (J) Differences in hypoxia scores between different types of binary response in IMvigor210 cohort. The Wilcoxon test was used to compare the statistical difference between different molecular subtypes (Wilcoxon test, P=0.006). (K) The proportion of immune microenvironment type in high and low hypoxia groups in IMvigor210 cohort. The statistical difference was measured with the Fisher’s exact test (P=0.02679). (L) Differences in hypoxia scores among different types of immune microenvironment in IMvigor210 cohort. The Kruskal-Wallis test was used to compare the statistical difference between different types of immune microenvironment (Kruskal-Wallis test, P=0.42). (M) ROC curves measuring the predictive value of the hypoxia scores, TMB, and combination of them in IMvigor210 cohort. LIHC, liver hepatocellular carcinoma; ns, not significant; TGF-β, transforming growth factor-β; PD, progressive disease; SD, stable disease; PR, partial response; CR, complete response; ROC, receiver operating characteristic curve; TMB, tumor mutation burden. *, P<0.05; **, P<0.01; ***, P<0.001.

Screening potential drugs for high hypoxia score HCC patients

Our results indicated that high hypoxia HCC patients tend have immune escape and a poor response to immunotherapy as well as a worse prognosis. Thus, screening drugs for high hypoxia HCC patients could have a profound impact on prognosis and immunotherapy. In this study, we constructed a drug response prediction model using the pRRophetic R package based on expression profile data and drug sensitivity data of CCLE. Next, we predicted the drug response for TCGA, GSE36376, GSE14520, GSE39791, and GSE57957 in five HCC cohorts. We screened drugs in which Spearman correlations between the AUC value and hypoxia scores were less than −0.3. Meanwhile, we conducted a differential drug response analysis between the high hypoxia score group (top decile) and the low hypoxia score group (bottom decile) to identify the potential drugs for the high hypoxia score group (log2FC <−0.2). A total of 16 compounds in five cohorts were screened (). Notably, VLX600 had lower estimated AUC values in the high hypoxia score group and a significant negative correlation with the hypoxia scores in most cohorts (GSE14520, GSE36376, and GSE39791). A compound with a lower AUC suggested that this compound may have high drug sensitivity. These results indicated that the VLX600 may exhibit high drug sensitivity in high hypoxia score patients. To further evaluate the association between VLX600 and hypoxia score, we focused on the mechanisms and effects of VLX600. VLX600 is a novel iron chelator that can disrupt intracellular iron metabolism, leading to inhibition of mitochondrial respiration and bioenergetic catastrophe with resultant tumor cell death (45). VLX600 showed enhanced cytotoxic activity under conditions of nutrient starvation and tumor growth inhibition in vivo (46). Therefore, VLX600 showed good performance in vitro and in silico evidence and may be a promising treatment.

Figure 11

Screening of potential drugs with higher drug sensitivity in high hypoxia patients. (A) The potential sensitive drugs for patients with high hypoxia scores in 5 cohort. Blue block represents sensitive in high hypoxia patients, while red block represents not sensitive in high hypoxia patients. Black box represents the sensitive value is significant different in high hypoxia group compared with low hypoxia group. *, the sensitive value is significant correlated with hypoxia scores. (B) The binding mode of VLX600 in USP14 binding pocket. VLX600 existed hydrogen bond with ASN-109. (C) The protein-protein interaction network of VLX600 sensitive related genes. The size of the node positively correlated with the degree of the node. (D) The KEGG enrichment analysis outcome of VLX600 sensitive related genes. TCGA, The Cancer Genome Atlas; PD, progressive disease; SD, stable disease; PR, partial response; CR, complete response; KEGG, Kyoto Encyclopedia of Genes and Genomes.

The annotation information downloaded from the PRISM website showed that USP14 was the target of VLX600. To explore how VLX600 binds with its target, we carried out molecule docking of VLX600 and USP14 using Autodock vina. Based on the binding mode of VLX1570 (an analog of VLX600) with USP14 in a previous study (47), we roughly identified the binding position of VLX600 with USP14. The docking outcome showed that VLX600 binds with USP14 in a similar position as compared with VLX1570 (). In particular, VLX600 binds with ASN109 via hydrogen bonds. This indicated that ASN109 may play an important role in the binding of VLX600 and USP14. To further explore how VLX600 works, we used the WGCNA method to screen gene modules that were significantly correlated with the AUC of VLX600. The blue module (GSE14520), brown module (GSE36376), and turquoise module (GSE39791) were most positively correlated with the AUC of VLX600 (Figure S4), and there were no shared genes among these modules. The yellow module (GSE14520), green module (GSE36376), and red module (GSE39791) were most negatively correlated with the AUC of VLX600 (Figure S4). A total of 516 genes were selected as the overlap for three modules (). The 516 genes were enriched in 44 pathways (). Notably, most of these pathways were associated with metabolism. In general, VLX600, which showed robust in vitro and in silico evidence, was considered to have the most promising therapeutic potential in HCC patients with high hypoxia scores.

Discussion

Accumulating evidence suggests that hypoxia has an excellent effect in tumorigenesis and cancer treatment (48,49). However, how hypoxia affects prognosis, immune escape, and drug responses in HCC is still poorly understood. This study aimed to analyze the relationship between hypoxia of HCC and the effects of prognosis, immune escape, and drug responses in five different cohorts. Initially, we screened hypoxia signatures and calculated hypoxia scores for 947 HCC patients in five cohorts. Next, we conducted an in-depth analysis of the association between hypoxia scores and molecular subtypes, molecular features, and metabolism pattern. We observed that hypoxia scores were significantly positively correlated with advanced stage, angiogenesis, and EMT but negatively correlated with metabolism pathways. The hypoxia scores showed significant association with molecular subtypes. We identified also identified significant genomics alterations related to hypoxia. Then, the survival analysis outcome showed that hypoxia scores were significantly correlated with patients’ prognosis. Furthermore, the nomogram and calibration curves showed that hypoxia scores combining stage perform well in predicting patients’ prognosis. Through analyzing the immune factors systematically, we observed that high hypoxia score patients had a significantly higher infiltration of CD4+ T cells and MDSC. Moreover, T cells in high hypoxia TME showed a significantly higher level of dysfunction. HCC cells in high hypoxia conditions may sculpt develop their immune escape ability by enhancing genomic instability. Additionally, high hypoxia attenuated sensitivity to immunotherapy. Eventually, we screened VLX600 as the potential drug for high hypoxia HCC patients.

In this study, we observed that HCC in the hypoxic microenvironment inhibited metabolism and enhanced the ability of invasion and metastasis of tumor cells. Amino acid metabolism, carbohydrate metabolism, energy, lipids, nucleotides, tricarboxylic acid cycle, and vitamin and cofactor metabolism were all inhibited in the hypoxic microenvironment in HCC. In addition, HCC patients with high hypoxia scores had a significantly higher level of angiogenesis, EMT, and the Wnt/β-catenin pathway. Research has demonstrated that intra-tumoral hypoxia is caused by an imbalance of O2 availability due to insufficient blood supply from poor tumor vasculature and increased O2 consumption by metabolically active HCC cells (50). Reduced blood supply and high metabolic consumption may lead to lack of nutrients and O2, inhibiting metabolism. This appears to penalize the survival of HCC cells. However, HCC patients with high hypoxia scores showed worse survival and a more advanced stage. A previous study showed that the hypoxic microenvironment can trigger the expression of HIF1A, which can stimulate Wnt/β-catenin signaling in HCC (51). This will also enhance the EMT, invasion, and metastasis of HCC (52). This was in accordance with our outcome. Altogether, we hypothesize that hypoxia weakens HCC cell metabolism and promotes the survival, invasion, and metastasis of HCC cells, eventually leading to poor survival for patients.

Associations between hypoxia and the tumor immune microenvironment have previously been observed in some contexts (53,54). How hypoxia affects immune escape in HCC has not been well characterized. In this work, we systematically explored how hypoxia affected immune escape for HCC patients from extrinsic and intrinsic viewpoints. In extrinsic immune escape, we observed some anti-tumor characteristics [highly adaptive immune cells such as activated CD4 T cells, higher cytolytic scores, and higher expression level of TAP1 (an MHC I class molecule)] that may contribute to better survival outcome in the high hypoxia group; however, both a previous study and the current outcome have shown that HCC patients with high hypoxia had a worse survival outcome (55). We found that the high hypoxia group had more infiltration of immunosuppressive cells [tumor associated macrophages (TAMs) and MDSC]. Meanwhile, immunotherapy analysis showed that high hypoxia groups had higher T-cell dysfunction scores and TGF-β scores. These immunoinhibitory factors might contribute to the extrinsic immune escape and lead to worse survival.

From the intrinsic viewpoint, although non-silent mutation load and neo-antigen load were not significantly different between the high hypoxia group and low hypoxia group, the high hypoxia group showed a significantly higher silent mutation load, HRD scores, aneuploidy scores, fraction of LOH altered, and CNV load. These outcomes indicated that the high hypoxia group have greater chromosomal instability, which may contribute to producing more tumor antigen. Meanwhile, the high hypoxia group had a significantly higher level of MHC molecules, CD4 T cells, and cytolytic score. Hypoxia scores were also positively correlated with cytolytic activity. These results indicate an activated immune state in the high hypoxia group. However, we also found that both the hypoxia scores and cytolytic activity were significant correlated with immunoinhibitor and immune checkpoint molecules (such as TIGIT, CTLA4, and CD96). This indicated that T-cell dysfunction and immunosuppression may be related in high hypoxia. The significantly higher level of T-cell dysfunction scores in the high hypoxia group is consistent with this finding. The higher expression level of immune checkpoints and more T-cell exhaustion contribute to immune escape and decrease survival.

A previous study showed that a high level of hypoxia in HCC indicated a poor immunotherapy response (54). Our outcome was consistent with this finding. In this study, we analyzed the potential mechanisms of how hypoxia attenuates immunotherapy. We observed that high hypoxia patients had significantly higher TGF-β scores, higher T-cell dysfunction scores, and more fibroblasts. Moreover, high hypoxia scores were significantly correlated with a high level of immuoinhibitors. These features may explain the unsatisfactory immunotherapeutic response to some extent. These outcomes further demonstrated that high hypoxia patients with increased T-cell infiltration, high cytolytic activity, and active anti-tumor immune status are more likely to receive immunotherapy but do not show high responsiveness because of increased immune escape and T-cell dysfunction. We further explored the association between hypoxia and immunotherapy. We found that patients with high hypoxia had a higher proportion of CR/PR. Also, higher hypoxia scores were strikingly associated with desert immune phenotypes in which it was difficult for immunotherapy to exert an antitumor effect.

We observed that high hypoxia patients had worse survival outcome than low hypoxia patients. Therefore, screening drugs in high hypoxia patients could have a profound impact on prognosis and treatment. Our prediction outcome showed VLX600 had significant sensitivity in high hypoxia patients in most cohorts. VLX600 is a mitochondrial oxidative phosphorylation (OxPhos) inhibitor which the primary mechanism is inhibition of Cox 1, 2, and 4 of the electron transport chain (46). Notably, VLX600 showed enhanced cytotoxic activity under conditions of nutrient starvation and displayed tumors growth inhibition in vivo (46). Coincidentally, our outcome showed that patients with high hypoxia had significant inhibition of metabolism. This supported that VLX600 may be more sensitive for high hypoxia patients on the other side. Moreover, we observed that the gene modules that were positively correlated with the AUC value of VLX600 were significantly enriched in metabolism-related pathways. This indicated that a high expression level of genes in these modules may lead to active metabolism related pathways and a high AUC value for VLX600. Incidentally, a high AUC value represents low sensitivity of VLX600. Consequently, we inferred that active metabolism may attenuate the sensitivity of VLX600. We were curious as to why VLX600 was more sensitive in an inactive metabolism microenvironment. We inferred that VLX600 can inhibit the enzymes Cox I, II, and IV, which are the key mitochondrial electron transport enzymes. This limits mitochondrial electron transport and reduces the production of ATP. This will result in low ATP/ADP ratio and contribute to aerobic glycolysis. Finally, this will lead to the Warburg phenotype and contribute to proliferation of cancer cells that may be resistant to VLX600. Cancer cells lack the materials for aerobic glycolysis in conditions of nutrient starvation, leading to a lower level of aerobic glycolysis and Warburg phenotype. Also, this may be unfavorable for the survival and proliferation of cancer cells. HCC with a high hypoxia level often lacks blood vessels to provide nutrients and oxygen. This may explain why VLX600 was more sensitive in the high hypoxia microenvironment.

Conclusions

This study provides further evidence of the link between hypoxia and prognosis and immune escape in HCC patients. Moreover, our research screened VLX600 as potential sensitive drug for HCC patients with high hypoxia level and elucidate the potential mechanism.

The article’s supplementary files as