Hypoxia-Related Immune Gene Signature Predicting Prognosis in Patients With Hepatocellular Carcinoma

Lingshan Zhou Lanzhou University First A liated Hospital Yuan Yang Lanzhou University First A liated Hospital Min Liu Lanzhou University First A liated Hospital Rong Liu Lanzhou University First A liated Hospital Man Ren Lanzhou University First A liated Hospital Ya Zheng Lanzhou University First A liated Hospital Yuping Wang Lanzhou University First A liated Hospital Yongning Zhou (  yongningzhou@sina.com ) Lanzhou University First A liated Hospital


Abstract
Background Hepatocellular carcinoma (HCC) remains a global health challenge. Increasing evidence indicates that hypoxia is crucial in the evolution and progression of HCC by regulating the tumor immune microenvironment. The present study aimed to construct a prognostic relevant hypoxia-related immune gene (HRIG) signature.

Methods
We analyzed the expression pro le of the 163 HRIGs and clinical information of 371 patients with HCC obtained from The Cancer Genome Atlas (TCGA). Then, consensus clustering analysis was performed to divide HCC patients into clusters 1 and 2 based on the HRIG expression. Subsequently, A multigene signature was constructed by Least absolute shrinkage and selection operator (LASSO) Cox regression analyses. Then, we evaluated the prognostic capability of this signature by Kaplan-Meier analysis, univariate Cox regression and multivariate Cox regression. The prognostic value of the signature was validated in the International Cancer Genome Consortium (ICGC) database. Furthermore, the functional enrichment analyses were preformed to elucidate their biological signi cance. Finally, we evaluated the in ltration of immune cells and the sensitivity of administrating chemotherapeutic agents.

Results
A total of 37 prognosis-related HRIGs were obtained. Subsequently, we constructed an 8-gene signature on the basis of prognosis-related HRIGs, which had a good performance in predicting the overall survival of patients with HCC. In addition, the signature expressed robust when validated in ICGC. The results revealed that these genes involved in some of the HCC-related pathways and was associated with the in ltration of immune cell subtypes. More importantly, the identi ed model was linked to the sensitivity of some chemotherapeutic agents.

Conclusions
HRIG signature is an effective predictor for the prognosis of patients with HCC.

Background
Primary liver cancer is the sixth most common cancer in the world and the third leading cause of cancer-related mortality, accounting for approximately 906,000 new cases and 830,000 deaths in 2020 [1]. Globally, hepatocellular carcinoma (HCC) is the most common form of primary liver cancer and comprises 75%-85% of cases, which pose a serious health burden [1,2]. Many efforts have been made to develop surveillance systems and treatment strategies that include liver resection, radiofrequency ablation, liver transplantation, radiotherapy, systemic treatment and immunotherapy [3]. However, the prognosis of HCC remains poor (5-year survival rate of 5%) due to the high recurrence rate and the lack of diagnostic tools for early detection [4]. More than 75% of HCC patients are diagnosed at advanced stage and curative treatments are not recommended [5]. Furthermore, the recurrence of primary HCC is as high as 70% after clinical resection [6]. Hence, it is imperative to search for novel biomarkers for early detection, predicting the prognosis, and prediction and monitoring for treatment response of HCC.
Many solid tumors exist in an oxygen-starved environment due to the relative inadequate blood supply. Hypoxia occurs in about 50-60% of advanced solid tumors and the potential roles in the pathophysiological process of cancers have attracted lots of interest [7]. To adapt to the low oxygen conditions, the hypoxia inducible factor (HIF) stabilization and downstream signaling are induced, resulting in the numerous changes in gene expression [8]. The related molecular mechanisms of hypoxia play important roles in regulating cell proliferation and survival, induction of angiogenesis, driver of epithelial-tomesenchymal (EMT), increasing therapeutic resistance and reprograming metabolism in cancer [9]. And hypoxia has become one attractive therapeutic target in cancer [10]. In recent years, more and more evidences emphasize the signi cant of hypoxia in the prognostic evaluation of cancers.
Hepatocarcinogenesis is a multistep process driven by liver in ammation and tissue damage. Liver damage disrupts the microvasculature and decreases of the blood ow, leading to the hypoxia situation [11]. Hypoxia microenvironment is common in HCC and has an intimate association with hepatocarcinogenesis, the progression and the recurrence of HCC [12,13]. The key molecular mechanism induced by hypoxia is through the stabilization of transcription factor, hypoxia-inducible factors (HIFs) [14]. And HIF-1α and HIF-2α are highly expressed in HCC and signi cantly associated with poor prognosis of HCC 15. Many researches reveal that HIFs have an intimate link with tumor immunity. HIFs can act as central regulators to damage and inhibit the innate immune response of liver cancer, resulting in immune evasion and creation of immunosuppressive microenvironment [11]. Moreover, HIFs have effect on the adaptive immune system by transcriptionally activate their target genes that are important for the exhaustion of T cell, NK cell and M1 macrophages, and the enrichment of Treg and M2 macrophages [11,16]. However, the relationship between hypoxia-related immune genes and the prognosis of HCC is still far from clear.
In this study, we analyzed the datasets of HCC patients from The Cancer Genome Atlas (TCGA) and identi ed the differentially expressed hypoxia-related immune genes (DEHRIGs) related to the overall survival (OS). Then, a prognostic model with DEHRIGs was constructed to predict the survival potential and validated in the ICGC cohort. We further explored the underlying mechanisms by performing functional enrichment analyses and evaluating the characteristics of tumor immune in the high-and low-risk groups. Last, we screened the candidate agents by exploring the sensitivity of chemotherapeutic agents. The ow diagram of the study was shown in Fig .1.

Data source
The RNA-seq transcriptome data and corresponding clinical data regarding 371 HCC were downloaded from the TCGA project (https://portal.gdc.cancer.gov).

Consensus clustering analysis based on hypoxia-related immune genes (HRIGs)
We assessed the correlation between HRGs and immune-related genes. And the hypoxia-related immune gene (HRIGs) were identi ed, according to the thresholds of correlation coe cients of > 0.4 and P-value of<0.001. In order to elucidate the biological function of HRIGs in HCC, we classi ed the HCC samples into different subtypes with the "ConsensusClusterPlus" (50 iterations and resample rate of 80%).

Construction and validation of a HRIG prognostic signature
The R package limma was performed to screen the differentially expressed HRIGs with the thresholds of adjusted false discovery rate (FDR) P-value of <0.05.
Sequentially, the univariate Cox regression analysis of OS was performed to obtain prognosis-related HRIGs in HCC by the STRING database (https://stringdb.org/). And the genes were identi ed as prognostic-associated genes with P-value of <0.05. Then, we performed the least absolute shrinkage and selection operator (LASSO) algorithm to lter these prognosis-related HRIGs most correlated with OS by the "glmnet" R package. Subsequently, these ltered genes were rolled in the LASSO regression analysis to construct a prognostic signature with penalty parameter (λ) determined by 10-round cross-validation. The risk score of HRIG signature for each patient was calculated: Risk score=Σ n i=1 (expression level of gene x regression coe cient). HCC patients were divided into high-and low-risk groups according to the median risk score determined by using "survminer" R package. Furthermore, we used R package stats and Rtsne to explore the characteristics of different groups by PCA and t-SNE, respectively. Finally, "timeROC" R package was used to assess the predictive potential of the gene signature by measuring the area under the receiver operating characteristic curve (AUC-ROC) [18].

Functional enrichment analyses
Subsequently, we screened the differentially expressed genes (DEGs) between high-risk and low-risk groups with the thresholds of absolute fold change (log2) of >1.5 and FDR of <0.05. To explore the underlying mechanism, we conducted the Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) analyses by using the "clusterPro ler" R package [19].
Evaluation of tumor-in ltrating immune cells single-sample gene set enrichment analysis (ssGSEA) was performed to uncover the characteristics of immune cell in ltration by R package gsva [20]. We calculated the status of immune cell in ltration and relevant immune-related pathways among the HCC samples.

Exploration of the sensitivity of chemotherapeutic agents
To evaluate the sensitivity of chemotherapeutic drugs, we used the "pRRophetic" package in R to measure the half-maximal inhibitory concentration (IC50) of HCC patients in the high-and low-risk groups by ridge regression [21]. Based on the AJCC guidelines, cisplatin, mitomycin.C, vinblastine and sorafenib were selected as candidate antitumor drugs. Finally, we compared the IC50 in different groups by Wilcoxon signed-rank test.

Correlation of consensus clustering for HRIGs with OS of HCC patients
We conducted Pearson correlation coe cient analysis to identify HRIGs. A total 163 HRIGs were screened and rolled in the subsequent differentialexpressionanalysis. Sequentially, consensus clustering analysis was performed based on the HRIG expression pro les. Based on the results with K=2 to 9, K=2 was identi ed as the optimal clustering value, which expressed the good clustering stability ( Fig.S1a-b). Patients with HCC from TCGA cohort were dived into two molecular subgroups (cluster 1 and 2) (Fig.S1c). As shown in the Fig.1, cluster 1 exhibited more favorable prognosis than cluster 2. The results revealed that clustering cancer subtypes charactered by the HRIG expression were involve in the heterogeneity of HCC.

Analysis of differentially expressed HRIGs in HCC and identi cation of prognostic DE HRIGs in TCGA cohort
Among HRIGs, 89 DEHRIGs were identi ed ( Fig.2a-b). Based on the survival data of HCC tumor samples, univariate Cox regression was conducted to identify the prognosis-related DEHRIGs. The results indicated that 37 DEHRIGs were correlated to OS in HCC. GHR, COLEC10 and VIPR1 were downregulated in tumor samples, while others were upregulated (Fig.2c-d). The result of the 37-gene interaction network generated by the STRING database showed that SRC, HSP90AB1, and SHC1 were Hub genes (Fig.2e). The correlation network of these genes was presented in Fig.2f.
Construction and validation of a HRIG prognostic signature By LASSO Cox regression analysis, we built an 8-HRIG prognostic signature from the above 37 OS-related DEHRIGs. The risk score of every patient was computed based on the following formula. Risk score= 0.017176 × S100A10 + 0.226828 × MAPT + 0.021352 × CACYBP + 0.139515 × BIRC5 + 0.064381 × KITLG + 0.043464 × SPP1+ 0.145751 × STC2 + ( -0.02423) × GHR. Based on the median risk score, 365 patients with HCC (after excluding 7 patients without time of follow-up) were divided into high-risk group (N=182) and low-risk group (N=183) (Fig.3a-b). PCA and t-SNE analyses indicated the patients were distributed in two directions on the basis of risk score (Fig.3c-d). The Kaplan-Meier survival analysis showed that the patients in the high-risk group exhibited signi cant shorter OS compared with those in the low-risk group (Fig.3e). The AUC-ROC of the prognostic signature calculated from TCGA was 0.787 at 1 year, 0.724 at 3 years, 0.688 at 5 years, indicating the HRIG prognostic signature had a good performance in monitoring survival (Fig.3f). To veri ed the robustness of the HRIG prognostic signature, we evaluated the predictive value of the identi ed HRIG signature in the ICGC dataset. The risk score of every patient with HCC was calculate based the above formula and the patients were also categorized into high-risk and low-risk groups by the median value ( Fig.4a-b). Similar to the results above, patients in two subgroups from the ICGC cohort were also distributed in two directions ( Fig.4c-d). As shown in Fig.4e, the patients in the high-risk group had a shorter survival time compared with those in the low-risk group, which was consistent with the result from the TCGA dataset. Besides, The AUC-ROC of the prognostic signature calculated from ICGC was 0.740 at 1 year, 0.758 at 3 years, 0.752 at 5 years (Fig.4f).
Review and comparison of recent immune-related gene prognostic signatures constructed for HCC To evaluated the predictive performance and application value of our signature, we further reviewed four recently published prognostic signatures associated with immune-related genes and compared them with this study. Notably, as shown in Table.1, our signature exhibited its superiority. Firstly, it was built by the hypoxia-related immune genes. Compared with other prognostic signatures with single immune genes, it was better able to reveal the complex interaction of various parameters of tumor microenvironment affecting the development and progress of HCC. Additionally, our signature exhibited better performance in predicting OS than the signature built by Pan et al [22]. Compared with other three signatures [23][24][25], HRIG signature expressed more robust when validated in anther cohort. Thus, the study provided a more precise predictor for the prognosis of HCC patients and could provide support in clinical decision-making.

Functional analyses and evaluation of Tumor Immune In ltration
Pathways in KEGG and biological processes in GO were performed to further explore the potential in uence of the classi er based on the risk score in HCC. The analyses indicated that differentially expressed genes between high-and low-risk groups were mainly enriched in immune-related signaling pathways, cancer-related pathways, metabolism-related pathways. (Fig.6a, c), which were validated in the ICGC cohort (Fig.6b, d).
To further explore the in uence of risk score on the tumor immune microenvironment, we evaluated immune in ltration status in HCC samples. In both cohorts, the in ltration of aDCs, DCs, Macrophages, Th2_cells and Treg were signi cantly increased in the high-risk group (P <0.05), suggesting great signi cance of these in ltrating immune cells in the progression of HCC. At the same time, the in ltration of NK cells was signi cantly decreased (P <0.05) ( Fig.7a-b). Subsequently, we investigated the potential correlation between the risk score and immune function. The score of Check-point, HLA and T_cell_costimulation were higher in the high-risk group, while the score of type_I_IFN_Reponse and type_II_IFN_Reponse were the opposite (Fig.7c-d).

Exploration of correlation between the risk model and the e cacy of chemotherapeutic agents
To evaluate the clinical value of the identi ed risk model, we investigated the relationship between risk score and the sensitivity of common administrating chemotherapeutics. The present study demonstrated that the higher risk score was linked to higher chemosensitivity in mitomycin.c and cisplatin, and lower chemosensitivity in vinblastine, respectively (Fig.8a-c). However, there is no connection between the risk score and chemosensitivity in sorafenib (Fig.8d). Thus, these results suggested that the model had potential predictive performance for chemosensitivity.

Discussion
Liver cancer is an extraordinarily heterogeneous malignant tumor, which is one of main obstacles to the implementation of precision medicine [26]. Although the histopathological classi cation of liver cancer has been modi ed and re ned, it still cannot solve the problem that predict patient prognosis or response to therapy accurately [27]. With the development of high throughput sequencing technology, more and more genes with oncogenic functions or tumor suppressive functions were identi ed. The expression patterns of these genes are recurrently altered in HCC, which give tumor cells corresponding abilities to in uence malignant behaviors, such as recurrence, metastasis and drug resistance, leading to different clinical outcomes [28]. One of most important objectives of precision medicine is to provide each cancer patient with most accurate and effective treatment based on the genomic characteristics of the cancer [29]. Without doubt, accurate prognostic biomarkers need to be identi ed for executing precision medicine in a personalized manner. The heterogeneity of tumor exists not only in the genotypes but also in the tumor microenvironment. The accumulation of genomic alterations may result in changes of microenvironment of liver cancer [27]. Immune cells and tissue hypoxia, the elements of tumor microenvironment, have been veri ed to play an important role in the development of HCC and treatment response [30]. In this study, therefore, we focus on the effect of tumor immunity and hypoxia on the prognosis of HCC, which is convincingly bene cial to understanding the mechanism of tumorigenesis and tumor development, and the improvement of personalized therapeutic approaches. Using the RNA-seq transcriptome data and corresponding clinical data of HCC obtained from the TCGA dataset, we constructed the 8hypoxia-related immune gene prognostic signature. The 8-gene signature exhibited a good performance to predict OS of HCC in the training set and expressed robust when validated in another cohort.
Our prognostic signature contained eight genes (S100A10, MAPT, CACYBP, BIRC5, KITLG, SPP1, STC2, GHR). These genes were correlated with OS of HCC separately and showed a better performance when combined, for one reason that gene signature analysis can re ect the complex interaction of various genes affecting hypoxia-related tumor immunity in cancer. S100A10, a member of S100 family, mediates the process of transforming plasminogen to an active protease and has been reported to link to the HCC tumorigenesis and progression [31]. As a binding partner of S100 family proteins, CACYBP was highly expressed in HCC, which was signi cantly associated with elevated AFP level, increased dead and recurrence events, and reduced OS [32]. BIRC5, also known as survivin, is the most potent member of the inhibitor of apoptosis protein (IAP) family and a vital promoter of development and progression of HCC. The overexpression of survivin affects the prognosis of patients with HCC by promoting cell proliferation, inhibiting cell apoptosis and inducing tumor stromal angiogenesis [33]. The osteopontin protein is encoded by the SPP1 gene. Osteopontin overexpression increases proliferation, stem-like properties, glycolysis, and resistance to chemotherapy of HCC cells, which correlates with poor prognosis of HCC patients [34]. Importantly, osteopontin can associate with HIF-1α to promote tumorigenesis [34]. In addition, it plays a critical role in the formation of immunosuppressive microenvironment of HCC [35]. Likewise, STC2 is a proliferation-facilitating gene and correlates with occurrence, development, and prognosis of HCC [36]. Summarily, these genes above have been veri ed to be associated with the poor prognosis of HCC, which was consistent with the result of our study. Moreover, in this model, GHR was positively correlated with OS of the patients with HCC. One study demonstrated that HCC patients with a signi cant down-regulation of GHR expression showed higher incidence of recurrence and poor survival rates [37], which also supported our results. Among these genes, MAPT had the highest coe cient and was considered to be the main contributing gene. Many researches have demonstrated that MAPT was correlated with the prognosis of many cancers. The controversial roles of MAPT in both oncogenesis and tumor suppression were reported in these cancers, depending on the cell type and context. In addition, MAPT is also involved in resistance to paclitaxel, platinum, and bicalutamide, which contributed to the poor prognosis [38,39]. Several lines of evidence have also demonstrated the effect of KITLG in tumorigenesis [40]. However, the biological functions of MAPT and KITLG in HCC have not been elucidated and need further research. Thus it can be seen that the gene signature was closely related to the malignant behaviors of HCC and could be an accurate predictor.
As previously mentioned, multiple pathways were signi cantly enriched in the high-risk group, including ECM-receptor interaction pathway. The signature contained two ECM-related genes, S100A10 and SPP1 [35,41]. As one of the critical components in the tumor microenvironment, dysregulated ECM has been involved in the development and progression of HCC [42]. Hypoxia is one of the key inducers of the dysregulated ECM, by affecting ECM deposition, remodeling and degradation [43]. As shown in our analyses, HCC samples in the high-risk group had higher proportions of macrophages. Interestingly, Hypoxia can recruitment macrophages and induce the alteration of macrophages toward protumor phenotype in the tumor microenvironment, which contribute to ECM remodeling in favor of tumorigenesis and cancer progression [43,44]. At the same time, ECM stiffness through induction of hypoxia is also a barrier for delivery of chemotherapeutic drugs to the tumor site, which weakens the e cacy of antitumor drugs [44]. The high resistance to antitumor drugs is one of the main contributors to the poor outcome of patients with advanced HCC who have no opportunity for surgical resection. The eight genes included in the signature are all linked to the drug resistance. Additionally, in the study, the patients in the high-risk group expressed the lower sensitivity to vinblastine. Apart from the mechanism mentioned above, hypoxia also contributes to the drug resistance by regulating the glycolytic metabolism. As a master regulator of glycolytic metabolism, HIF-1α activation can enhance the glycolytic metabolism at the transcriptional level, which in uences the drug sensitivity in HCC [45]. The glycometabolism-related pathways were enriched in the high-risk group in the study. And glycolytic metabolism also takes roles for the regulation of proliferation, immune evasion, invasion, metastasis, angiogenesis in the liver cancer [46]. More importantly, metabolites released from cancer cells have an impact on the function of immune system cells, such as tumor-associated macrophages, NK cells, dendritic cells (DCs), regulatory T cells (Tregs) and cytotoxic T lymphocytes, leading to inhibition of the immune response [47,48]. Notably, in addition to predicting the OS, the present results suggest that the eight-gene signature has a good performance in predicting the response of HCC patients to antitumor drugs. Therefore, the identi ed signature is a practical tool to help making clinical decision.
In addition to the above-mentioned macrophages, high-risk patients had signi cantly higher proportions of Tregs. Tregs belong in one category of tumourin ltrating lymphocytes and play an important role in immune suppression [49]. Hypoxia has been identi ed in promoting the recruitment of Tregs to the tumour microenvironment by inducing the expression of CCL28, leading to promote tumour tolerance and angiogenesis [50]. In addition, these cells can suppress the functions of T cells, NK cells and DC [49]. As shown in the results, the score of type_II_IFN_Reponse was signi cantly lower in the high-risk group. T cell stimulation under hypoxia conditions inhibits the IFN-γ production [51]. In general, Hypoxia participates in accentuating the immunosuppressive characteristics of immune cells in the tumor environment of HCC, which may have important clinical signi cance for cancer immunotherapy and predicting the prognosis.
To our knowledge, this was the rst study constructing a prognostic signature based on the hypoxia-related immune genes in HCC. Compared with the models based on the genes with single property, the signature contained higher information context. Additionally, the robustness and high predictive value strictly guaranteed the applicability of the signature. However, there were still some limitations in this study. We did the research by utilizing published retrospective data sets. Therefore, our ndings should be validated in prospective studies. Furthermore, some important immune genes with prognostic potential were excludes in this research only based on hypoxia-related immune genes.

Conclusions
We identi ed the hypoxia-related immune genes and constructed an 8-gene prognostic signature through integrated analyses. The present results demonstrated that the HRIG signature was a reliable predictive marker for OS and response to the antitumor drugs in patients with HCC. Thus, the 8-gene signature might be a practical tool for facilitating personalized medicine and improving survival of HCC.

Declarations
Ethics approval and consent to participate Not applicable.

Consent for publication
Not applicable.

Availability of data and materials
The data of this study are available in the TCGA and ICGC.

Competing interests
No con icts of interest exist in this work.       Validation of the HRIG prognostic signature in the ICGC cohort. a. The permutation of risk scores of patients with HCC and identi ed the median risk score. b.
The scatter plot of survival time, survival state and risk score in the ICGC cohort. c. PCA plot showing the result of principal component analysis of the two groups (high-risk group and low-risk group) in the ICGC cohort. d. t-SNE analysis of the two groups (high-risk group and low-risk group) in the ICGC cohort. e.
Kaplan-Meier survival curve of HRIG prognostic signature in the ICGC cohort. f. The AUCs for the risk score was calculated according to the ROC curve in the ICGC cohort. Cox regression analyses for identifying independent predictors of OS in the TCGA cohort and the ICGC cohort. Univariate Cox regression analyses in the TCGA cohort (a) and the ICGC cohort (b). Multivariate Cox regression analyses in the TCGA cohort (c) and the ICGC cohort (d).