Construction of a 3-mRNA hypoxia prognostic model to evaluate immune microenvironment in hepatocellular carcinoma.

BACKGROUND: Hypoxia is a key factor in the development of hepatocellular carcinoma (HCC), which is the most common primary liver cancer with poor prognosis. The current study aimed to identify the potential prognostic biomarkers of the hypoxia-associated gene signature in patients with HCC, and to further explore the relationship between hypoxia and immune infiltration.
METHODS: After the determination of differentially expressed genes (DEGs) using the HCC transcriptome data of The Cancer Genome Atlas database and hypoxia-related gene set, the prognosis-associated genes were identified using univariate Cox regression analysis. Then, the hypoxia prognosis model was established via multivariate Cox regression analysis, with functional annotation conducted using Gene Set Enrichment Analysis. CIBERSORT was utilized to analyze the degree of tumor immune invasion, and an International Cancer Genome Consortium cohort to verify the reliability of the prognosis model. Expression levels of hypoxia-associated genes were detected by real-time quantitative polymerase chain reaction in HCC samples.
RESULTS: 3 genes (ENO1, SAP30, and STC2) constructed the hypoxia prognosis model. The patients were subdivided into 2 groups based on median risk score, with a high hypoxic score indicating poor prognosis of HCC. The hypoxia signature could be employed as an independent prognostic factor in HCC. In addition, the proportion of macrophages was higher in the high-risk group.
CONCLUSION: The hypoxia-associated signature could be a potential prognostic marker of HCC and provides a different perspective for immunotherapy of HCC.

1. Introduction

Liver cancer, 1 of the most common digestive system malignant tumors, ranks the third cause of cancer death in the world, causing >700,000 deaths every year.[ Hepatocellular carcinoma (HCC) accounts for 90% of primary liver cancer, which has become a global public health challenge.[ Despite many treatment methods for HCC, its 5-year survival rate is <19%.[ Low early diagnosis rate, rapid tumor progression, and high recurrence rate may be associated with poor prognosis of HCC patients. However, the mechanism of HCC progression is still unclear. At present, bioinformatics has become 1 of the principal means of cancer research, and it is of profound significance to look for reliable biomarkers to predict the diagnosis, progression, and prognosis of HCC for the prevention and treatment of HCC.

Hypoxia is 1 of the characteristics of solid tumor microenvironment, caused by the insufficient oxygen supply capacity of tumor blood vessels to meet the needs of rapid proliferation of cancer cells.[ Hypoxia has a significant effect on tumor growth, invasion, and metastasis.[ Hypoxia of tumor tissue can induce abnormal vascular regeneration, leading to vascular dysfunction.[ In addition, hypoxia affects the invasion and migration of cancer cells through epithelial–mesenchymal transition (EMT),[ and can also lead to slow proliferation and cell cycle arrest of tumor cells, which reduces the sensitivity of tumor cells to radiotherapy and chemotherapy.[ Studies have shown that high expression of hypoxia inducible factor-1α (HIF-1α) is related to poor prognosis of HCC patients,[ yet the underlying mechanism remains unclear. Immune cells in tumor microenvironment are known as rather important in tumor development,[ while hypoxia leads to immune resistance and immunosuppression, thus helping tumor cells to get rid of immune surveillance.[ Moreover, tumor hypoxia can inhibit the maturation of T cells and rabbit cells and the production of cytokines.[ Therefore, it is necessary to explore the interaction between tumor hypoxia and immunity, and develop a new therapeutic regimen for HCC.

The Cancer Genome Atlas (TCGA) and International Cancer Genome Consortium (ICGC) are 2 cancer research projects, which contain a variety of cancer key genome change profiles, aiming to find potential biomarkers for cancer prevention, diagnosis and treatment through genome-wide data analysis. In this study, TCGA and ICGC databases were used to analyze the genes related to hypoxia in HCC, and the prognosis prediction model of HCC was established. The correlation between hypoxia and immune infiltration in HCC was discussed as well, with a view to providing new ideas for remedy in HCC.

2. Materials and Methods

2.1. Datasets and clinical samples

RNA-seq transcriptome data of 371 HCC patients and the corresponding clinicopathological information were obtained as the training set from TCGA (https://portal.gdc.cancer.gov/). Similarly, the data of 231 HCC patients were downloaded as the validation set from ICGC (https://icgc.org/). Samples without complete clinical information and overall survival <90 days were rejected. Data collection and application conform to TCGA and ICGC publication guidelines and data access policies. Six clinical samples of HCC and adjacent tissues were collected from patients treated in the First Affiliated Hospital of Guangxi Medical University. Six adjacent tissues were taken about 1 cm from the tumor edge. The collected samples were stored in the −80°C refrigerator. The study has obtained the patients’ written informed consent and been approved by the ethics committee of First Affiliated Hospital of Guangxi Medical University. The research was conducted in accordance with the Declaration of Helsinki.

2.2. Identification of differentially expressed genes between HCC and adjacent tissues

Hypoxia-related genomes were extracted from the Molecular Signatures Database (MsigDB) v7.2 (HALLMARK_HYPOXIA M5891, http://www.broadinstitute.org/gsea/msigdb/index.jsp), which contains 200 genes up-regulated in response to hypoxia. Hypoxia-related differentially expressed genes (DEGs) were identified by limma (an R package). DEGs with |log2FC | > 1 and false discovery rate < 0.05 were considered for further analysis.

2.3. Construction and validation of hypoxia model

Prognosis-related hypoxia genes were screened by univariate Cox regression analysis (P < .001). Then, the hypoxia risk signature was established using stepwise regression multivariate Cox analysis. The hypoxia risk signature formula was constructed as follows[:

Ei was the expression of prognostic genes and βi was the regression coefficient of prognostic genes. According to the median risk score, the patients were divided into the low- and high-risk groups. The Kaplan–Meier (KM) curve was plotted to reflect the prognosis of the 2 groups, and was further validated on the ICGC data set.

2.4. Independent prediction analysis and nomogram construction

In order to judge the clinical independence of the prognostic model, Cox regression analysis was performed in combination with other clinical characteristics, such as incorporating the survival time and clinical data integrated into multivariate Cox regression analysis to determine independent prognostic risk factors. The variable with hazard ratio >1 was considered as an inferior prognostic factor. In addition, all independent prognostic factors were used to establish nomograms to assess the 1-, 3-, and 5-year survival rates of HCC patients.

2.5. Relationship of prognostic gene characteristics, hypoxia, and immune cell infiltration

Two gene sets (H: hallmark gene sets; C5.BP: subset of GO) were selected from the Molecular Signatures Database for Gene Set Enrichment Analysis (GSEA) to detect different functional phenotypes between low- and high-risk groups. 1000 genome permutations were performed for each analysis. Phenotypic markers were used as a risk score. Tumor-infiltrating immune cells were estimated by TIMER2.0 data bank (https://cistrome.shinyapps.io/timer/).[ We downloaded the TCGA estimation file from TIMER database, which contains the infiltration level of 6 types of immune cells in different samples. Then, we combined the estimation data with the risk score, and calculated the correlation between the risk score and immune invasion using the cor.test function in R. CIBERSORT (https://cibersort.stanford.edu/),[ a deconvolution algorithm based on gene expression, was used to evaluate 22 immune cell types’ proportions in diverse risk groups. The histogram was used to evaluate the distribution of immune cells in different samples, and the violin chart was used to show the proportion of immune cells in different risk groups.

2.6. Real-time quantitative PCR

According to the manufacturer’s protocol, the total RNA was extracted with RNAiso Plus reagent (Takara, Japan), and then was reverse-transcribed into complementary DNA (cDNA) using PrimeScriptTMRT reagent Kit (Takara, Japan). TB Green Premix Ex TaqTM II Kit (Takara, Japan) was used for Q-PCR in ABI7500 real-time polymerase chain reaction system (Applied Biosystems). PCR was carried out as follows: 30 seconds at 95°C for 1 cycle, then 5 seconds at 95°C, and 34 seconds at 60°C for 40 cycles. The primer sequence is shown in Table 1.

Table 1

The primer sequences of 3 hypoxia-related genes.

ENO1	F primer (5′-3′)	AATGGCGGTTCTCATGCT
	R primer (5′-3′)	ACCTCTGCTCCAATGCG
SAP30	F primer (5′-3′)	CCGCTGTCTAACTTGGTGT
	R primer (5′-3′)	GAAGCCGTTCATGTCTCC
STC2	F primer (5′-3′)	AACTGGGGAAGCCTGTG
	R primer (5′-3′)	GGCTCTGGGAGGTGATG
β-Actin	F primer (5′-3′)	CTACCTCATGAAGATCCTCACCGA
	R primer (5′-3′)	TTCTCCTTAATGTCACGCACGATT

F primer = forward primer; R primer = reverse primer.

2.7. Statistical analysis

R (v.3.6.0) software was used for statistical analysis. Continuous variables were summarized as mean ± SD. Wilcoxon test was performed for comparing differences among the groups. The qualitative variable was tested by Pearson chi-square test. The survival curves were compared by log-rank test. KM curve was drawn for displaying the difference of overall survival (OS).

3. Results

3.1. Identification of hypoxia-associated DEGs in hepatocellular carcinoma

Differentially expressed genes (|log2FC|>1, P < .05) were screened based on the transcriptome data of HCC and noncancerous liver samples from TCGA database. In total, 200 hypoxia marker genes were screened via the MSigDB v7.0. Subsequently, 72 DEGs were identified, which contained 16 down regulated genes and 56 up regulated genes.

3.2. Construction of hypoxia-related prognosis model

Cox regression analysis was performed on 72 differentially expressed hypoxia marker genes in TCGA data set for the establishment of a hypoxic prognostic model of HCC patients. Expression profiles of DEGs were merged with survival data. Univariate Cox regression analysis was performed for revealing 10 hypoxia genes associated with OS. After analyzing the interaction of these genes, a 3-mRNA (ENO1, SAP30 and STC2) model was identified as the best prognostic model for predicting OS by multivariate Cox regression analysis (Table 2). The prognostic model formula was as follows: risk score = (0.0015 × expression level of ENO1) + (0.0960 × expression level of SAP30) + (0.0282 × expression level of STC2). Analysis of the differential expression of 3 genes in the TCGA cohort showed that all genes were highly expressed in tumor, and the same was true in the ICGC cohort.

Table 2

The result of multivariate Cox regression analysis.

Ensemble ID	Gene name	Hazard ratio	P value	Coefficient
ENSG00000074800	ENO1	1.001484	.006108	0.001
ENSG00000164105	SAP30	1.100767	.018702	0.096
ENSG00000113739	STC2	1.028639	.013638	0.028

3.3. Prognostic value of hypoxia signature in hepatocellular carcinoma

Since hypoxia often promotes the malignant phenotype of tumors, we analyzed the relationship between hypoxia signature and the prognosis of tumors. In TCGA and ICGC databases, the risk score of 0.846 divided patients into high- and low-risk groups (Figure S1, Supplemental Digital Content 1, http://links.lww.com/MD/H325). The survival rate of the high-risk group was significantly lower. The data also showed that the expression of 3 hypoxia-related genes increased with the increase of the risk score, which means that high-risk patients tend to form hypoxia microenvironment. In addition, the KM curve displayed prognostic effect of the risk model in HCC. In the TCGA cohort, a high hypoxia score had correlation with poor overall survival (Fig. 1A), which was further confirmed by the ICGC cohort.

Figure 1.

Prognosis evaluation of hypoxic risk signature for overall survival in TCGA and ICGC cohort. (A) Kaplan–Meier curve for the low- and high-risk group. The difference of OS between various risk groups was evaluated via log-rank test. (B) ROC curve analysis for 1-, 3- and 5-year survival prediction. (C) The nomogram constructed to predict 1-, 3- and 5-year survival. ICGC = Cancer Genome Consortium, OS = overall survival, ROC = receiver operating characteristic, TCGA = The Cancer Genome Atlas.

3.4. Prognosis evaluation of hypoxic risk signature

To assess the predictive capability of hypoxia signature in 1-, 3-, and 5-year survival rates, receiver operating characteristic curve was performed using data from TCGA and ICGC datasets. The area under the receiver operating characteristic curve is 0.773 in 1 year, 0.66 in 3 years, and 0.625 in 5 years (Fig. 1B), indicating that the model had good sensitivity and specificity, which is further validated by the ICGC data set.

The prognostic independence of hypoxia signature in HCC was verified by Cox regression analysis. Univariate analysis demonstrated that histological grade (P = .019), pathological stage (P = .001), vascular tumor invasion (P = .032) and 3-mRNA hypoxia signature (P < .001) were significantly correlated with overall survival. After multivariate analysis, hypoxia signature (P < .001) was still an independent factor related to overall survival (Table 3), which confirmed the reliability of our model. All these have been verified by ICGC database (Table 4). Independent prognostic factors were included in nomogram analysis (Fig. 1C), contributing to a more intuitive understanding of the predictive capability of hypoxia signature. Different independent factors were matched with the corresponding points, and the total points were obtained by adding them. Finally, the 1-, 3-, and 5-year survival probabilities of each patient were obtained. A higher total point indicates a worse prognosis in the nomogram.

Table 3

Univariate and multivariate analyses of overall survival in hepatocellular carcinoma patients of TCGA.

Variables	Univariate analysis		Multivariate analysis
	Hazard ratio (95% CI)	P value	Hazard ratio (95% CI)	P value
Age	1.017 (0.995–1.039)	.141	1.010 (0.987–1.033)	.416
Gender	0.698 (0.397–1.227)	.212	0.610 (0.323–1.151)	.127
Histologic grade	1.581 (1.077–2.322)	.019	1.885 (1.151–3.086)	.012
Pathologic stage	1.626 (1.210–2.184)	.001	1.447 (1.024–2.044)	.036
Weight	1.001 (0.988–1.015)	.841	1.009 (0.993–1.025)	.282
Vascular invasion	1.867 (1.054–3.307)	.032	1.619 (0.880–2.980)	.121
AFP	0.974 (0.515–1.842)	.935	0.678 (0.346–1.330)	.258
Prognostic model	1.488 (1.321–1.677)	<.001	1.427 (1.239–1.643)	<.001

CI = confidence interval, TCGA = The Cancer Genome Atlas.

Table 4

Univariate and multivariate analyses of overall survival in hepatocellular carcinoma patients of ICGC.

Variables	Univariate analysis		Multivariate analysis
	Hazard ratio (95% CI)	P value	Hazard ratio (95% CI)	P value
Age	0.999 (0.967–1.032)	.952	0.992 (0.957–1.028)	.673
Gender	0.416 (0.218–0.797)	.008	0.304 (0.155–0.595)	<.001
Pathologic stage	2.159 (1.460–3.195)	<.001	2.083 (1.402–3.095)	<.001
Prognostic model	5.080 (1.799–14.344)	.002	4.093 (1.425–11.757)	.009

CI = confidence interval, ICGC = International Cancer Genome Consortium.

3.5. Correlation between hypoxia gene expression and clinicopathological features of HCC

Given the pathophysiological significance of hypoxia in tumorigenesis and invasion, the relationship between 3 hypoxia genes and the clinicopathological characteristics of HCC was analyzed, including pathological stage, histological grade and vascular tumor invasion. Based on the analysis of TCGA and ICGC data sets, the expression levels of ENO1 and SAP30 increased gradually in later pathological stage (Fig. 2), and the expression level of ENO1 was significantly higher in HCC with more advanced histological grade and vascular tumor invasion.

Figure 2.

The relationship between 3 hypoxia-related mRNAs and clinicopathological features. Data were analyzed using the Kruskal–Wallis test. ** represents P < .01; *** represents P < .001. G = histological grade, Macro = macrovascular invasion, Micro = microvascular invasion.

3.6. Hypoxia-related signal pathways

GSEA was performed for identifying pathways that were significantly enriched in the high-risk group. We discovered that some pathways promoting tumor progression and inhibiting apoptosis were significantly enriched in the high-risk group, like DNA repair, hypoxia and PI3K-Akt-mTOR signaling (Fig. 3A). This indicated that a high hypoxia score was beneficial to boost the biological characteristics of tumor cells. These are further verified in ICGC database (Fig. 3B). The schematic diagram of hypoxia-related signaling pathway is shown in Figure S2, Supplemental Digital Content 2, http://links.lww.com/MD/H326,

Figure 3.

GSEA enrichment analysis. (A) GSEA displayed prominent enrichment of tumor progression related pathways in the high-risk group based on the TCGA cohort. (B) These pathways were verified in the ICGC cohort. GSEA = gene set enrichment analysis, ICGC = Cancer Genome Consortium, ns = not significant, TCGA = The Cancer Genome Atlas.

3.7. Immune landscape of patients with low and high risk of hypoxia in HCC

We also predicted the relation of hypoxia and immune-associated prognostic models and the degree of immune cell infiltration in TCGA HCC cases to test whether the risk score partly reflects the state of TIMER. Interestingly, the findings of this study showed that the levels of macrophages, neutrophils, CD8+ T cells and dendritic cells (DC) were significantly positively correlated with risk score (Fig. 4). Furthermore, CD4+ T cells and B cells also showed a weak relationship with risk score. In addition, CIBERSORT method combined with LM22 characteristic matrix was employed to estimate the invasion of 22 immune cells in different risk groups. Analysis outcomes of TCGA and ICGC cohorts were exhibited by bar chart (Fig. 5A). The TCGA cohort shows that 8 immune cell types (resting CD4+ memory T cells, activated CD4+ memory T cells, follicular helper T cells, activated NK cells, monocytes, M0 macrophases, resting mast cells and neutrophils) are significantly different between high- and low-risk groups. While the ICGC cohort indicates significant differences in 6 immune cell types (naive B cells, naive CD4+ T cells, activated CD4+ memory T cells, regulatory T cells (Tregs), M0 macrophages and resting dendritic cells) between different risk groups (Fig. 5B). Compared with the low-risk group, the proportion of M0 macrophages was higher, whereas the proportion of activated CD4+ memory T cells was lower in the high-risk group.

Figure 4.

Relationships between the hypoxia-related mRNAs signature and infiltration abundance of immune cells. Data were analyzed using Pearson’s correlation analysis.

Figure 5.

Immune landscape of patients with low and high risk of hypoxia in HCC. (A) Proportion of immune cell infiltration in high- and low-risk groups. (B) Violin plot showing the differential expression of immune cells in high- and low-risk groups. The blue represents the low-risk group and the red represents the high-risk group. Data were analyzed using the Wilcoxon test. HCC = hepatocellular carcinoma.

Then, we investigated the relationship between hypoxia signature and immune regulation. GSEA was carried out for analyzing the immune-related enrichment pathways in different risk groups. The results indicated that some negative regulatory pathways of the immune process were significantly enriched in the high-risk group, for instance, negative regulation of CD4+ αβ T cell activation, immune system process, T cell differentiation and αβ T cell activation (Figure S3, Supplemental Digital Content 3, http://links.lww.com/MD/H327).

Hence, targeted hypoxia could have important clinical significance in ameliorating immunotherapy.

3.8. Immunosuppressive microenvironment in the high-risk group

The cancer immune cycle has become an important focus of tumor immunotherapy. It involves a series of antitumor immune response procedures, from the release of tumor cell antigens to the killing of tumor cells by T-cell. In this study, we focused on genes that negatively regulate this process in different risk groups. The gene signature was downloaded from TIP website (http://biocc.hrbmu.edu.cn/TIP/index.jsp). Most of negative regulatory genes in tumor immune response were up-regulated in the high-risk group (Figure S4A, Supplemental Digital Content 4, http://links.lww.com/MD/H328), which indicated the immune process of this group was strongly inhibited.

In view of past research that hypoxia promotes the expression of immunosuppressive cytokine and checkpoint, the expression of these genes in different risk groups was included in our study. The outcome demonstrated TIM-3, which had a positive correlation with the hypoxia signature, was up-regulated in high-risk patients (Figure S4B, Supplemental Digital Content 4, http://links.lww.com/MD/H328). In addition, the expression of checkpoints (such as CTLA-4, TIGIT, and programmed cell death 1) was more augmented in patients with high hypoxia scores (Fig. 6A), as did immunosuppressive cytokines (Fig. 6B).

Figure 6.

Immunosuppressive microenvironment in the high-risk group. (A) Expression of checkpoints in high- and low-risk groups. (B) Expression of immunosuppressive cytokines in high- and low-risk groups. Data were analyzed using an unpaired t test. (C) Verification of mRNA levels of 3 hypoxia prognostic genes. Data were analyzed using a paired t test.* represents P < .05; ** represents P < .01; *** represents P < .001. ns = not significant.

3.9. Clinical validation based on mRNA levels of 3 genes

We analyzed 6 pairs of HCC and para-cancerous tissues to verify the mRNA levels of 3 genes. The results showed that all 3 mRNAs were highly expressed in tumor tissues (P < .05), which was consistent with our data analysis results (Fig. 6C).

4. Discussion

As a significant feature of solid malignant tumors, hypoxia tumor microenvironment is considered to be an independent prognostic factor, which is associated with poor prognosis in a variety of malignant tumors, including HCC.[ Hypoxia was reported to be related to the infiltration and metastasis of HCC,[ 1 of the most hypoxic tumor types with an average oxygen level of 0.8%.[ The complexity of HCC tumor microenvironment has a significant correlation with the development of HCC. Activation of immunosuppressive cells and immune checkpoints promotes the immune escape of HCC.[ In addition, hypoxia can further stimulate angiogenesis and immunosuppressive cell expression, thus promoting the progression of HCC.[ However, the role of hypoxia in the progression and prognosis of HCC has not been fully studied. In recent years, transcriptome sequencing has become an important means in biomedical exploration to identify biomarkers in predicting tumor prognosis.[

The HCC prognostic model established in our research consists of 3 hypoxia-related genes, that is, ENO1, SAP30, and STC2. ENO1, a key glycolytic enzyme closely related to the “Warburg effect” of cancer cells,[ has been reported to be overexpressed in a variety of cancers.[ ENO1 promotes the proliferation, metastasis and diffusion of cancer cells by playing the role of plasminogen receptor,[ and can also induce the immune response in cancer patients.[ STC2, as a HIF-1α target gene, promotes the proliferation of human breast cancer and ovarian cancer cells under hypoxia.[ Some reports indicate that STC2 is up-regulated in many cancers, including renal cell carcinoma, neuroblastoma, and HCC.[ In addition, a significant increase in SAP30 mRNA level was observed in clear cell renal cell carcinoma (ccRCC).[ We used 6 pairs of HCC tissues and para-cancerous tissues to test the mRNA levels of the 3 genes. The results showed that all 3 mRNAs were highly expressed in tumor tissues. Our experimental results were basically consistent with those of TCGA and ICGC databases, which also confirmed the reliability of our model.

At present, a few prognostic models composed of multiple genes could predict the prognosis of HCC, such as the prognostic model composed of 8 immune-related mRNAs developed by Xu et al.[ However, it seems that this model is not widely used in the study of hypoxia and immune microenvironment of HCC. Our results show that the risk model composed of these 3 genes can distinguish high-risk population from the low-risk population, and accurately predict prognosis. We found the prognosis of HCC patients with high hypoxia scores was worse based on TCGA and ICGC data. It was worth noting that the risk model was an independent prognostic risk factor in HCC. Based on this signature, we constructed a prognostic nomogram to help patients formulate short-term treatment strategies. In terms of clinical relevance, ENO1was significantly associated with the pathological stage, histological grade and vascular tumor invasion of HCC patients in the TCGA cohort. The same trend was observed in the ICGC queue, demonstrating the clinical predictive value of our risk model. In the enrichment analysis, GSEA confirmed that there were more hypoxia-related reactions in the high-risk group. Therefore, this prognostic model could be involved in the development of HCC, and may become a reliable clinical biomarker.

Hypoxia plays an important role in anticancer immune response, including reducing the activity of effector cell like CD4+ cell and natural killer cell, promoting the activity of immunosuppressive cells such as regulatory T cell and M2 macrophages, promoting immunosuppressive cytokine’ synthesis, and heightening the expression of checkpoint inhibitors to protect tumor from the influence of antitumor immune response.[ Under hypoxia, tumor cells can promote the secretion of Semaphorin 3A (Sema3A) and the recruitment of macrophages.[ In addition, hypoxia recruits tumor-associated neutrophils through chemokines CXCL1 and CXCL2.[ The hypoxia signature we established was found to reveal a positive correlation with infiltration of macrophages and neutrophils. Consistent with the result of vitro experiments that hypoxia could inhibit the proliferation and activation of T cells,[ it was found that the proportion of activated CD4+ memory T cells was less in patients with high hypoxic scores, while M0 macrophages’ infiltration was more obvious by CIBERSORT. Tregs cells promote immunosuppression of tumor microenvironment by producing TGF-β and suppressor effector T cells. HIF-1α promotes Tregs recruitment in tumor microenvironment by up-regulating CCL28 under hypoxia.[ However, our study showed that only the ICGC cohort showed significant differences in Tregs cells between the high- and low-risk groups, which may be related to differences in populations from different regions.

Cytokines are essential in adjusting tumor immunity. In tumor microenvironment, TGF-β inhibits the immune response by reducing the activity of Th1 T lymphocytes.[ In addition, hypoxia can activate TGF-β and weaken the function of NK cells.[ Up-regulation of IL-10 under hypoxia conditions can inhibit the differentiation and maturation of DC, thereby inhibiting the activation of T cells.[ In our study, the expression of immunosuppressive cytokines increased in high-risk patients, further inhibiting immune response.

Cancer cells can avoid the recognition and destruction of the immune system through immune checkpoints. Research has suggested that hypoxia could up-regulate the expression of programmed cell death 1 (PD-L1) in DC and tumor-infiltrating macrophages.[ In addition, HIF-1α promotes the expression of immune checkpoint PD-1, LAG-3, and CTLA-4 in CD8+ T cells, leading to T cell failure.[ This was similar to our results that the expression of immune checkpoints like PD-L1, CTLA4, TIGIT, and TIM-3 was significantly enhanced in patients with high hypoxia scores.

The advantage of our research is that a new hypoxia signature was constructed using the TCGA cohort for predicting the prognosis of HCC patients, which has been verified in ICGC database. The model presented powerful predictive capability in the overall survival of HCC patients, reflecting the immune status of HCC tumor microenvironment. However, there are still some limitations in this study. First of all, although the external database has been verified, further verification by proteomics and immunohistochemistry is required. Second, our retrospective analysis still needs to be verified in a prospective analysis. Finally, the interaction between hypoxia and immunity in tumor microenvironment needs to be further explored. In the future, the prognostic value of the 3-mRNA signatures will need to be further verified in larger clinical patients.

5. Conclusion

In conclusion, this study identified 72 hypoxia-related DEGs through differential analysis. Then, we constructed a prognostic model based on the expression of 3 prognostic mRNAs. This model can significantly distinguish the high- and low-risk groups of HCC patients, and the prognosis of patients in the high-risk group is worse. Cox regression analysis showed that the prognostic model was an independent prognostic factor for HCC and could effectively predict the prognosis of HCC patients. In addition, this prognostic model is related to the clinical progression of HCC. Functional enrichment analysis showed that hypoxia-related pathway and PI3K-Akt-mTOR signaling pathway were significantly correlated with the risk score, indicating that prognostic signature may be involved in regulating metabolism and proliferation of HCC cells. In addition, the outcomes suggested high hypoxia risk group could inhibit the immune microenvironment by up-regulating immunosuppressive cytokines and immune checkpoints. These results will help to provide new insights for hypoxia targeted and immunotargeted therapy of HCC.

Acknowledgements

The authors are grateful to the study participants.

Author contributions

Conceptualization: Wang Jue, Jin Zongrui.

Data curation: Wu Guolin.

Formal analysis: Guo Ya.

Funding acquisition: Wang Jilong, Wen Zhang.

Investigation: Wu Guolin, Zhu Hai.

Methodology: Wang Jue, Jin Zongrui, Deng Zhenfeng.

Project administration: Zhu Hai, Guo Ya.

Resources: Wang Jilong, Guo Ya.

Software: Deng Zhenfeng.

Supervision: Wen Zhang.

Validation: Xu Banghao, Zhu Hai.

Visualization: Xu Banghao.

Writing – original draft: Wang Jue.

Writing – review & editing: Wen Zhang.