A Novel hepatocellular carcinoma specific hypoxic related signature for predicting prognosis and therapeutic responses

Hypoxia is an important feature of the tumor microenvironment(TME) and is closely associated with cancer metastasis, immune evasion, and drug resistance. However, the precise role of hypoxia in hepatocellular carcinoma(HCC), as well as its influence on the TME, and drug sensitivity remains unclear. We found the excellent survival prediction value of Hypoxia_DEGs_Score model. In hypoxic HCC, somatic mutation, copy number variation, and DNA methylation were closely related to hypoxic changes and affected tumorigenesis, progression, metastasis, and drug resistance. In HCC, aggravated hypoxic stress was found to be accompanied by an immune exclusion phenotype and increased infiltration of immunosuppressive cells. In the validation cohort, patients with high Hypoxia_DEGs_Score were found to have worse immunotherapeutic outcomes and prognoses, and may benefit from drugs against cell cycle signaling pathways rather than those inhibiting the PI3K/mTOR pathway. Hypoxia_DEGs_Score has an excellent predictive capability of changes in the TME, the efficacy of immunotherapy, and the response of drugs. Therefore, Hypoxia_DEGs_Score can help develop personalized immunotherapy regimens and improve the prognosis of HCC patients.

Introduction

Hepatocellular carcinoma(HCC), an aggressive malignancy with a dismal prognosis, is the fourth most common cause of cancer-related deaths worldwide (1). Although the emerging immunotherapy and targeted therapy approaches offer certain benefits, the patients who benefit still account for only a small subset of the total patient population (2, 3). Even if patients have similar clinicopathological characteristics and treatment options, they may still have a completely different prognosis because of the heterogeneity of HCC (4). Therefore, it is imperative to find a signature with high predictive value to guide treatment and improve the prognosis.

Hypoxia is a status of imbalance between the tumor’s oxygen demand and the circulating oxygen supply, mainly because the growth rate of blood vessels cannot keep up with the oxygen demand of the tumor (5–7). To adapt to this status, the tumor will derive a new adaptive dynamic balance to ensure its survival and development (8). Surprisingly, such extreme conditions stimulate the formation of blood vessels, promote invasion and metastasis, increase the instability of the genome, and impair the antitumoral immune response (7, 9). Hypoxia is often particularly obvious and common in advanced solid tumors and is highly correlated with rapid tumor progression, drug resistance, and poor prognosis (8).

The changes in the tumor microenvironment caused by hypoxia involve multiple levels. Hypoxia activates the AKT and ERK pathways through different mechanisms to induce the epithelial-to-mesenchymal transition (EMT) and expression of matrix metalloproteinases (MMPs), which promote the invasion and metastasis of HCC (10). In patients with HCC, the resistance of sorafenib is considered to be related to the hypoxia caused by its anti-angiogenetic effect (11). After reviewing recent research on the characteristics of hypoxia across different cancer types (12–15), we selected 15 genes in the study of Buffa et al. to represent the characteristic genes of hypoxia (12). It is well known that hypoxic stress can drive changes in multiple biological pathways, playing a key role in the occurrence and development of HCC.

Although there have been many advances in HCC and hypoxia research, there are still certain difficulties in the evaluation of the hypoxic status. Direct and accurate measurement of HCC hypoxic status and oxygen amount is still difficult to implement, and hypoxia-inducible factor-1α(HIF-1α) expression cannot accurately reflect all the phenotypes of hypoxia (16). Therefore, it is crucial to find a robust signature that reflects the hypoxic status of HCC to guide treatment decisions.

This study used four HCC datasets for research, encompassing the TCGA-LIHC, ICGC-LIRI-JP, Fudan Zhongshan, and Tongji cohorts, and investigated their hypoxia status. We found that hypoxia was not only correlated with immune cell infiltration, but with genetic instability and epigenetic modifications in HCC. We then constructed a hypoxia score model based on significantly different genes to quantify the hypoxic stress of each sample, which enabled us to assess hypoxia-related changes of tumor biology. Finally, we validated the survival prediction value of Hypoxia_DEGs_Score and confirmed its ability to guide targeted therapy and immunotherapy.

Materials and methods

Dataset collection and preprocessing

The workflow diagram of our study is shown in Figure S1A. The RNA-seq and clinicopathological data of HCC samples were obtained from four different institutions: TCGA (17), ICGC (18), Zhongshan Hospital (19), and Tongji Hospital(this study). The cohorts of immunotherapy interventions were further included to evaluate the relationship between the Hypoxia_DEGs_Score and the benefits of immunotherapy. The patients in the IMvigor210 cohort had advanced urinary tract transitional cell carcinoma and had received anti-PD-L1 therapy (20). GSE78220 is a dataset of an interventional regimen in metastatic melanoma with anti-PD-1 antibodies (21). The basic information of the datasets is summarized in Table S1.

Unsupervised clustering for 15 characteristic hypoxia-related genes

A total of 15 characteristic hypoxia-related genes (ACOT7, ADM, ALDOA, CDKN3, ENO1, LDHA, MIF, MRPS17, NDRG1, P4HA1, PGAM1, SLC2A1, TPI1, TUBB6, and VEGFA) were selected based on the exploration of gene function and co-expression patterns in vivo, after cross-comparison of signature studies related to hypoxic stress (12–15). Unsupervised clustering analysis was used to distinguish the HCC samples according to their hypoxic status. The consensus clustering algorithm in the ConsensusClusterPlus package was used for unsupervised clustering (22, 23).

Gene set variation analysis

To evaluate signature level discrepancies between different hypoxic statuses and states of diverse functions in individual HCC samples, we applied the GSVA package to assess the differences in various functional gene sets, including Kyoto Encyclopedia of Genes and Genomes (KEGG), Gene ontology(GO), HALLMARK, and the 15 hypoxia hallmark genes (24).

Evaluation of immune function status and immune cell infiltration

The TIDE algorithm was used to evaluate the immune function status (25). The CIBERSORT algorithm and the ssGSEA algorithm were used to quantify the abundance of infiltrating immune cell (26–28). Immune cell infiltration signatures from the study of Charoentong et al. were used to estimate the infiltration of myeloid-derived suppressor cells(MDSCs) in ssGSEA (26, 27).

Constructing the Hypoxia_DEGs_Score scoring system to evaluate individual HCC cases

To identify hypoxia-related differentially expressed genes(DEGs) in HCC, we applied the Wilcoxon rank-sum test to determine whether there was a difference in gene expression between the two statuses. The criteria FDR < 0.05, absolute fold change > 1.5, and the same direction of change among the different hypoxia statuses based on TCGA-LIHC and ICGC-LIRI-JP cohorts were used to identify the DEGs. A univariate regression model was used to analyze the hazard ratio of each DEG separately. Then, based on the survival risk of DEGs, a part was selected to be used to construct a scoring system. Based on the univariate survival regression analysis, we defined the Hypoxia_DEGs_Score of each patient through a method analogous to Gene expression grade index (GGI):

$$Hypoxia\_DEGs\_Score={\displaystyle \sum _{i=1}^n\left(bet{a}_i\times Ex{p}_i\right)},$$

where i is the selected hypoxia-related DEGs, and beta is the prognostic value of each gene signature score (29). The enrichment analysis and functional annotation of DEGs was conducted using the clusterProfiler package (30).

Analysis of tumor-associated somatic mutations and CNVs

Mutational landscapes of tumors with different hypoxia statuses were exhibited in waterfall plots using the maftools package (31). The tumor-associated mutated genes demonstrated in the waterfall plot were downloaded from the cBioPortal website (32). We performed enrichment analysis for their corresponding genes for the 40 differential CNV loci that satisfy FDR< 0.05 for differential expression, and the expression changes were consistent with the changes of CNV.

Changes of DNA methylation

The analysis of DNA methylation differences in tumors with different hypoxia status was performed using the ChAMP package (33). Before analysis, the low-quality probes were removed and the missing values were imputed. The transcriptional repression caused by DNA methylation is mainly due to CpG methylation in the promoter region, after which the binding of CPG-binding proteins in the promoter region obstructs the binding of transcription factors to the promoter (34). Therefore, we selected candidate sites based on fluctuations in the beta value of individual promoter CpG islands. Differentially methylated positions (DMPs) were selected within 2000 bases upstream of the transcription start site, based on the criteria FDR< 0.01 and |delta beta| >0.15.

Correlation of the Hypoxia_DEGs_Score with drug response

We obtained transcriptome data of 1000 cancer cell lines and the drug response data from Genomics of Drug Sensitivity in Cancer (GDSC) (35). To further analyze the association between hypoxia and drug sensitivity in HCC, we used the pRRophetic package to predict drug effects in HCC patients (36). The transcriptome data of 25 HCC cell lines from GDSC were merged with HCC transcriptome data from TCGA-LIHC after batch effects were removed using the SVA package (37), followed by principal component analysis (PCA) to assess the hypoxic status of different HCC cell lines.

Cell proliferation assay

The HCC cell lines Huh7 and HepG2 were obtained from the China Center for Type Culture Collection (Wuhan, China). Both cell lines were cultured in Dulbecco’s Modified Eagle Medium (HyClone, USA) supplemented with 10% fetal bovine serum (Gibco). HCC cell lines were seeded into 96-well plates at 4000 cells per well(4 replicates). On the 4next day, the medium in the wells was replaced after diluting the drug in equal proportions and placed into a normoxic incubator (37°C, 5% CO2, and 21% O2) or anoxic incubator (37°C, 5% CO2 and 1% O2). AKT_inhibitor_VIII, FH535, BI_2536, and RO_3306 were purchased from MedChemExpress(USA). After 72 hours, the medium was changed to 10% Cell Counting Kit 8 solution and incubated for one hour, after which the absorbance at 450 nm was recorded using an ELx 800 Universal Microplate Reader(BIO-TEK, USA).

Protein–protein interaction network of hypoxia-related DEGs

STRING was used to evaluate the associations between 251 hypoxia-related DEGs (38), and the PPI network was visualized using Cytoscape (39). The degree method in the CytoHubba package was used to distinguish the hub-genes among 251 hypoxia-related DEGs, with a threshold of 12 (40).

RNA-sequencing data from the Tongji Hospital cohort

Tumor samples were obtained from 30 HCC patients. All patients underwent one-stage radical resection in Tongji Hospital from May 2015 to November 2015, without prior anticancer treatment. This study was approved by the Research Ethics Committee of Tongji Hospital, and all patients provided written informed consent forms. RNA library construction and sequencing (Tables S8 and S9) were performed by Novogene Corporation (China).

Results

Hypoxic status in HCC stratified based on the hypoxia signature

We chose a 15-gene hypoxia signature, which exhibited the best robustness in identifying hypoxic status in various cancers in recent studies (12–14), to distinguish the hypoxic status of HCC. The unsupervised clustering of TCGA-LIHC and ICGC-LIRI-JP cohorts showed that the hypoxia status could be distinguished among different HCC samples using the 15-gene hypoxia signature (Figures 1A and S2A). Further PCA of the clustering results showed that the hypoxia feature in the three-dimensional plot also distinguishes the different hypoxia statuses of the HCC samples (Figure 1D). The high hypoxia group showed the shortest overall survival in the TCGA-LIHC and ICGC-LIRI-JP cohorts (Figures 1B, C). We found that the 15 genes exhibited significant differences in expression between the adjacent tissue and different groups of tumor samples (Figures S2B, C). The hypoxia signature network delineated a comprehensive landscape of hypoxia signature genes, hypoxia-related pathways, and prognostic factors (Figures 1E and S3A). We found a significant correlation between the expression of 15 genes in the hypoxia signature and a positive correlation between the enrichment scores (ES) of the 15 gene hypoxia signature with hallmark_hypoxia. In the TCGA-LIHC cohort, CDKN3 and TUBB6 had two missense mutations, while NDRG1 and P4HA1 had one missense mutation (Figure S3B). We explored the association between CNV and the expression of 15 signature genes, and found that ENO1, TPI1, ACOT7, NDRG1, MIF, ADM, PGAM1, P4HA1, and VEGFA were positively correlated with the high expression of the hypoxia-related genes (Figure S3C). These results indicated that the 15-gene hypoxia signature could distinguish different hypoxia statuses of HCC samples at the transcriptomic level.

FIGURE 1

Figure 1 Unsupervised clustering analysis distinguishes hypoxic phenotypes of HCC. (A) Unsupervised clustering analysis of the 15-gene hypoxia signature with a comparison of clinicopathological features between three distinct clusters for TCGA-LIHC (ICGC-LIRIJP in Figure S2A). (B, C) Kaplan-Meier curves demonstrated that hypoxia phenotypes were highly associated with the overall survival of 366 patients in the TCGA-LIHC cohort (B) and 232 patients in the ICGC-LIRI-JP cohort (C). (D) Principal component analysis for the expression profiles of three hypoxia clusters in the TCGA-LIHC cohort exhibiting a remarkable difference in expression between different hypoxia clusters. (E) The interplay between 15 hypoxia signature genes and hypoxia-related signatures in hepatocellular carcinoma. Green dots in the circle, favorable prognostic factors; yellow dots in the circle, risk factors. The lines linking the signature genes indicate their interplay, and thickness reflects the strength of correlation between the signature genes.

Distinct hypoxia clusters associated with hypoxia-related biological processes and the TME

To compare the biological behavior of tumors from high and low hypoxia clusters, we performed GSVA enrichment analysis of Hallmark and KEGG gene sets in the TCGALIHC and ICGC-LIRI-JP cohorts (Figure 2A, Table S2). Under hypoxia, we found that the glucose, fatty acid, and amino acid metabolism was upregulated; upregulated MYC targets and mTORC1 signaling pathway enhanced proliferation; downregulated PPAR promoted tumor migration (41). In the GO term-based enrichment analysis (Figure S4A, Table S3), we found that upregulation of the Rho pathway changed the cell morphology under hypoxia, thereby promoting the expansion and metastasis of HCC (42). Hypoxia promoted HCC metastasis through the WNT signaling pathway (43). Hypoxia also inhibited the differentiation of immature B cells, while also downregulating complement activation and apoptotic signaling.

FIGURE 2

Figure 2 Biological characteristics and tumor immune microenvironment of distinct hypoxia groups. (A) Heatmap visualizing the GSVA enrichment analysis based on KEGG and HALLMARK terms showing the activation states of biological pathways in hypoxia high and low groups. Red, activated pathways; blue, inhibited pathways. KEGG, Kyoto Encyclopedia of Genes and Genomes (GO in Figure S4A). (B) Enrichment score of hypoxia-related pathways across distinct hypoxia groups and normal tissue group in the TCGA-LIHC cohort. Upper panel, boxplot showing the ES distribution and overall variance. Bottom panel, heatmap showed between-group differences for pairwise comparisons. P < 0.05 in the Kruskal-Wallis (upper panel) and Wilcoxon tests (bottom panel) was considered statistically significant. ****P < 0.0001. (C) Scatter plots showing the significant association between the hypoxia-related 15-gene signature (ES) and tumor weight according to Spearman’s rank correlation analysis in the TCGA-LIHC cohort. Rs is the coefficient of rank correlation. (D) Box plots showing differences in the TIDE, Exclusion, and Disfunction scores between distinct hypoxia groups in the TCGALIHC cohort. (E, F) Box plots showing differences in the infiltration of Tregs, M1 macrophages, M2 macrophages, TANs (E), and MDSCs (F) between distinct hypoxia groups in the TCGA-LIHC cohort. The infiltration of Tregs, M1 macrophages, M2 macrophages, TANs, and MDSCs was estimated using CIBERSORT (E) and ssGSEA (F), respectively. The Kruskal-Wallis test was used to determine the statistical significance of the difference.

We applied GSVA enrichment scoring to understand the glycolysis and oxidative phosphorylation changes in HCC tumors from different hypoxia groups (Figure 2B). The change trend of hypoxia signature ES from Hallmark gene sets was consistent with that of the hypoxia clustering results. We scored ES for the 15-gene hypoxia signature and found that the results were in excellent agreement with the hypoxia cluster groupings. The ES of the hypoxia-related 15-gene set was significantly positively correlated with tumor weight (Figure 2C).

We analyzed the differences of TME between the different hypoxia statuses based on TIDE scoring (Figure 2D, Table S5). TIDE and Exclusion scores increased significantly with higher hypoxia levels. Analysis of tumor infiltration by immunosuppressive cells showed an increase of infiltration by regulatory T cells (Tregs), tumor-associated neutrophils (TANs), and MDSCs, which was correlated with aggravation of hypoxia. The increased infiltration of M1 tumor-associated macrophages(TAM) was meaningfully associated with a reduction of hypoxia (Figures 2E, F, and S4B).

Construction of a specific hypoxic signature of HCC

To further explore the changes brought about by hypoxia in HCC, we identified 251 differentially expressed signature genes (Table S4). Considering the complexity and breadth of the impact of hypoxia on tumor biology, we quantified the hypoxia signature of individual HCC patients based on a scoring model of these DEGs, called Hypoxia_DEGs_Score(Hypoxia Differential Expression Gene Score; see methods). We performed an unsupervised clustering analysis based on the expression of 251 DEGs. This analysis divided patients with different hypoxia statuses into Hypoxia_DEGs_Cluster_1 and Hypoxia_DEGs_Cluster_2 (Figure 3A). We evaluated the Hypoxia_DEGs_Score of the high, mid, and low hypoxia groups. The Hypoxia_DEGs_Score in the hypoxia group also showed an increasing trend (Figure 3B). The Hypoxia_DEGs_Score of Cluster_1 was remarkably higher than that of Cluster_2 (Figure 3C). In the Sankey diagram, the high and low grouping results of Hypoxia_DEGs_Score were highly consistent with the results of the hypoxia-related DEG clustering (Figure 3D, Table S7). And Hypoxia_DEGs_Score positively correlated with the transcript level of HIF1A (Figures S4C, D).

FIGURE 3

Figure 3 Construction of an HCC-specific hypoxia signature. (A) Unsupervised clustering of the hypoxia-related differentially expressed genes in the TCGA-LIHC and ICGC-LIRI-JP cohorts. Red, high expression; blue, low expression. (B, C) Differences in the Hypoxia_DEGs_Score between hypoxia clustered groups (B) and hypoxia DEG clusters (C) in the TCGA-LIHC cohort. (D) Sankey diagram exhibiting the shifts of 15-gene hypoxia signature clusters, hypoxia DEG clusters, and Hypoxia_DEGs_Score groups based on the TCGA-LIHC cohort. The distinction between high and low Hypoxia_DEGs_Score groups is based on the optimal cutoff value of the survival analysis. (E) Differences in the Hypoxia_DEGs_Score between various clinicopathological features. The Wilcoxon test was used to determine the statistical significance of the differences. **P < 0.01, ***P < 0.001, **** P < 0.0001.

We analyzed the 251 DEGs using the STRING database, and identified 19 hub genes with values greater than the threshold of 12 (Figure 4A). The enrichment analysis of 251 hypoxia DEGs found that the main effects of the signature of hypoxia-related DEGs were still in energy metabolism-related pathways (Figure 4B). In addition, secretion of exosomes, complement activation, and PPAR signaling were also found to be involved in the biological effects of hypoxia in HCC.

FIGURE 4

Figure 4 Validation of the Hypoxia_DEGs_Score survival prediction value. (A) Protein-protein interaction network of the hypoxia-related DEGs. The thickness of the lines represents the strength of the association, yellow circles represent the hub gene in the PPI network, and blue represents non-hub genes. DEGs with a Degree value greater than the threshold of 12 were considered hub genes. (B) Enrichment analysis for the 251 hypoxiarelated DEGs based on KEGG and GO term analysis. (C, E, G, I) Kaplan-Meier survival curves showing differences in overall survival between high and low Hypoxia_DEGs_Score groups in the TCGA-LIHC (C), ICGC-LIRI-JP (E), Zhongshan (G), and Tongji (I) cohorts. (D, F, H, K) ROC analysis of 0.5-, 1-,2- and 3-year OS demonstrating the outstanding prognostic value of Hypoxia_DEGs_Score in the TCGALIHC (D), ICGC-LIRI-JP (F), Zhongshan (H), and Tongji (K) cohorts. (J) Scatterplots of survival status and risk curves sowing the prognosis and Hypoxia_DEGs_Score values of patients in the Tongji cohort.

Clinicopathological characteristics associated with the Hypoxia_DEGs_Score

We compared the changes of Hypoxia_DEGs_Score stratified by different pathological features. A high Hypoxia_DEGs_Score was found in patients with elevated alpha-fetoprotein, higher tumor stage, worse pathological stage and grade, as well as a higher tumor recurrence rate (Figure 3E).

Initially, we analyzed the survival predicted capability of the Hypoxia_DEGs_Score model in the TCGA-LIHC and ICGC-LIRI-JP cohorts. Patients with high Hypoxia_DEGs_Score had a significantly poorer prognosis in both cohorts (Figures 4C, E, S5, B). Hypoxia_DEGs_Score also had excellent survival prediction capability in both cohorts in the ROC analysis (Figures 4D, F). Since the TCGA-LIHC and ICGCLIRI-JP cohorts were equivalent to the experimental dataset in our study, the Zhongshan cohort and our cohort from Tongji Hospital were selected as the validation datasets. Further analysis showed that the Hypoxia_DEGs_Score model also had an outstanding capability to predict the clinical outcomes in the Zhongshan and Tongji cohorts (Figures 4G–K and S5C). Therefore, the Hypoxia_DEGs_Score model exhibited excellent prognostic robustness among the results from four different cohorts from independent sources.

In multivariate Cox regression analysis, Hypoxia_DEGs_Score was also a strong independent predictor of clinical outcomes in both cohorts (Figure S5D). A nomogram was constructed for readers to predict the survival probability (Figure S5E). The 2- and 3-year calibration curves of the nomogram showed a high consistency between predicted survival probability and actual observation (Figures S6A, B).

Somatic mutations and CNVs under distinct hypoxic circumstances

Somatic mutations not only contribute to the occurrence and development of tumors, but can also change the state of various biological functions of cells (44). We analyzed the mutational profiles and found that mutation probability was lower in the low Hypoxia_DEGs_Score group (Figure 5A).

FIGURE 5

Figure 5 The somatic mutation landscape of tumor-associated genes in groups with different hypoxic statuses. (A) Waterfall plots showing the distribution of mutation among the top 20 most frequently mutated tumor-associated genes in the two Hypoxia_DEGs_Score groups from the TCGA-LIHC cohort. Each column represents an individual patient. The upper bar graph shows TMB; the number on the right indicates the mutation frequency in each tumor-associated mutated gene. The right bar graph shows the proportion of each variant type. Mutation types are indicated in the legend at the bottom. (B) The bar plot compares mutation rates for tumor-associated mutated genes in the two Hypoxia_DEGs_Score groups. (C) Forest plot showing the top 5 most significantly differentially mutated genes between the two groups. The Fisher test was used to determine the statistical significance of the differences. P values were corrected using the Benjamini-Hochberg method. *FDR < 0.05. **FDR < 0.01. ***FDR < 0.001 (D) Lollipop plots showing the amino acid residues corresponding to the mutated sites of the top 5 mutated genes. (E) Kaplan-Meier survival curves showing differences in overall survival between patients with wild-type or mutant SETD2 in the TCGA-LIHC cohort. (F) Kaplan-Meier survival curves demonstrating distinctions in overall survival stratified by both Hypoxia_DEGs_Score grouping and SETD2 mutation status in the TCGA-LIHC cohort. Mut, mutant SETD2; WT, wild-type SETD2.

Further analysis showed that the mutation frequencies of tumor suppressor genes (TP53, TSC2, CDH11, and SETD2) in the high Hypoxia_DEGs_Score group representing severe hypoxia were much higher than in the low-score group. At the same time, the mutation frequency of the proto-oncogene CTNNB1 was lower in the high-score group than in the low-score group (Figures 5B, C). Mutations of TP53 lead to the rapid proliferation of HCC and promote the progression toward hypoxia (45). Similarly, TSC2 mutations weaken mTOR regulation, leading to tumor cell proliferation and hypoxia (46). The loss of SETD2 is involved in the occurrence and development of HCC (47). The loss of CDH11 function will change the adhesion status resulting in the invasion and metastasis of HCC (48). The profiles of TP53, TSC2, CDH11, SETD2, and CTNNB1 mutation sites are shown in Figure 5D.

In further exploration, patients with SETD2 mutations had a significantly worse prognosis (Figure 5E). Notably, SETD2 gene mutations accompanied by a high hypoxia status suggested a very poor prognosis. By contrast, there was no difference in the prognosis of patients bearing tumors with mutated and wild-type SETD2 in the low hypoxia group (Figure 5F).

CNV is widely present in the human genome and plays a pivotal role in the occurrence and development of tumors (49). We further explored the association between hypoxia and CNV. Figures S6C, D shows the landscape of CNVs in the high and low Hypoxia_DEGs_Score groups. Hypoxia_DEGs_Score showed an excellent positive correlation with the frequency of CNVs of both loss and gain type (Figures 6A, B). Further analysis identified 40 differential CNV loci between different hypoxia cohorts. The CNV profiles of the top 20 and all loci are shown in Figure 6C and Table S10, respectively. We performed an enrichment analysis of possible target genes among 40 significantly different CNV loci (Figure 6D). The results suggested that the drug resistance of HCC tumors under hypoxia may be caused by a loss of genome stability, which leads to an increased frequency of CNVs, thereby altering the metabolism of cytochrome P450s and promoting drug resistance.

FIGURE 6

Figure 6 Copy number variation, DNA methylation and N6-methyl-adenosine modification in groups with different hypoxic statuses. (A, B) Scatterplots showing the robust correlation of the Hypoxia_DEGs_Score with loss (A) and gain (B) CNVs according to Spearman’s rank correlation analysis. Rs is Spearman’s coefficient of rank correlation. (C) A heatmap exhibiting the top 20 most significantly different CNV loci in distinct Hypoxia_DEGs_Score groups. The Fisher test was used to determine the statistical significance of the differences. Red, deletions; blue, gene number amplifications. (D) Enrichment analysis for the associated genes of 40 significantly different CNV loci based on KEGG gene sets. (E) Circle diagram showing the enrichment analysis for the targeted genes of DMPs based on KEGG gene sets. FDR< 0.05 for enriched pathways was considered statistically significant and was shown. Red, increased level of DNA methylation in Hypoxia_DEGs_Score high group vs. low group; Blue, decreased level of DNA methylation. (F) Boxplots showing the expression of 28 N6-methyl-adenosine regulators in groups with different hypoxic statuses. The boxes indicate the median ± 1 quartile, with the whiskers extending from the hinge to the smallest or largest value within 1.5× IQR from the box boundaries. *P < 0.05, **P < 0.01. ***P < 0.001. ****P < 0.0001.

General DNA methylation and N6-methyl-adenosine modification of HCC stratified by hypoxic status

Hypoxia can alter DNA methylation status, thereby affecting the transcription and translation of functional proteins, and altering downstream biological behaviors (50). We analyzed differentially methylated positions(DMPs) between high and low Hypoxia_DEGs_Score groups. A total of 1104 DMPs were identified, including 552 upregulated and 552 downregulated DMPs (Figure S7A, Table S11). The landscape of the top 50 DMPs is shown in Figure S7B. According to the enrichment analysis of DMP associated genes, methylation altered bile acid metabolism and possibly also inhibited the p53 signaling pathway, which significantly promoted tumor proliferation (Figure 6E). Alterations in the tight junction protein complex affect junction assembly, barrier regulation, and cell polarity, thereby promoting tumor progression and metastasis (51).

The N6-methyladenosine(m6A) DNA modification plays a crucial role in many biological functions, and its dysregulation is associated with cancer progression in HCC (52). We analyzed differences in the transcript levels of 28 m6A regulators (53) between high and low Hypoxia_DEGs_Score groups (Figure 6F). Notably, 12 of the 14 m6A readers were significantly different. The reader YTHDF1 promotes tumor progress by influencing ATG2A, ATG14, and HIF-1α (54). At the same time, there were large differences in 7 of the 11 m6A writers. High expression of WTAP promotes tumor progression, by influencing HuR, p21/27, and Ets−1 (55). Interestingly, the expression levels of the FTO and ALKBH3 m6A erasers exhibited a sharp decrease under hypoxia. In general, hypoxia promotes the occurrence and progression of HCC by increasing the abundance of the m6A DNA modification.

Potential predictive value of Hypoxia_DEGs_Score for immunotherapy response

Numerous studies have demonstrated that hypoxia can be the core driver of transformation in the tumor immune microenvironment (56). We applied the TIDE method (25) to compare the TME of tumors with different hypoxia statuses (Tables S5 and 6). TIDE and Exclusion scores were positively correlated with the Hypoxia_DEGs_Score. Notably, Dysfunction scores were not associated with hypoxia levels (Figures 7A and S8A).

FIGURE 7

Figure 7 The relationship between the Hypoxia_DEGs_Score and efficacy of immunotherapy. (A) Correlation of the Hypoxia_DEGs_Score with the TIDE score, infiltration of immunosuppressive cells, and expression of TGF-β. Rs is the coefficient of rank correlation. (B, H) Kaplan-Meier survival curves showing differences in overall survival between the high and low Hypoxia_DEGs_Score groups in the IMviogor210 (B) and GSE78220 (H) cohorts. (C, I) The proportion of patients in the IMvigor210 (C) and GSE78220 (I) cohorts with different PD-L1 and PD-1 blockade immunotherapy responses. (D, J) Difference in the Hypoxia_DEGs_Score between distinct clinical outcomes of anti-PD-L1 and PD-1 treatment in the IMvigor210 (D) and GSE78220 (J) cohorts. (E) Kaplan-Meier survival curves demonstrating distinctions in overall survival stratified by both the Hypoxia_DEGs_Score and tumor neoantigen burden. (F) Differences in PD-L1 expression between low and high Hypoxia_DEGs_Score groups in the IMviogor210 cohort. (G) Differences in the Hypoxia_DEGs_Score among immune phenotypes, including the inflamed (yellow), excluded (blue), and desert (red) immune types in the IMvigor210 cohort. CR, complete response; PR, partial response; SD, stable disease; PD, progressive disease; NEO, tumor neoantigen burden; L, low; H, high.

Further analysis of infiltration by tumor immunosuppressive cells revealed that the Hypoxia_DEGs_Score was significantly positively correlated with infiltration by Tregs, TANs, and MDSCs (Figures 7A and S8B). High infiltration levels of Tregs lead to attenuated tumor cell killing (57). TANs can promote tumor proliferation and metastasis by releasing elastase (58). The increased abundance of MDSCs in TME increases immune evasion, angiogenesis, and metastasis (59). Patients with high Hypoxia_DEGs_Score had high tumor infiltration of immunosuppressive cells, which promoted tumor progression and resulted in a poor prognosis. M1 macrophages are not only capable of secreting pro-inflammatory mediators, but can kill tumor cells by driving type 1 T helper cell(Th1) responses to produce cytotoxic effects (60). HCC samples with higher Hypoxia_DEGs_Score values had significantly less M1 macrophage infiltration (Figures 7A and S8B), suggesting an attenuated M1-mediated inflammatory state under hypoxia. TGF-β can increase the number of Tregs, suppress effector T cell activity and attenuate dendritic cell function (61). A positive correlation of the Hypoxia_DEGs_Score with the expression of TGFB1 and TGFB2 (Figure 7A and S8C) suggests that the hypoxic status of HCC results in a tumor immunosuppressive state with increased infiltration of immunosuppressive cells.

Immunotherapy based on PD-1 and PD-L1 blockade has dramatically improved the prognosis of patients. Considering that the Hypoxia_DEGs_Score can predict the TME phenotype of HCC, we selected an anti-PD-L1 treated cohort(IMvigor210) (20) and an anti-PD-1 treated cohort(GSE78220) (21) for further analysis. The results showed that patients with low Hypoxia_DEGs_Score values exhibited significant clinical benefits and markedly prolonged overall survival compared with those with a high Hypoxia_DEGs_Score (Figures 7B–D). In the IMvigor210 cohort, we found that the Hypoxia_DEGs_Score of the patients in the complete response (CR) group was lower than that of other groups, and the scores in the progressive disease(PD) group were higher than in the other groups (Figure 7D). Further analysis revealed that patients with high Hypoxia_DEGs_Score values and a low tumor neoantigen burden had the poorest prognosis (Figure 7E). PD-L1 expression was significantly elevated in patients with low Hypoxia_DEGs_Score values, suggesting a potential benefit of anti-PD-1/L1 immunotherapy (Figure 7F). In the Hypoxia_DEGs_Score analysis of different immune subtypes, we found that the immune desert type had higher scores than the immune excluded and inflammatory types (Figure 7G). Consequently, the analysis indicates that the Hypoxia_DEGs_Score model can guide ICIs drug selection and predict the efficacy of anti-PD-1/PD-L1 immunotherapy.

The Hypoxia_DEGs_Score model predicts drug responses

To further investigate the association between drug resistance and hypoxic stress, we assessed the correlation between the Hypoxia_DEGs_Score and drug responses in various cancer cell lines based on the GDSC database (35). The results revealed a significant correlation between the responses to 39 drugs and the Hypoxia_DEGs_Score (Figure S9A, Table S12). Among them, the sensitivity of the tumors to 14 drugs increased as hypoxia progressed, and conversely, tumor resistance to 25 drugs increased as hypoxia progressed. In addition, drug sensitivity in the high Hypoxia_DEGs_Score tumors was mainly associated with drugs targeting the EGFR signaling pathway. By contrast, drug tolerance was primarily associated with drugs targeting the PI3K/Akt and histone acetylation signaling pathways (Figure S9B, Table S12).

For HCC, we applied a ridge regression model to analyze the correlation of Hypoxia_DEGs_Score with the drug response based on the TCGA-LIHC and ICGC-LIRIJP datasets (Figure 8A, Table S13). Further analysis identified 20 drugs that were effective and 20 drugs to which the tumors were resistant under hypoxia. In hypoxic HCC, most effective drugs were found to target cell cycle-related signaling pathways. By contrast, hypoxic HCC appears to be resistant to most drugs targeting the PI3K/Akt signaling pathway. These results were consistent with the analysis of various cancer types.

FIGURE 8

Figure 8 Relationship between the Hypoxia_DEGs_Score and drug response. (A) Correlation between the Hypoxia_DEGs_Score and drug sensitivity assessed using Spearman correlation analysis in the TCGA-LIHC and ICGCLIRI-JP cohort. The upper left triangle shows the analysis of the TCGA-LIHC cohort; the lower right triangle shows the analysis of the ICGC-LIRI-JP cohort. Red indicates ineffective drugs (tumor is resistant), and blue indicates effective drugs (tumor is sensitive). Rs is the coefficient of rank correlation. Each small rectangle on the right represent the targeted pathways corresponding to the drugs. *P < 1.0e-04; •P < 1.0e-12 (B) Principal component analysis for the transcriptome profiles of 251 hypoxia-related DEGs based on the combination between HCC cell lines and patient samples from the GDSC dataset and the TCGA-LIHC cohort. Red represents HCC samples and cell lines with high Hypoxia_DEGs_Score values; blue Represents low Hypoxia_DEGs_Score. (C) Dose-response curves for the the mean value of cell viability for BI-2536, RO-3306, FH535, and Akt inhibitor-VIII in the HCC cell lines Huh7 (C) and HepG2 (Figure S12C) incubated under hypoxic (red, n = 4) and normoxic (blue, n = 4) conditions. Cell viability was normalized to that of cells mock-treated with dimethyl sulfoxide (vehicle control). The error bars indicate the mean ± SD. The sigmoidal dose-response curves were fitted to the data.

To directly confirm the association between the Hypoxia_DEGs_Score and the drug response, we conducted a PCA analysis. The results showed that HepG2, Huh7, and Huh1 were the HCC cell lines with the lowest hypoxia levels, while HLF, JHH2, and SNU387 had the highest hypoxia levels (Figure 8B). Additionally, we cultured HepG2 and Huh7 cells under hypoxic and normoxic conditions, and treated them with representative drugs. The results indicated that under hypoxic conditions, HCC exhibited higher drug sensitivity to BI-2536 and RO-3306 but higher resistance to FH535 and AKT inhibitor VIII (Figures 8C and S9C). These results were consistent with our computational prediction. Thus, the Hypoxia_DEGs_Score is a potential biomarker for establishing appropriate therapeutic regimens.

Discussion

Due to differences in organ function, metabolic background, and blood supply, the hypoxia status varies significantly in different tumors. As a heterogeneous disease, HCC is subject to comprehensive and diverse biological processes, resulting in a variable and inconsistent prognosis. The existing biomarkers used to guide diagnosis and treatment cannot fully meet the clinical needs. Therefore, establishing new biomarkers with higher predictive value is critical for maximizing the clinical benefits of a personalized regimen and accurate prognostic assessment. We distinguished three HCC hypoxia phenotypes based on 15 hypoxia signature genes and constructed the HCC-specific hypoxia scoring model Hypoxia_DEGs_Score. The robustness of the Hypoxia_DEGs_Score model was validated in four independent HCC cohorts from different sources (Figures 4C–K, S5A–C, and Table S1). Encouragingly, the Zhongshan and Tongji cohorts used as validation datasets also had a prominent area under curve (AUC) of the ROC, and patients with low Hypoxia_DEGs_Score values had a significantly better prognosis.

Hypoxia is a powerful driving force that promotes the heterogenization and genomic evolution of cancer. However, the understanding its association with specific mutational processes and CNVs remains limited (62). In our analysis, the mutation frequency of tumor driver genes generally increased under hypoxia (Figure 5A), and the frequency of CNVs also surged with rising Hypoxia_DEGs_Score values (Figures 6A, B). These results are consistent with the widespread elevation of the mutational burden across all mutational types related to hypoxia in pan-cancer studies (62). TP53 mutations prevent apoptosis and lead tumors into the dynamic cycle of hypoxia (63, 64). The high mutation rate of TP53 in the high hypoxia group indicated that TP53 mutations play a pivotal role in the development of HCC under hypoxia (Figures 5B, C). It has been reported that SETD2 mutations were more frequently found in the hypoxic group of kidney tumors (65, 66), and the hypoxic group of HCC also had more patients with mutant SETD2 in our study (Figures 5B, C). It is worth noting that patients with both mutant SETD2 and a hyper-hypoxic status had the worst prognoses (Figures 5E, F). An enriched drug metabolism-cytochrome P450 gene set based on the target genes of hypoxic differential CNV loci indicated that CNVs of cytochrome P450 genes may enhance drug resistance in hypoxic HCC, in agreement with previous literature (67). In our study, the PPAR signaling pathway was enriched in three analyses, respectively in the GSVA of high and low clustering groups, 251 DEG enrichment analysis, and CNV loci-related genes enrichment analysis (Figures 2A, 4B, and 6D). Downregulation of the PPAR signaling pathway has been reported to promote tumor growth and the EMT resulting in metastasis. Therefore, CNVs may promote the tumor growth and metastasis of HCC by downregulating the PPAR signaling pathway under hypoxia.

Aberrant methylation of genes affects multiple biological functions of tumors, such as cell cycle, DNA repair, toxin catabolism, cell adherence, apoptosis, and angiogenesis (68). In our study, hypoxia inhibited or promoted the expression of various genes through methylation modification of DNA, including genes related to cholesterol metabolism, p53 signaling, tight junction protein function and ATP-binding cassette (ABC) transporters, thereby altering tumor metabolism, cell cycle, DNA repair, cell adherence, and drug catabolism (Figure 6E). Previous reports indicate that increased expression of ABCA1 and ABCC1 is positively correlated with chemotherapy resistance (69–72). In our analysis, ABCA1 and ABCC1 were methylated under hypoxia, suggesting that hypoxia increases the number of ABC transporters and promotes drug resistance in HCC by downregulating the methylation of ABCA1 and ABCC1.

Hypoxia is an essential feature of the TME and is closely linked to cell proliferation, angiogenesis, metabolism, and tumor immunity (73). Our study showed that a high Hypoxia_DEGs_Score was positively correlated with the TIDE score (Figures 7A and S8A). A high TIDE score indicates that the patient’s tumor is immunosuppressed, which predicts a poor response to ICIs treatment (25). The positive correlation between the Hypoxia_DEGs_Score and the Exclusion score indicated that hypoxic immunosuppression in HCC was dominated by immune cell reduction rather than immune dysfunction (Figures 7A and S8A). The increased infiltration of TAN, MDSC, and Treg cells indicates the suppression of immune function and promotion of tumor progression in HCC under hypoxia (Figures 7A and S8B). Furthermore, TGFB1 and TGFB2 were positively correlated with the Hypoxia_DEGs_Score, further illustrating the immunosuppressive state induced by hypoxia. We evaluated the predictive power of Hypoxia_DEGs_Score for responses to anti-PD-L1 and anti-PD-1 immunotherapy. Although the data for these assessments were not derive from HCC cohorts, the results showed that predicting survival time and treatment response was significantly associated with the Hypoxia_DEGs_Score (Figures 7B–J). Consistently, patients with high scores had a poor prognosis in four independent cohorts (Figures 4C–K, S5A–C). At the same time, the immunotherapy outcome was also compatible with the TIDE score and TME status of HCC under aggravated hypoxia (Figures 7A and S8). These results indicated that patients with high Hypoxia_DEGs_Score values have a poorer prognosis and may not benefit from ICIs therapy.

Hypoxia promotes drug resistance through apoptosis, autophagy, DNA damage, mitochondrial activity, p53, drug efflux, and an acidic microenvironment (56). Therefore, it is vital to find a regime that can effectively eliminate tumors in a hypoxic environment. The pan-cancer analysis was not limited to HCC, so we analyzed the drug response of HCC under hypoxic stress based on TCGA-LIHC and ICGC-LIRI-JP. The primary targeted pathway of effective drugs was cell cycle-related, and the main targeted pathway of drugs to which the tumors were resistant was PI3K/Akt signalling (Figure 8A). To confirm these findings, we placed HCC cell lines in hypoxic and normoxic environments to assess their drug sensitivity (Figures 8C and S9C), and the results were consistent with the predictions. Notably, the response to AKT inhibitor VIII and RO-3306, which respectively target the PI3K/mTOR and cell cycle pathways, was consistent with the pan-cancer drug prediction results, but FH535 and BI-2536 did not appear in the results. These analyses suggest that patients with higher Hypoxia_DEGs_Score values may benefit from drugs against cell cycle signaling pathways rather than those inhibiting the PI3K/mTOR pathway. It is suggested that Hypoxia_DEGs_Score could serve as a potential biomarker to guide drug selection for patients.

The advancement of sequencing technology and understanding of the tumor microenvironment, as well as the combination of molecular and genetic features, and clinicopathological factors, offer new hope for precise prediction and personalized therapy. Therefore, genetic signatures bridging hypoxia and prognosis offer promising guidelines for exploring and treating HCC. Our study provides new possibilities for improving therapy outcomes in hepatocellular carcinoma by identifying distinguishing hypoxic states and enabling personalized regimes.

In conclusion, we performed a systematic, comprehensive analysis of the effects of hypoxia on multiple HCC cohorts based on the hypoxia signature. The Hypoxia_DEGs_Score prediction model was constructed, and exhibited an excellent survival prediction value. It was revealed that hypoxia not only has a complex effect on the TME, but also has an extensive impact on genome instability and DNA methylation. A comprehensive evaluation of the hypoxia status of individual HCC cases will help us deepen our understanding of the altered biological functions of HCC and reveal the key role of hypoxia in targeted therapy and immunotherapy. Hypoxia_DEGs_Score can greatly aid the development of optimal personalized immunotherapy regimens for HCC patients, with the potential to greatly improve the prognosis.

Data availability statement

The Tongji cohort's RNA sequencing data presented in the study are deposited in the Gene Expression Omnibus repository (https://www.ncbi.nlm.nih.gov/geo/), accession number GSE201868. The other RNA sequencing, Whole Genome Sequencing, methylation chip, and clinicopathological data of HCC patients, or patients with immune checkpoint were described in method section “Data collection and processing” and Supplementary Material. The resources, tools and codes used in our analyses were described in each method section in the methods.

Author contributions

GC, JZ and DN designed this study. YZ, YX, JL, and GL extracted the information from the databases and analyzed the data. GC completed the verification experiments. CY, XC, HL, and ZD supervised the entire study. GC wrote the manuscript. All authors revised the manuscript.

Funding

This work was supported by the National Nature Science Foundation of China (No’s. 81874065 and 81874149), National Basic Research Program of China (2020YFA0710700), Tongji Hospital (HUST) Fundation for Excellent Young Scientist (No. 2020YQ05), the first level of the public health youth top talent project of Hubei province.

Acknowledgments

The authors would like to sincerely thank the TCGA, ICGC, Zhongshan hosipital, GEO, IMvigor, and GDSC for data sharing, as well as TIDE, cibersort, and, cBioPortal for providing data processing. In particular, the author Xiaoping Chen would like to express deepest thanks to Dr. Ivan Hajal for his support in revising this manuscript.

Conflict of interest

ZD served as a speaker and consultant for Bayer, Eisai, Roche, MSD, Astra-Zeneca, Innovent, Hengrui, and BeiGene. XC served as a consultant of Bayer, Eisai, and Hengrui. XC conducts studies for Hengrui and CTTQ.

The remaining authors declare that the research was conducted in the absence of any commercial or financial relationships that could be constructed as a potential conflict of interest.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fimmu.2022.997316/full#supplementary-material

Supplementary Figure 1 | Overview of study design. Flowchart of the steps in the performed analyses.

Supplementary Figure 2 | Unsupervised clustering and expression of 15 hypoxia signature genes. (A) Unsupervised clustering analysis of the 15 genes hypoxia signature with the comparison of clinicopathological features between three distinct clusters for ICGC-LIRI-JP. (B, C) The expression of 15 hypoxia signature genes across distinct hypoxia groups and normal tissue group in TCGA-LIHC (B) and ICGC-LIRI-JP (C) cohorts. *, P< 0.05; **, P< 0.01; ***, P< 0.001; ****, P< 0.0001.

Supplementary Figure 3 | Correlation analysis and mutation status of 15 hypoxia signature genes. (A) Heatmap displays the correlation between 15 hypoxia signature genes and hypoxia-related signatures. Exhibited values represent correlation coefficients and were satisfied with FDR< 0.05. (B) The mutation frequency of 15 hypoxia signature genes in 364 patients with HCC from TCGA-LIHC cohort. Each column represented individual patients. (C) The expression of 15 hypoxia signature genes between normal, gain, and loss of CNV. *, P < 0.05; **, P < 0.01; ***, P < 0.001; ****, P < 0.0001.

Supplementary Figure 4 | Biological characteristics and immune cell infiltration of distinct hypoxia groups. (A) Heatmap visualizing the GSVA enrichment analysis based on GO terms shows the activation states of biological pathways in hypoxia high and low groups. (B) The infiltration of 22 immune cells between distinct hypoxia groups in TCGA-LIHC cohorts. * P < 0.05; ** P < 0.01; *** P < 0.001; **** P < 0.0001. Correlation analysis of Hypoxia_DEGs_Score with HIF1A transcript levels in TCGA (C) and ICGC cohorts (D).

Supplementary Figure 5 | Survival status in the different hypoxic statuses and multivariate Cox regression analysis and comprehensive nomogram for the Hypoxia_DEGs_Score. (A–C) Scatter plots of survival status and risk curves exhibited prognosis and Hypoxia_DEGs_Score in patients from TCGA-LIHC (A), ICGC-LIRI-JP (B), and Zhongshan (C) cohorts. (D) Forest plots showing multivariate Cox regression model analysis in the TCGA-LIHC and ICGC-LIRI-JP cohorts. The horizontal line represents the 95% confidence interval (CI) for each group. The vertical dotted line represents all patients’ hazard ratios (HR). (E) The comprehensive nomogram predicting the clinical outcomes of HCC patients with 2- and 3-year OS based on the TCGA-LIHC cohort.

Supplementary Figure 6 | Calibration of comprehensive nomogram and landscape of CNV of different hypoxic statuses. (A, B) Calibration plots for predicting the 2-year (A) and 3-year (B) OS of HCC patients in the TCGA-LIHC cohort. Nomogram-predicted probability of survival was plotted on the x-axis, and actual survival was plotted on the y-axis. (C, D) The landscape of CNV in high (A) and low (B) Hypoxia_DEGs_Score groups.

Supplementary Figure 7 | The beta value of DNA methylation differential sites in the different hypoxic statuses. (A) Heatmap exhibiting the beta value of 1104 DNA methylation differential sites in high and low Hyposia_DEGs_Score groups. (B) Heatmap displayed the beta value of top 50 DNA methylation differential sites in high and low Hyposia_DEGs_Score groups.

Supplementary Figure 8 | Changes of the tumor immune microenvironment in groups with different hypoxic statuses. (A–C) Scatterplots showing the correlation of the Hypoxia_DEGs_Score with the TIDE score (A), infiltration of immunosuppressive cells (B), and expression of TGF-β (C) according to Spearman’s rank correlation analysis in the TCGA-LIHC and ICGCLIRI-JP cohorts. Rs is the coefficient of rank correlation.

Supplementary Figure 9 | Relationship between the Hypoxia_DEGs_Score and drug response. (A) Correlation between the Hypoxia_DEGs_Score and drug sensitivity based on Spearman correlation analysis in GDSC. Each column represents a drug. The brightness of the column indicates the significance of the correlation. Rs is the coefficient of rank correlation. The height of the column indicates the correlation between Hypoxia_DEGs_Score and drug resistance (Rs > 0) or drug sensitivity (Rs< 0). (B) Signaling pathways targeted by ineffective (red) or effective (blue) drugs according to the Hypoxia_DEGs_Score. Drug names were listed on the horizontal axis and the signaling pathway targeted by the drug on the vertical axis. The bar graph on the right shows the number of drugs targeting each signaling pathway. (C) Dose-response curves for the the mean value of cell viability for BI-2536, RO-3306, FH535, and Akt inhibitor-VIII in the HCC cell lines HepG2 under hypoxic (red, n = 4) and normoxic (blue, n = 4) conditions.

References

1. Bray F, Ferlay J, Soerjomataram I, Siegel RL, Torre LA, Jemal A. Global cancer statistics 2018: Globocan estimates of incidence and mortality worldwide for 36 cancers in 185 countries. CA Cancer J Clin (2018) 68(6):394–424. doi: 10.3322/caac.21492

PubMed Abstract | CrossRef Full Text | Google Scholar

2. Al-Salama ZT, Syed YY, Scott LJ. Lenvatinib: A review in hepatocellular carcinoma. Drugs (2019) 79(6):665–74. doi: 10.1007/s40265-019-01116-x

PubMed Abstract | CrossRef Full Text | Google Scholar

3. Li H, Li CW, Li X, Ding Q, Guo L, Liu S, et al. Met inhibitors promote liver tumor evasion of the immune response by stabilizing Pdl1. Gastroenterology (2019) 156(6):1849–61 e13. doi: 10.1053/j.gastro.2019.01.252

PubMed Abstract | CrossRef Full Text | Google Scholar

4. Craig AJ, von Felden J, Garcia-Lezana T, Sarcognato S, Villanueva A. Tumour evolution in hepatocellular carcinoma. Nat Rev Gastroenterol Hepatol (2020) 17(3):139–52. doi: 10.1038/s41575-019-0229-4

PubMed Abstract | CrossRef Full Text | Google Scholar

5. Graham K, Unger E. Overcoming tumor hypoxia as a barrier to radiotherapy, chemotherapy and immunotherapy in cancer treatment. Int J Nanomed (2018) 13:6049–58. doi: 10.2147/IJN.S140462

CrossRef Full Text | Google Scholar

6. Hockel M, Vaupel P. Tumor hypoxia: Definitions and current clinical, biologic, and molecular aspects. J Natl Cancer Inst (2001) 93(4):266–76. doi: 10.1093/jnci/93.4.266

PubMed Abstract | CrossRef Full Text | Google Scholar

7. Wigerup C, Pahlman S, Bexell D. Therapeutic targeting of hypoxia and hypoxia-inducible factors in cancer. Pharmacol Ther (2016) 164:152–69. doi: 10.1016/j.pharmthera.2016.04.009

PubMed Abstract | CrossRef Full Text | Google Scholar

8. Bertout JA, Patel SA, Simon MC. The impact of O2 availability on human cancer. Nat Rev Cancer (2008) 8(12):967–75. doi: 10.1038/nrc2540

PubMed Abstract | CrossRef Full Text | Google Scholar

9. Hiraga T. Hypoxic microenvironment and metastatic bone disease. Int J Mol Sci (2018) 19(11):3523. doi: 10.3390/ijms19113523

CrossRef Full Text | Google Scholar

10. Liu Z, Wang Y, Dou C, Xu M, Sun L, Wang L, et al. Hypoxia-induced up-regulation of vasp promotes invasiveness and metastasis of hepatocellular carcinoma. Theranostics (2018) 8(17):4649–63. doi: 10.7150/thno.26789

PubMed Abstract | CrossRef Full Text | Google Scholar

11. Mendez-Blanco C, Fondevila F, Garcia-Palomo A, Gonzalez-Gallego J, Mauriz JL. Sorafenib resistance in hepatocarcinoma: Role of hypoxia-inducible factors. Exp Mol Med (2018) 50(10):1–9. doi: 10.1038/s12276-018-0159-1

PubMed Abstract | CrossRef Full Text | Google Scholar

12. Buffa FM, Harris AL, West CM, Miller CJ. Large Meta-analysis of multiple cancers reveals a common, compact and highly prognostic hypoxia metagene. Br J Cancer (2010) 102(2):428–35. doi: 10.1038/sj.bjc.6605450

PubMed Abstract | CrossRef Full Text | Google Scholar

13. Fox NS, Starmans MH, Haider S, Lambin P, Boutros PC. Ensemble analyses improve signatures of tumour hypoxia and reveal inter-platform differences. BMC Bioinf (2014) 15:170. doi: 10.1186/1471-2105-15-170

CrossRef Full Text | Google Scholar

14. Thienpont B, Steinbacher J, Zhao H, D'Anna F, Kuchnio A, Ploumakis A, et al. Tumour hypoxia causes DNA hypermethylation by reducing tet activity. Nature (2016) 537(7618):63–8. doi: 10.1038/nature19081

PubMed Abstract | CrossRef Full Text | Google Scholar

15. Ye Y, Hu Q, Chen H, Liang K, Yuan Y, Xiang Y, et al. Characterization of hypoxia-associated molecular features to aid hypoxia-targeted therapy. Nat Metab (2019) 1(4):431–44. doi: 10.1038/s42255-019-0045-8

PubMed Abstract | CrossRef Full Text | Google Scholar

16. Walsh JC, Lebedev A, Aten E, Madsen K, Marciano L, Kolb HC. The clinical importance of assessing tumor hypoxia: Relationship of tumor hypoxia to prognosis and therapeutic opportunities. Antioxid Redox Signal (2014) 21(10):1516–54. doi: 10.1089/ars.2013.5378

PubMed Abstract | CrossRef Full Text | Google Scholar

17. Tomczak K, Czerwinska P, Wiznerowicz M. The cancer genome atlas (Tcga): An immeasurable source of knowledge. Contemp Oncol (Pozn) (2015) 19(1A):A68–77. doi: 10.5114/wo.2014.47136

PubMed Abstract | CrossRef Full Text | Google Scholar

18. Fujimoto A, Furuta M, Totoki Y, Tsunoda T, Kato M, Shiraishi Y, et al. Whole-genome mutational landscape and characterization of noncoding and structural mutations in liver cancer. Nat Genet (2016) 48(5):500–9. doi: 10.1038/ng.3547

PubMed Abstract | CrossRef Full Text | Google Scholar

19. Gao Q, Zhu H, Dong L, Shi W, Chen R, Song Z, et al. Integrated proteogenomic characterization of hbv-related hepatocellular carcinoma. Cell (2019) 179(2):561–77 e22. doi: 10.1016/j.cell.2019.08.052

PubMed Abstract | CrossRef Full Text | Google Scholar

20. Mariathasan S, Turley SJ, Nickles D, Castiglioni A, Yuen K, Wang Y, et al. Tgfbeta attenuates tumour response to pd-L1 blockade by contributing to exclusion of T cells. Nature (2018) 554(7693):544–8. doi: 10.1038/nature25501

PubMed Abstract | CrossRef Full Text | Google Scholar

21. Hugo W, Zaretsky JM, Sun L, Song C, Moreno BH, Hu-Lieskovan S, et al. Genomic and transcriptomic features of response to anti-Pd-1 therapy in metastatic melanoma. Cell (2016) 165(1):35–44. doi: 10.1016/j.cell.2016.02.065

PubMed Abstract | CrossRef Full Text | Google Scholar

22. Hartigan JA, Wong MA. Algorithm as 136: A K-means clustering algorithm. J R Stat Soc (1979) 28(1):8. doi: 10.2307/2346830

CrossRef Full Text | Google Scholar

23. Wilkerson MD, Hayes DN. Consensusclusterplus: A class discovery tool with confidence assessments and item tracking. Bioinf (Oxford England) (2010) 26(12):1572–3. doi: 10.1093/bioinformatics/btq170

CrossRef Full Text | Google Scholar

24. Hanzelmann S, Castelo R, Guinney J. Gsva: Gene set variation analysis for microarray and rna-seq data. BMC Bioinf (2013) 14:7. doi: 10.1186/1471-2105-14-7

CrossRef Full Text | Google Scholar

25. Jiang P, Gu S, Pan D, Fu J, Sahu A, Hu X, et al. Signatures of T cell dysfunction and exclusion predict cancer immunotherapy response. Nat Med (2018) 24(10):1550–8. doi: 10.1038/s41591-018-0136-1

PubMed Abstract | CrossRef Full Text | Google Scholar

26. Charoentong P, Finotello F, Angelova M, Mayer C, Efremova M, Rieder D, et al. Pan-cancer immunogenomic analyses reveal genotype-immunophenotype relationships and predictors of response to checkpoint blockade. Cell Rep (2017) 18(1):248–62. doi: 10.1016/j.celrep.2016.12.019

PubMed Abstract | CrossRef Full Text | Google Scholar

27. Jia Q, Wu W, Wang Y, Alexander PB, Sun C, Gong Z, et al. Local mutational diversity drives intratumoral immune heterogeneity in non-small cell lung cancer. Nat Commun (2018) 9(1):5361. doi: 10.1038/s41467-018-07767-w

PubMed Abstract | CrossRef Full Text | Google Scholar

28. Newman AM, Liu CL, Green MR, Gentles AJ, Feng W, Xu Y, et al. Robust enumeration of cell subsets from tissue expression profiles. Nat Methods (2015) 12(5):453–7. doi: 10.1038/nmeth.3337

PubMed Abstract | CrossRef Full Text | Google Scholar

29. Sotiriou C, Wirapati P, Loi S, Harris A, Fox S, Smeds J, et al. Gene expression profiling in breast cancer: Understanding the molecular basis of histologic grade to improve prognosis. J Natl Cancer Inst (2006) 98(4):262–72. doi: 10.1093/jnci/djj052

PubMed Abstract | CrossRef Full Text | Google Scholar

30. Yu G, Wang LG, Han Y, He QY. Clusterprofiler: An r package for comparing biological themes among gene clusters. OMICS (2012) 16(5):284–7. doi: 10.1089/omi.2011.0118

PubMed Abstract | CrossRef Full Text | Google Scholar

31. Mayakonda A, Lin DC, Assenov Y, Plass C, Koeffler HP. Maftools: Efficient and comprehensive analysis of somatic variants in cancer. Genome Res (2018) 28(11):1747–56. doi: 10.1101/gr.239244.118

PubMed Abstract | CrossRef Full Text | Google Scholar

32. Gao J, Aksoy BA, Dogrusoz U, Dresdner G, Gross B, Sumer SO, et al. Integrative analysis of complex cancer genomics and clinical profiles using the cbioportal. Sci Signal (2013) 6(269):l1. doi: 10.1126/scisignal.2004088

CrossRef Full Text | Google Scholar

33. Tian Y, Morris TJ, Webster AP, Yang Z, Beck S, Feber A, et al. Champ: Updated methylation analysis pipeline for illumina beadchips. Bioinf (Oxford England) (2017) 33(24):3982–4. doi: 10.1093/bioinformatics/btx513

CrossRef Full Text | Google Scholar

34. Greenberg MVC, Bourc'his D. The diverse roles of DNA methylation in mammalian development and disease. Nat Rev Mol Cell Biol (2019) 20(10):590–607. doi: 10.1038/s41580-019-0159-6

PubMed Abstract | CrossRef Full Text | Google Scholar

35. Yang W, Soares J, Greninger P, Edelman EJ, Lightfoot H, Forbes S, et al. Genomics of drug sensitivity in cancer (Gdsc): A resource for therapeutic biomarker discovery in cancer cells. Nucleic Acids Res (2013) 41(Database issue):D955–61. doi: 10.1093/nar/gks1111

PubMed Abstract | CrossRef Full Text | Google Scholar

36. Geeleher P, Cox N, Huang RS. Prrophetic: An r package for prediction of clinical chemotherapeutic response from tumor gene expression levels. PloS One (2014) 9(9):e107468. doi: 10.1371/journal.pone.0107468

PubMed Abstract | CrossRef Full Text | Google Scholar

37. Leek JT, Johnson WE, Parker HS, Jaffe AE, Storey JD. The sva package for removing batch effects and other unwanted variation in high-throughput experiments. Bioinf (Oxford England) (2012) 28(6):882–3. doi: 10.1093/bioinformatics/bts034

CrossRef Full Text | Google Scholar

38. Szklarczyk D, Franceschini A, Wyder S, Forslund K, Heller D, Huerta-Cepas J, et al. String V10: Protein-protein interaction networks, integrated over the tree of life. Nucleic Acids Res (2015) 43(Database issue):D447–52. doi: 10.1093/nar/gku1003

PubMed Abstract | CrossRef Full Text | Google Scholar

39. Shannon P, Markiel A, Ozier O, Baliga NS, Wang JT, Ramage D, et al. Cytoscape: A software environment for integrated models of biomolecular interaction networks. Genome Res (2003) 13(11):2498–504. doi: 10.1101/gr.1239303

PubMed Abstract | CrossRef Full Text | Google Scholar

40. Chin CH, Chen SH, Wu HH, Ho CW, Ko MT, Lin CY. Cytohubba: Identifying hub objects and Sub-networks from complex interactome. BMC Syst Biol (2014) 8 Suppl 4:S11. doi: 10.1186/1752-0509-8-S4-S11

PubMed Abstract | CrossRef Full Text | Google Scholar

41. Silva VR, Santos LS, Dias RB, Quadros CA, Bezerra DP. Emerging agents that target signaling pathways to eradicate colorectal cancer stem cells. Cancer Commun (Lond) (2021) 41(12):1275–313. doi: 10.1002/cac2.12235

PubMed Abstract | CrossRef Full Text | Google Scholar

42. Wong CC, Wong CM, Au SL, Ng IO. Rhogtpases and rho-effectors in hepatocellular carcinoma metastasis: Rock n'rho move it. Liver Int (2010) 30(5):642–56. doi: 10.1111/j.1478-3231.2010.02232.x

PubMed Abstract | CrossRef Full Text | Google Scholar

43. Xu C, Xu Z, Zhang Y, Evert M, Calvisi DF, Chen X. Beta-catenin signaling in hepatocellular carcinoma. J Clin Invest (2022) 132(4):e154515. doi: 10.1172/JCI154515

PubMed Abstract | CrossRef Full Text | Google Scholar

44. Garcia-Lezana T, Lopez-Canovas JL, Villanueva A. Signaling pathways in hepatocellular carcinoma. Adv Cancer Res (2021) 149:63–101. doi: 10.1016/bs.acr.2020.10.002

PubMed Abstract | CrossRef Full Text | Google Scholar

45. Calderaro J, Ziol M, Paradis V, Zucman-Rossi J. Molecular and histological correlations in liver cancer. J Hepatol (2019) 71(3):616–30. doi: 10.1016/j.jhep.2019.06.001

PubMed Abstract | CrossRef Full Text | Google Scholar

46. Zucman-Rossi J, Villanueva A, Nault JC, Llovet JM. Genetic landscape and biomarkers of hepatocellular carcinoma. Gastroenterology (2015) 149(5):1226–39 e4. doi: 10.1053/j.gastro.2015.05.061

PubMed Abstract | CrossRef Full Text | Google Scholar

47. Li XJ, Li QL, Ju LG, Zhao C, Zhao LS, Du JW, et al. Deficiency of histone methyltransferase set domain-containing 2 in liver leads to abnormal lipid metabolism and hcc. Hepatol (Baltimore Md) (2021) 73(5):1797–815. doi: 10.1002/hep.31594

CrossRef Full Text | Google Scholar

48. Chen X, Xiang H, Yu S, Lu Y, Wu T. Research progress in the role and mechanism of cadherin-11 in different diseases. J Cancer (2021) 12(4):1190–9. doi: 10.7150/jca.52720

PubMed Abstract | CrossRef Full Text | Google Scholar

49. Lauer S, Gresham D. An evolving view of copy number variants. Curr Genet (2019) 65(6):1287–95. doi: 10.1007/s00294-019-00980-0

PubMed Abstract | CrossRef Full Text | Google Scholar

50. Choudhry H, Harris AL. Advances in hypoxia-inducible factor biology. Cell Metab (2018) 27(2):281–98. doi: 10.1016/j.cmet.2017.10.005

PubMed Abstract | CrossRef Full Text | Google Scholar

51. Zeisel MB, Dhawan P, Baumert TF. Tight junction proteins in gastrointestinal and liver disease. Gut (2019) 68(3):547–61. doi: 10.1136/gutjnl-2018-316906

PubMed Abstract | CrossRef Full Text | Google Scholar

52. Sivasudhan E, Blake N, Lu ZL, Meng J, Rong R. Dynamics of M6a rna methylome on the hallmarks of hepatocellular carcinoma. Front Cell Dev Biol (2021) 9:642443. doi: 10.3389/fcell.2021.642443

PubMed Abstract | CrossRef Full Text | Google Scholar

53. Deng LJ, Deng WQ, Fan SR, Chen MF, Qi M, Lyu WY, et al. M6a modification: Recent advances, anticancer targeted drug discovery and beyond. Mol Cancer (2022) 21(1):52. doi: 10.1186/s12943-022-01510-2

PubMed Abstract | CrossRef Full Text | Google Scholar

54. Li Q, Ni Y, Zhang L, Jiang R, Xu J, Yang H, et al. Hif-1alpha-Induced expression of M6a reader Ythdf1 drives hypoxia-induced autophagy and malignancy of hepatocellular carcinoma by promoting Atg2a and Atg14 translation. Signal transduct Targeted Ther (2021) 6(1):76. doi: 10.1038/s41392-020-00453-8

CrossRef Full Text | Google Scholar

55. Chen Y, Peng C, Chen J, Chen D, Yang B, He B, et al. Wtap facilitates progression of hepatocellular carcinoma via M6a-Hur-Dependent epigenetic silencing of Ets1. Mol Cancer (2019) 18(1):127. doi: 10.1186/s12943-019-1053-8

PubMed Abstract | CrossRef Full Text | Google Scholar

56. Jing X, Yang F, Shao C, Wei K, Xie M, Shen H, et al. Role of hypoxia in cancer therapy by regulating the tumor microenvironment. Mol Cancer (2019) 18(1):157. doi: 10.1186/s12943-019-1089-9

PubMed Abstract | CrossRef Full Text | Google Scholar

57. Nishikawa H, Koyama S. Mechanisms of regulatory T cell infiltration in tumors: Implications for innovative immune precision therapies. J Immunother Cancer (2021) 9(7):e002591. doi: 10.1136/jitc-2021-002591

PubMed Abstract | CrossRef Full Text | Google Scholar

58. Wu L, Zhang XH. Tumor-associated neutrophils and macrophages-heterogenous but not chaotic. Front Immunol (2020) 11:553967. doi: 10.3389/fimmu.2020.553967

PubMed Abstract | CrossRef Full Text | Google Scholar

59. Mortezaee K. Myeloid-derived suppressor cells in cancer immunotherapy-clinical perspectives. Life Sci (2021) 277:119627. doi: 10.1016/j.lfs.2021.119627

PubMed Abstract | CrossRef Full Text | Google Scholar

60. Mills CD. M1 and M2 macrophages: Oracles of health and disease. Crit Rev Immunol (2012) 32(6):463–88. doi: 10.1615/critrevimmunol.v32.i6.10

PubMed Abstract | CrossRef Full Text | Google Scholar

61. Lind H, Gameiro SR, Jochems C, Donahue RN, Strauss J, Gulley JM, et al. Dual targeting of tgf-beta and pd-L1 via a bifunctional anti-Pd-L1/Tgf-Betarii agent: Status of preclinical and clinical advances. J Immunother Cancer (2020) 8(1):e000433. doi: 10.1136/jitc-2019-000433

PubMed Abstract | CrossRef Full Text | Google Scholar

62. Bhandari V, Li CH, Bristow RG, Boutros PC, Consortium P. Divergent mutational processes distinguish hypoxic and normoxic tumours. Nat Commun (2020) 11(1):737. doi: 10.1038/s41467-019-14052-x

PubMed Abstract | CrossRef Full Text | Google Scholar

63. Greijer AE, van der Wall E. The role of hypoxia inducible factor 1 (Hif-1) in hypoxia induced apoptosis. J Clin Pathol (2004) 57(10):1009–14. doi: 10.1136/jcp.2003.015032

PubMed Abstract | CrossRef Full Text | Google Scholar

64. Graeber TG, Osmanian C, Jacks T, Housman DE, Koch CJ, Lowe SW, et al. Hypoxia-mediated selection of cells with diminished apoptotic potential in solid tumours. Nature (1996) 379(6560):88–91. doi: 10.1038/379088a0

PubMed Abstract | CrossRef Full Text | Google Scholar

65. Kovac M, Navas C, Horswell S, Salm M, Bardella C, Rowan A, et al. Recurrent chromosomal gains and heterogeneous driver mutations characterise papillary renal cancer evolution. Nat Commun (2015) 6:6336. doi: 10.1038/ncomms7336

PubMed Abstract | CrossRef Full Text | Google Scholar

66. Dalgliesh GL, Furge K, Greenman C, Chen L, Bignell G, Butler A, et al. Systematic sequencing of renal carcinoma reveals inactivation of histone modifying genes. Nature (2010) 463(7279):360–3. doi: 10.1038/nature08672

PubMed Abstract | CrossRef Full Text | Google Scholar

67. Ingelman-Sundberg M, Sim SC, Gomez A, Rodriguez-Antona C. Influence of cytochrome P450 polymorphisms on drug therapies: Pharmacogenetic, pharmacoepigenetic and clinical aspects. Pharmacol Ther (2007) 116(3):496–526. doi: 10.1016/j.pharmthera.2007.09.004

PubMed Abstract | CrossRef Full Text | Google Scholar

68. Simo-Riudalbas L, Esteller M. Cancer genomics identifies disrupted epigenetic genes. Hum Genet (2014) 133(6):713–25. doi: 10.1007/s00439-013-1373-5

PubMed Abstract | CrossRef Full Text | Google Scholar

69. Chou JL, Huang RL, Shay J, Chen LY, Lin SJ, Yan PS, et al. Hypermethylation of the tgf-beta target, Abca1 is associated with poor prognosis in ovarian cancer patients. Clin Epigenet (2015) 7:1. doi: 10.1186/s13148-014-0036-2

CrossRef Full Text | Google Scholar

70. Munoz M, Henderson M, Haber M, Norris M. Role of the Mrp1/Abcc1 multidrug transporter protein in cancer. IUBMB Life (2007) 59(12):752–7. doi: 10.1080/15216540701736285

PubMed Abstract | CrossRef Full Text | Google Scholar

71. Pan ST, Li ZL, He ZX, Qiu JX, Zhou SF. Molecular mechanisms for tumour resistance to chemotherapy. Clin Exp Pharmacol Physiol (2016) 43(8):723–37. doi: 10.1111/1440-1681.12581

PubMed Abstract | CrossRef Full Text | Google Scholar

72. Zappe K, Cichna-Markl M. Aberrant DNA methylation of abc transporters in cancer. Cells (2020) 9(10):2281. doi: 10.3390/cells9102281

CrossRef Full Text | Google Scholar

73. Li Y, Zhao L, Li XF. Hypoxia and the tumor microenvironment. Technol Cancer Res Treat (2021) 20:15330338211036304. doi: 10.1177/15330338211036304

PubMed Abstract | CrossRef Full Text | Google Scholar