Dissecting a hypoxia-related angiogenic gene signature for predicting prognosis and immune status in hepatocellular carcinoma

Background: Hypoxia and angiogenesis, as prominent characteristics of malignant tumors, are implicated in the progression of hepatocellular carcinoma (HCC). However, the role of hypoxia in the angiogenesis of liver cancer is unclear. Therefore, we explored the regulatory mechanisms of hypoxia-related angiogenic genes (HRAGs) and the relationship between these genes and the prognosis of HCC.

Methods: The transcriptomic and clinical data of HCC samples were downloaded from public datasets, followed by identification of hypoxia- and angiogenesis-related genes in the database. A gene signature model was constructed based on univariate and multivariate Cox regression analyses, and validated in independent cohorts. Kaplan-Meier survival and time-dependent receiver operating characteristic (ROC) curves were generated to evaluate the model’s predictive capability. Gene set enrichment analysis (GSEA) was performed to explore signaling pathways regulated by the gene signature. Furthermore, the relationships among gene signature, immune status, and response to anti-angiogenesis agents and immune checkpoint blockade (ICB) were analyzed.

Results: The prognostic model was based on three HRAGs (ANGPT2, SERPINE1 and SPP1). The model accurately predicted that low-risk patients would have longer overall survival than high-risk patients, consistent with findings in other cohorts. GSEA indicated that high-risk group membership was significantly associated with hypoxia, angiogenesis, the epithelial-mesenchymal transition, and activity in immune-related pathways. The high-risk group also had more immunosuppressive cells and higher expression of immune checkpoints such as PD-1 and PD-L1. Conversely, the low-risk group had a better response to anti-angiogenesis and ICB therapy.

Conclusions: The gene signature based on HRAGs was predictive of prognosis and provided an immunological perspective that will facilitate the development of personalized therapies.

Introduction

Globally, liver cancer was the sixth-most common cancer (accounting for 4.5% of all new tumors) and the third leading cause of cancer-related death (accounting for over 8.3% of all such deaths) in 2020 (1). Among the primary types of liver cancer, hepatocellular carcinoma (HCC) accounts for approximately 75–85% of all cases, and the 5-year survival rate of HCC in China is only 14.1% (1, 2). Ablative therapies, transarterial chemoembolization (TACE), and surgery are routine treatments for HCC (3). Because of the complex etiology and high heterogeneity of HCC, its treatments and prognosis are unsatisfactory, and prognosis prediction is challenging (4–6). To improve prognosis, identifying biomarkers for treatment and prognostic prediction of HCC is critical.

Due to tumor neovascularization and high metabolism, hypoxia is present in approximately half of solid tumors, including HCC (7–10). Hypoxia promotes tumor cell proliferation and tumor progression by initiating multiple adaptive behaviors, such as angiogenesis, proliferation, and invasion. Hypoxia reconstructs the tumor immune microenvironment (TIM) by promoting the recruitment of innate immune cells, and interfering with the differentiation and function of adaptive immune cells (11–13). As a solid tumor rich in blood vessels, tumor-driven hypoxia increases the expression of proangiogenic factors leading to abnormal vascular proliferation of HCC, which contributes to tumor growth, invasion, metastasis, and immunosuppression. Tumor susceptibility to angiogenesis has been a hot spot of research recently, and targeted drugs are in clinical use (14, 15). During cancer progression, several proangiogenic cytokines—such as vascular endothelial growth factor A (VEGFA), fibroblast growth factor (FGF) and hypoxia inducible factor-1 (HIF-1)—which contribute to neovasculature sprouting and formation in the tumor microenvironment (TME), are implicated in immune TME remodeling and direct or indirect immune cell regulation (16–18). Nevertheless, the role of hypoxia in liver cancer angiogenesis is unclear. Closely related to angiogenesis, hypoxia regulates a series of genes involved in tumor angiogenesis, leading to the epithelial-mesenchymal transition (EMT) and immune escape, rendering tumor cells more tolerant to the hypoxic microenvironment, and enhancing their proliferation, metastasis, and invasion. Therefore, we analyzed the effect of angiogenesis-related genes under hypoxic conditions on the survival and immune microenvironment of HCC.

We first downloaded mRNA expression profiles and the corresponding clinical data of patients with HCC from public databases. Second, we constructed a prognostic multigene signature based on hypoxia-related angiogenic genes (HRAGs) in The Cancer Genome Atlas (TCGA) cohort, and validated it in an independent cohort. Third, we performed functional enrichment analysis and immune infiltration, as well as anti-angiogenesis and immune checkpoint blockade (ICB) therapy response predictions, to provide guidance for precise and effective HCC treatment.

Materials and methods

Patients and datasets

The mRNA expression profiles and corresponding clinical information of HCC patients were obtained from The Cancer Genome Atlas-Liver Hepatocellular Carcinoma dataset (TCGA-LIHC, https://portal.gdc.cancer.gov/, 370 HCC and 50 normal tissue samples), International Cancer Genome Consortium (ICGC, https://dcc.icgc.org/projects/LIRI-JP, 229 HCC tissue samples), and Gene Expression Omnibus (GEO, https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE14520, 220 HCC tissue samples). The inclusion criteria for the follow-up analysis were HCC confirmed by pathology, available RNA expression data, and complete clinical data (> 30 days of follow-up). In total, 337 hypoxia-related genes (HRGs) and 201 angiogenesis-related genes (ARGs) were acquired by gene set enrichment analysis (GSEA) (hallmark-hypoxia or hallmark-angiogenesis), as well as from the GeneCards database (using the terms hypoxia and angiogenesis; relevance scores > 4) and from previous reports.

Development and validation of a prognostic gene signature

The R package limma was used to identify differentially expressed genes (DEGs) between tumor tissues and adjacent nontumorous tissues, according to the criteria of | log 2 (fold-change) | > 1 and false discovery rate (FDR)< 0.05. Next, hypoxia-related angiogenic DEGs (HRAGs) were identified. The Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) functional annotations of these HRAGs were analyzed and visualized using the R package clusterProfiler.

Univariate Cox regression analysis was performed using the R package survival to determine the prognostic value of the DEGs for OS; P< 0.05 was considered statistically significant. An interaction network for the prognostic DEGs was generated using the STRING database (https://www.string-db.org/) and entered into a stepwise multivariate Cox regression analysis to identify covariates with independent prognostic value for OS. The risk score was based on the expression of predictive genes and the multivariate Cox regression risk model coefficients, and was calculated as follows:

$$\text{Risk\hspace{0.17em}score\hspace{0.17em}}={\displaystyle {\sum }_{i=1}^n(GeneExpressioni\times Coefi)}$$

Based on the median risk score or optimal cut-off value, HCC patients were divided into high- and low-risk groups. A Kaplan–Meier (K-M) survival analysis was performed to compare the high- and low-risk groups according to predictive signatures. The utility of prognostic prediction models was evaluated by calculating the areas under the curve (AUC) values of the receiver operator characteristic (ROC) curve using the R package survivalROC. Principal component analysis (PCA) and t-distributed stochastic neighbor embedding (t-SNE) were used to examine the clustering of the signature genes with the prcomp function of the R package stats and Rtsne. The ICGC and GSE14520 cohorts were analyzed to verify the results.

Prognostic value of the gene signature and construction of a predictive nomogram

Univariate and multivariate Cox regression analyses were performed to identify factors independently associated with OS in the TCGA cohort. To predict the survival of HCC patients, a nomogram was established based on the risk score and other clinical parameters. Calibration curve and time-dependent ROC analyses were performed to assess the accuracy of the nomogram. The nomogram and calibration curves were plotted using the R package rms.

Molecular characteristics and biological function analysis

GSEA of the MSigDB Collection (h.all.v7.4.symbols.gmt) was performed using GSEA software (version 4.1.0) to detect the set of genes expressed difference between the high- and low-risk groups. For each analysis, 1,000 gene set permutations were performed.

Estimation of tumor immune microenvironment

xCell (https://xcell.ucsf.edu/), which is based on a deconvolution algorithm, was used to infer immune cell infiltration from RNA-sequencing data. The infiltration scores of 17 immune cells and activities of 13 immune-related pathways were evaluated by single sample gene set enrichment analysis (ssGSEA) using the R package gsva. The normalized gene expression data of the TCGA and ICGC cohorts were uploaded into Sangerbox tools (http://www.sangerbox.com/tool) for bioinformatics analysis. Next, infiltration of cancer-associated fibroblasts (CAFs) was estimated using the MCPcounter and EPIC algorithms.

Anti-angiogenesis and ICB therapy response prediction and validation

Sensitivity and resistance to anti-angiogenesis drugs were evaluated using the R package pRRophetic. ImmuCellAI was used to predict the response to Immune checkpoint blockade (ICB) therapy (http://bioinfo.life.hust.edu.cn/ImmuCellAI#!/) (19). Immunohistochemical staining was performed as described previously. Briefly, HCC sample slides were deparaffinized and dehydrated by serial immersion in xylene, ethanol, and distilled water. Antigen retrieval was performed using citrate buffer (10 mM, pH 9.0) in a microwave on medium power. After blocking with goat serum at room temperature for 1 h, the tissues were sequentially incubated with an anti-secreted phosphoprotein-1 (SPP1) antibody (1:100, 22952-1-AP; Proteintech, Rosemont, IL, USA) at 4°C overnight and a secondary antibody at room temperature for 1 hour. A slide scanner (3DHistech Ltd., Budapest, Hungary) was used to capture images. Staining intensity, expressed as the H-score (range: 0–300) was automatically quantified using Pannoramic Viewer software (3DHistech Ltd.). Approval was obtained from the Institutional Review Boards of the Research Institute and Hospital National Cancer Center and The First Affiliated Hospital, Sun Yat-Sen University.

Statistical analysis

Data management and statistical analysis were conducted using R (version 4.1.0; R Development Core Team, Vienna, Austria) and GraphPad Prism software (version 8.3.0; GraphPad Software Inc., San Diego, CA, USA). Student’s t-test or the Wilcoxon rank-sum test was used to compare gene expression between two groups. The chi-squared test was used to compare differences in proportions. Survival curves were plotted using the K-M method and compared by log-rank test. A value of P< 0.05 was taken to indicate statistical significance.

Results

Datasets

The study flow chart is shown in Figure 1. A total of 343 HCC patients from TCGA cohort, 229 from the ICGC cohort, and 220 from the GSE14520 cohort were included in this study. The clinical characteristics of the patients are listed in Supplementary Table S1.

FIGURE 1

Figure 1 Flow chart of data collection and analysis.

Identification of prognostic hypoxia-related angiogenic genes in the TCGA dataset

DEGs analysis of the 374 HCC samples and 50 normal liver samples from TCGA revealed 21 differentially expressed HRAGs (Figures 2A–C). GO and KEGG enrichment and functional analysis showed that these genes were enriched in epithelial cell migration, chemotaxis, PI3K-Akt signaling pathway, focal adhesion, MAPK signaling pathway, and HIF-1 signaling pathway (Figures 2D, E). The HRAGs were related to cancer proliferation and invasion. Univariate Cox regression analysis showed that nine prognostic HRAGs significantly correlated with OS were risk factors for a poor prognosis of HCC (Figure 3A). PPI and gene correlation networks suggested connections among the prognostic genes (Figures 3B, C).

FIGURE 2

Figure 2 Screening and enrichment analysis of hypoxia-related angiogenesis genes. (A) Venn diagram of 21 differentially expressed HRAGs in the TCGA cohort. (B) Deviation plot and (C) heatmap of 21 differentially expressed HRAGs in HCC and noncancerous tissues. (D) GO (E) and KEGG analyses revealed the most significantly enriched biological functions and pathways of the overlapping differentially expressed HRAGs.

FIGURE 3

Figure 3 Identification of candidate genes related to the prognosis of HCC. (A) Forest plots of univariate Cox regression analysis of gene expression and OS in the TCGA cohort. (B) Protein–protein interaction (PPI) network and (C) correlated regulation network of prognostic HRAGs. (D) Correlation analysis of three model HRAGs from GeneMANIA.

Construction and validation of the prognostic model

Three of the nine prognostic HRAGs were used to construct a prognostic predictive model. The corresponding coefficients and gene expression levels were used to calculate the risk score as follows: (0.389932838 × expression level of ANGPT2 + 0.108256217 × expression level of SERPINE1 + 0.121709881 × expression level of SPP1). Patients in the TCGA, ICGC, and GSE14520 datasets were divided into high- and low-risk groups based on the median risk score or optimal cut-off value. Interestingly, K-M survival curves based on the three genes showed that the predicted OS of the low-expression group was significantly longer than that of the high-expression group (Supplementary Figure S1). A PPI network was generated based on coexpression information for the three candidate genes from the GeneMANIA database (http://genemania.org/) (Figure 3D). The expression of prognostic genes and risk scores were significantly higher for tumor stage III/IV compared with stage I/II patients (Supplementary Figure S2). Therefore, a higher risk score is associated with more malignant HCC.

Furthermore, the high-risk group membership was significantly associated with higher tumor grade and advanced tumor, node, metastasis (TNM) stage (Table 1). Patients with a high risk score had a higher mortality rate, shorter OS, and higher expression of the three model genes (Figure 4A). K-M survival curves indicated significantly longer survival time of the low-risk group, including OS (Figure 4B) and disease-free survival (DFS) (Supplementary Figure S3), compared to the high-risk group. The AUC of the time-dependent ROC curves was 0.783 at 0.5 years, 0.746 at 1 year, 0.733 at 1.5 years, 0.668 at 2 years, 0.648 at 3 years, and 0.669 at 5 years (Figure 4C). PCA (Figure 4D) and t-SNE (Figure 4E) confirmed risk profile differences between the two groups. To test the robustness of the gene signature model constructed based on the TCGA cohort, patients from the ICGC and GSE14520 cohorts were categorized into high- and low-risk groups according to the optimal cut-off value. The results for the ICGC (Figures 4G–J) and GSE14520 cohort (Supplementary Figure S4) were similar to those of TCGA cohort. Therefore, the prognostic model was predictive of the prognosis and progression of HCC.

TABLE 1

Table 1 Baseline characteristics of the patients in different risk groups of three cohorts.

FIGURE 4

Figure 4 Survival analysis of HCC patients in the TCGA training and ICGC validation datasets. (A, B) Risk score distribution, survival status, and heatmap of the expression of the three HRAGs in the high- and low-risk groups in the training and validation cohorts. (B, G) Kaplan-Meier OS curve. (C, H) AUCs of time-dependent ROC curves. PCA (D, I) and t-SNE (E, J) analysis confirmed the clustering of the three genes comprising the HRAG signature. (A–E) TCGA cohort; (F–J) ICGC cohort.

Prognostic value of the gene signature

Univariate and multivariate Cox regression analyses showed that TNM stage and risk score were significantly associated with the prognosis of HCC in the TCGA cohort (Figures 5A, B). ROC curve analysis showed that the risk score was better for predicting prognosis than the other clinicopathological factors (Figure 5C). The risk scores and clinicopathological factors of 343 HCC patients with complete clinical information were used to create a prognostic nomogram for predicting survival (Figure 5D). The calibration curves of the prognostic nomogram showed good consistency between the predicted and actual 1-, 2-, 3-, and 5-year survival rates in the TCGA cohort (Figure 5E). The nomogram AUCs of the predicted 1-, 2-, 3-, and 5-year OS were 0.777, 0.712, 0.721, and 0.728, respectively (Figure 5F).

FIGURE 5

Figure 5 Independent prognostic power of the three genes comprising the signature in the TCGA cohort. (A) Univariate and (B) multivariate Cox regression analyses of the associations of the risk index (RI) and clinical parameters with OS. (C) ROC curve of the prognostic utility of the risk score, age, gender, grade, and TNM stage. (D) Nomogram for predicting 1-, 2-, 3-, and 5-year survival. (E) Calibration curves of the nomogram showed consistency in the predicted and observed 1-, 2-, 3-, and 5-year survival rates. (F) ROC curve analysis of the nomogram for 1-, 2-, 3-, and 5-year OS.

Molecular characteristics and biological function analysis

GSEA analysis revealed that high-risk group membership was significantly associated with angiogenesis, EMT, glycolysis and hypoxia. The immune-related pathways IL2/STAT5 and IL6/JAK/STAT3, as well as the inflammatory response, interferon (IFN)-γ/response and tumor necrosis factor (TNF)-α signaling via NFKB were significantly enriched in the high-risk group (P< 0.05, FDR< 0.25) (Figure 6).

FIGURE 6

Figure 6 Molecular characteristics and biological function analysis. (A, B) Gene set enrichment analysis of biological functions and pathways in the risk score groups (P< 0.05, FDR< 0.25).

Tumor immune microenvironment analysis

xCell was used to evaluate the proportions of infiltrating immune cells. Patients in the high-risk group had higher proportions of macrophages and T regulatory cells (Tregs) than those in the low-risk group, but there was no difference in the proportions of other immune cell types (Supplementary Figures S5A, B). Next, the enrichment scores of immune cell subpopulations, related functions, and pathways were calculated by ssGSEA (Figure 7). Interestingly, dendritic cells (DCs), antigen-presenting cells (APCs), human leukocyte antigen (HLA), and major histocompatibility complex class I (MHC class I), which are involved in antigen presentation, were significantly elevated in the high-risk group. Furthermore, the activities of cytokine-cytokine receptors and immune checkpoints, and the enrichment scores for macrophages, myeloid-derived suppressor cells (MDSCs), and Tregs in the high-risk group were higher than in the low-risk group, whereas the type II IFN response showed the opposite trend. The risk score was significantly positively correlated with the proportion of cancer-associated fibroblasts (CAFs). The expression of immunosuppressive genes in the high-risk group was higher than in the low-risk group, in both the TCGA and ICGC cohorts (Supplementary Figures S5C, E). The expression of the pro-angiogenic factors, VEGFA and VEGFB, was significantly upregulated in the high-risk group (Figure S5D, F). However, there was no difference in the tumor mutation burden (TMB) or microsatellite instability (MSI) score (Supplementary Figures S5G, H). Therefore, the TME of patients in the high-risk group is in an immunosuppressive state.

FIGURE 7

Figure 7 Tumor immune microenvironment analysis. (A, E) ssGSEA scores of 17 immune cells and 13 immune-related functions. ns: not significant; *P< 0.05; **P< 0.01; ***P< 0.001. (B, F) Spearman’s rank correlation analysis between risk score and immune checkpoints, macrophages, MDSCs, Tregs, and type II IFN response. (C, D, G, H) Spearman’s rank correlation analysis between risk score and CAFs according to the (C, G) EPIC and (D, H) MCPcounter algorithms. (A–D) TCGA cohort; (E–H) ICGC cohort.

Prediction and validation of the anti-angiogenesis and ICB responses

As shown in Figure 8A, the estimated IC₅₀ showed that the low-risk group in both cohorts had a better response to sorafenib (P< 0.05). ImmuCellAI showed that the response rate to ICB therapy was higher in the group with a lower risk score (Figure 8B). The correlation heatmap indicated that SPP1 expression was most relevant to the risk score, and thus has great potential for predicting the response to anti-angiogenesis and ICB therapy (Figure 8C). Next, 19 Barcelona Clinic Liver Cancer (BCLC) C stage-HCC samples subjected to anti-angiogenesis treatment and immunotherapy after resection were further subjected to IHC staining. Follow-up information was collected from January 2015 to April 2022. Based on the H-scores, we divided the 19 patients into high- and low-expression groups. Representative immunohistochemically stained images are shown in Figure 8D. K-M curves showed that the OS of the low-expression group was trend longer than that of the high-expression group though not statistically significant (P = 0.1519, Figure 8E). In summary, the gene signature was predictive of the response to anti-angiogenesis treatment and ICB therapy.

FIGURE 8

Figure 8 Therapeutic response prediction and immunohistochemistry of SPP1. (A) Sensitivity to sorafenib by risk group. (B) Response rate to ICB therapy. (C) Heatmap showing the correlation between risk score and the three model genes. (D) Representative immunohistochemically stained images of HCC tissues showing SPP1 expression (upper, 5.0×; lower, 20.0×). (E) K-M OS curves of 19 HCC patients with different SPP1 levels.

Discussion

Predicting the effect of HCC treatment is challenging due to the paucity of useful biomarkers. The American Joint Committee on Cancer (AJCC) TNM staging system has long been considered reliable for predicting the prognosis of patients with HCC. However, it is based on macroscopic information, which does not reflect the biological features and heterogeneity of HCC. Integration of prognostic gene signatures and traditional parameters has advantages for predicting the prognosis of HCC. The importance of predicting the prognosis of HCC and administering treatments in a timely manner highlights the need to identify robust prognostic biomarkers for risk stratification.

Few studies have focused on the predictive value of hypoxia- or angiogenesis-related signatures. The complex regulatory mechanisms and effects of hypoxia-related angiogenesis genes (HRAGs) give rise to an ambiguous relationship with the prognosis of HCC. To our knowledge, this study is the first to develop a novel scoring system based on HRAGs. We investigated the expression of HRAGs in HCC tumor tissues, and their associations with survival, and constructed a novel prognostic model based on three HRAGs (ANGPT2, SERPINE1, and SPP1) to stratify HCC patients according to their estimated survival. The prognostic model showed good predictive power in both the training cohort and two external validation cohorts. Univariate and multivariate Cox regression analyses revealed the independent prognostic value of the three genes comprising the gene signature. A nomogram showed that the gene signature was predictive of 1-, 2-, 3-, and 5-year OS, and so may be useful when planning short-term follow-up. However, only three databases were used in this study, and the signature needs to be validated in an independent cohort.

The prognostic gene signature comprised three HRAGs (ANGPT2, SERPINE1, and SPP1). ANGPT2, which belongs to the angiopoietin family, is highly expressed in diverse tumor cells, and is implicated in tumor angiogenesis and inflammation (20, 21). ANGPT2 is highly expressed in, and closely related to, the development and prognosis of, HCC (15, 22). Increased expression of plasminogen activator inhibitor type 1 (also known as SERPINE1) was associated with tumor cell migration and invasion via activation of the PI3K-Akt pathway (23, 24). High expression of SERPINE1 in cancer tissues predicts a poor clinical outcome. In this study, SERPINE1 expression was lower in HCC tumor tissues than adjacent normal tissues in the TCGA dataset (Figure 2), possibly because of the different roles SERPINE1 in tumor and normal tissues. SPP1, also called osteopontin, is an arginine-glycine-aspartate-containing phosphoprotein that is overexpressed in many cancers, including lung adenocarcinoma and HCC, and serves as a prognostic biomarker (25–27). SPP1 is involved in tumor immunosuppression and affects the TME (28). The expression of the three prognostic model genes, and the risk score, increased with increasing tumor stage, suggesting that they correlate with tumor malignancy and progression (Supplementary Figure S2).

To identify the pathway involved in hypoxia-related angiogenesis, HCC patients were divided into two groups according to the median risk score calculated by GSEA. Biological function analysis revealed greater activity not only in hypoxia- and angiogenesis-related pathways, but also in the EMT pathway (Figure 6). The EMT is a process of phenotypic plasticity that has roles in organ development, wound healing, tumor progression, and the response to therapeutics (29). Our findings suggest that HRAGs promote HCC progression by regulating the EMT pathway.

Because of its importance for immunotherapy and, potentially, precision therapy, the TME is a focus of research (30). In terms of predictive biomarkers, our high-risk group had higher proportions of macrophages, Tregs, and MDSCs (Figure 7). Higher proportions of tumor-associated macrophages and Tregs are associated with a poor prognosis of HCC (31). In addition, impairment of the type II IFN response and increased infiltration of MDSCs, Tregs, and macrophages, as in the high-risk group, is implicated in tumor immunological escape and tolerance, and impairs the antitumor T-cell response in HCC (32–34). Furthermore, the risk score was positively related to the proportion of CAFs in HCC. As a critical component of the TME, CAFs contribute to immune evasion and immunotherapy failure, and promote the proliferation and invasion of tumors, including HCC, by secreting growth factors and cytokines (35–37). The immune-related pathways IL2/STAT5 and IL6/JAK/STAT3, as well as the inflammatory response, IFN-γ response, and TNFα signaling via NFκB, were significantly enriched in the high-risk group (Figure 6). Cancer cells may drive the expression of immune checkpoints via these immune-related pathways (38). The expression of immune-related checkpoints (PD-L1, LAG3, CTLA4, TIM3, and TIGIT) in our high-risk group was higher than in the low-risk group, but the MSI and TMB, which indirectly reflect the ability of a tumor to produce new antigens and predict the efficacy of immunotherapy for a variety of tumors, were not significantly different between the groups (Supplementary Figure S5). The expression of the pro-angiogenic factors VEGFA and VEGFB was significantly upregulated in the high-risk group in the ICGC and TCGA databases; these cytokines promote TME remodeling and directly or indirectly regulate immune cells (17). Therefore, in the high-risk group, anti-tumor immunity is probably attenuated; also, a high risk score may be correlated with immunosuppression in HCC, possibly explaining its poor prognosis. This may arise because the expression of angiogenesis-related genes under hypoxic conditions induces the production of myeloid suppressor cells, Tregs, and immunosuppression-related checkpoints, thereby disrupting the immune balance. In summary, tumor-driven hypoxia related angiogenesis plays a crucial role in modulating the tumor immune microenvironment.

In view of the fact that angiogenesis plays a central role in immunosuppression and can lead to resistance to ICBs, there is growing evidence to support a strategy combining anti-angiogenesis and ICBs with promising clinical activity. Due to the encouraging efficacy and safety findings of the IMbrave150 trial for atezolizumab plus bevacizumab, this novel anti-angiogenesis combined with immunotherapy has become the first-line treatment for patients with unresectable HCC (39, 40). The combination therapy can render the TME conducive to immune cell function. However, because angiogenesis and the TME have multiple roles, the mechanism underlying the effect of anti-angiogenesis combined with immunotherapy on liver cancer is unclear.

Our findings suggest that HRAGs not only contribute to tumor growth, metastasis, progression, and poor prognosis, but also to the infiltration of immunosuppressive cells and high expression of immune checkpoints. This could explain why anti-angiogenesis combined with immunotherapy is effective for unresectable HCC. The anti-angiogenesis and ICB therapy response rate was higher in the low-risk group. And SPP1 expression was most relevant to the risk score. Previous studies have illustrated that tumor-driven hypoxia promotes the expression of SPP1, which in turn promotes tumor angiogenesis and immunosuppressive microenvironment (28, 41, 42). SPP1 may be considered as a general marker of cancer progression, would be valuable in combination with other biomarkers to guide patient stratification and treatment strategies, and would be an attractive therapeutic target due to its multiple roles in promoting tumor aggressiveness. Our present study interestingly found that the OS of the low-SPP1-expression group of HCC patients who received anti-angiogenesis combined with immunotherapy after resection was trend longer than that of the high-SPP1-expression group (Figure 8). These findings may aid the development of comprehensive therapeutic strategies for HCC.

This study also had some limitations. First, the gene signature model was constructed and validated based on retrospective data from public databases. Second, the small sample size hampered statistical analysis. Third, a further study should investigate the biological mechanisms underlying the signature.

Conclusion

A novel gene signature model and HRAGs-based prognostic biomarker were developed and validated in an independent cohort. The gene signature may help identify immune infiltration and predict sensitivity to anti-angiogenesis and immunotherapy, and the outcomes of HCC. The mechanism underlying the associations of HRAGs expression in HCC with the TME and sensitivity of anti-angiogenesis and immunotherapy are unclear, so further studies are needed.

Data availability statement

The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found below: TCGA-LIHC: https://portal.gdc.cancer.gov/; ICGC: https://dcc.icgc.org/projects/LIRI-JP; GEO: https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE14520.

Ethics statement

The studies involving human participants were reviewed and approved by the Institutional Review Boards of the Research Institute and Hospital National Cancer Center and The First Affiliated Hospital, Sun Yat-Sen University. The patients/participants provided their written informed consent to participate in this study.

Author contributions

JL, GZ and YX conceived and designed the project. JL supervised the study. GZ and YX conducted and performed data collection and data analysis and interpretation. XZ, WF, YZ, YW, and HW contributed to discussion and reviewed and edited the manuscript. All authors read and approved the final version of the manuscript.

Funding

This research was supported by the National Natural Science Foundation of China (NSFC) (82172036, 81971719 and 81171441), and the Key Research and Development Project of Guangzhou City (202103000021), and the major scientific and technological project of Guangdong Province (2020B0101130016), and the National Natural Science Foundation of China (Youth Project, 82102161).

Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed 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/fonc.2022.978050/full#supplementary-material

Abbreviations

AUC, Areas under the curve; CAFs, Cancer-associated fibroblasts; DEGs, Differentially expressed genes; DFS, disease free survival; FDG, false discovery rate; GEO, Gene-Expression Omnibus; GSEA, Gene Set Enrichment Analysis; HCC, hepatocellular carcinoma; HRAGs, Hypoxia-related angiogenesis genes; ICB, Immunological checkpoint blockade; ICGC, International Cancer Genome Consortium; MDSCs, Myeloid-derived suppressor cells; OS, overall survival; PCA, Principal component analysis; ROC, Receiver operator characteristic; ssGSEA, Single-sample gene-set enrichment analysis; TCGA, The Cancer Genome Atlas; TME, Tumor microenvironment; t-SNE, t-distributed stochastic neighborembedding.

References

1. Sung H, Ferlay J, Siegel RL, Laversanne M, Soerjomataram I, Jemal A, et al. Global cancer statistics 2020: Globocan estimates of incidence and mortality worldwide for 36 cancers in 185 countries. CA Cancer J Clin (2021) 71(3):209–49. doi: 10.3322/caac.21660

PubMed Abstract | CrossRef Full Text | Google Scholar

2. Allemani C, Matsuda T, Di Carlo V, Harewood R, Matz M, Niksic M, et al. Global surveillance of trends in cancer survival 2000-14 (Concord-3): Analysis of individual records for 37 513 025 patients diagnosed with one of 18 cancers from 322 population-based registries in 71 countries. Lancet (2018) 391(10125):1023–75. doi: 10.1016/S0140-6736(17)33326-3

PubMed Abstract | CrossRef Full Text | Google Scholar

3. Wang C, Vegna S, Jin H, Benedict B, Lieftink C, Ramirez C, et al. Inducing and exploiting vulnerabilities for the treatment of liver cancer. Nature (2019) 574(7777):268–72. doi: 10.1038/s41586-019-1607-3

PubMed Abstract | CrossRef Full Text | Google Scholar

4. Cancer Genome Atlas Research Network. Electronic address wbe, cancer genome atlas research n. comprehensive and integrative genomic characterization of hepatocellular carcinoma. Cell (2017) 169(7):1327–41.e23. doi: 10.1016/j.cell.2017.05.046

PubMed Abstract | CrossRef Full Text | Google Scholar

5. Lin DC, Mayakonda A, Dinh HQ, Huang P, Lin L, Liu X, et al. Genomic and epigenomic heterogeneity of hepatocellular carcinoma. Cancer Res (2017) 77(9):2255–65. doi: 10.1158/0008-5472.CAN-16-2822

PubMed Abstract | CrossRef Full Text | Google Scholar

6. Torrecilla S, Sia D, Harrington AN, Zhang Z, Cabellos L, Cornella H, et al. Trunk mutational events present minimal intra- and inter-tumoral heterogeneity in hepatocellular carcinoma. J Hepatol (2017) 67(6):1222–31. doi: 10.1016/j.jhep.2017.08.013

PubMed Abstract | CrossRef Full Text | Google Scholar

7. Brown JM, Wilson WR. Exploiting tumour hypoxia in cancer treatment. Nat Rev Cancer (2004) 4(6):437–47. doi: 10.1038/nrc1367

PubMed Abstract | CrossRef Full Text | Google Scholar

8. Jain RK. Normalization of tumor vasculature: An emerging concept in antiangiogenic therapy. Science (2005) 307(5706):58–62. doi: 10.1126/science.1104819

PubMed Abstract | CrossRef Full Text | Google Scholar

9. Wilson GK, Tennant DA, McKeating JA. Hypoxia inducible factors in liver disease and hepatocellular carcinoma: Current understanding and future directions. J Hepatol (2014) 61(6):1397–406. doi: 10.1016/j.jhep.2014.08.025

PubMed Abstract | CrossRef Full Text | Google Scholar

10. Lee SY, Jeong EK, Ju MK, Jeon HM, Kim MY, Kim CH, et al. Induction of metastasis, cancer stem cell phenotype, and oncogenic metabolism in cancer cells by ionizing radiation. Mol Cancer (2017) 16(1):10. doi: 10.1186/s12943-016-0577-4

PubMed Abstract | CrossRef Full Text | Google Scholar

11. Terry S, Buart S, Chouaib S. Hypoxic stress-induced tumor and immune plasticity, suppression, and impact on tumor heterogeneity. Front Immunol (2017) 8:1625. doi: 10.3389/fimmu.2017.01625

PubMed Abstract | CrossRef Full Text | Google Scholar

12. Zeng F, Zhang Y, Han X, Zeng M, Gao Y, Weng J. Employing hypoxia characterization to predict tumour immune microenvironment, treatment sensitivity and prognosis in hepatocellular carcinoma. Comput Struct Biotechnol J (2021) 19:2775–89. doi: 10.1016/j.csbj.2021.03.033

PubMed Abstract | CrossRef Full Text | Google Scholar

13. Palazon A, Goldrath AW, Nizet V, Johnson RS. Hif transcription factors, inflammation, and immunity. Immunity (2014) 41(4):518–28. doi: 10.1016/j.immuni.2014.09.008

PubMed Abstract | CrossRef Full Text | Google Scholar

14. Hilmi M, Neuzillet C, Calderaro J, Lafdil F, Pawlotsky JM, Rousseau B. Angiogenesis and immune checkpoint inhibitors as therapies for hepatocellular carcinoma: Current knowledge and future research directions. J Immunother Cancer (2019) 7(1):333. doi: 10.1186/s40425-019-0824-5

PubMed Abstract | CrossRef Full Text | Google Scholar

15. Morse MA, Sun W, Kim R, He AR, Abada PB, Mynderse M, et al. The role of angiogenesis in hepatocellular carcinoma. Clin Cancer Res Off J Am Assoc Cancer Res (2019) 25(3):912–20. doi: 10.1158/1078-0432.CCR-18-1254

CrossRef Full Text | Google Scholar

16. Rahma OE, Hodi FS. The intersection between tumor angiogenesis and immune suppression. Clin Cancer Res Off J Am Assoc Cancer Res (2019) 25(18):5449–57. doi: 10.1158/1078-0432.CCR-18-1543

CrossRef Full Text | Google Scholar

17. Tian Y, Li Y, Shao Y, Zhang Y. Gene modification strategies for next-generation car T cells against solid cancers. J Hematol Oncol (2020) 13(1):54. doi: 10.1186/s13045-020-00890-6

PubMed Abstract | CrossRef Full Text | Google Scholar

18. Zhen Z, Shen Z, Hu Y, Sun P. Screening and identification of angiogenesis-related genes as potential novel prognostic biomarkers of hepatocellular carcinoma through bioinformatics analysis. Aging (Albany NY) (2021) 13(13):17707–33. doi: 10.18632/aging.203260

PubMed Abstract | CrossRef Full Text | Google Scholar

19. Miao YR, Zhang Q, Lei Q, Luo M, Xie GY, Wang H, et al. Immucellai: A unique method for comprehensive T-cell subsets abundance prediction and its application in cancer immunotherapy. Adv Sci (Weinh) (2020) 7(7):1902880. doi: 10.1002/advs.201902880

PubMed Abstract | CrossRef Full Text | Google Scholar

20. Scholz A, Plate KH, Reiss Y. Angiopoietin-2: A multifaceted cytokine that functions in both angiogenesis and inflammation. Ann N Y Acad Sci (2015) 1347:45–51. doi: 10.1111/nyas.12726

PubMed Abstract | CrossRef Full Text | Google Scholar

21. Akwii RG, Sajib MS, Zahra FT, Mikelis CM. Role of angiopoietin-2 in vascular physiology and pathophysiology. Cells (2019) 8(5):471. doi: 10.3390/cells8050471

CrossRef Full Text | Google Scholar

22. Villa E, Critelli R, Lei B, Marzocchi G, Camma C, Giannelli G, et al. Neoangiogenesis-related genes are hallmarks of fast-growing hepatocellular carcinomas and worst survival. Results Prospective Study Gut (2016) 65(5):861–9. doi: 10.1136/gutjnl-2014-308483

CrossRef Full Text | Google Scholar

23. Balsara RD, Castellino FJ, Ploplis VA. A novel function of plasminogen activator inhibitor-1 in modulation of the akt pathway in wild-type and plasminogen activator inhibitor-1-Deficient endothelial cells. J Biol Chem (2006) 281(32):22527–36. doi: 10.1074/jbc.M512819200

PubMed Abstract | CrossRef Full Text | Google Scholar

24. Li R, Li Y, Liang X, Yang L, Su M, Lai KP. Network pharmacology and bioinformatics analyses identify intersection genes of niacin and covid-19 as potential therapeutic targets. Brief Bioinform (2021) 22(2):1279–90. doi: 10.1093/bib/bbaa300

PubMed Abstract | CrossRef Full Text | Google Scholar

25. Shin HD, Park BL, Cheong HS, Yoon JH, Kim YJ, Lee HS. Spp1 polymorphisms associated with hbv clearance and hcc occurrence. Int J Epidemiol (2007) 36(5):1001–8. doi: 10.1093/ije/dym093

PubMed Abstract | CrossRef Full Text | Google Scholar

26. Shen XY, Liu XP, Song CK, Wang YJ, Li S, Hu WD. Genome-wide analysis reveals alcohol dehydrogenase 1c and secreted phosphoprotein 1 for prognostic biomarkers in lung adenocarcinoma. J Cell Physiol (2019) 234(12):22311–20. doi: 10.1002/jcp.28797

PubMed Abstract | CrossRef Full Text | Google Scholar

27. Song Z, Chen W, Athavale D, Ge X, Desert R, Das S, et al. Osteopontin takes center stage in chronic liver disease. Hepatology (2021) 73(4):1594–608. doi: 10.1002/hep.31582

PubMed Abstract | CrossRef Full Text | Google Scholar

28. Shurin MR. Osteopontin controls immunosuppression in the tumor microenvironment. J Clin Invest (2018) 128(12):5209–12. doi: 10.1172/JCI124918

PubMed Abstract | CrossRef Full Text | Google Scholar

29. Dongre A, Weinberg RA. New insights into the mechanisms of epithelial-mesenchymal transition and implications for cancer. Nat Rev Mol Cell Biol (2019) 20(2):69–84. doi: 10.1038/s41580-018-0080-4

PubMed Abstract | CrossRef Full Text | Google Scholar

30. Binnewies M, Roberts EW, Kersten K, Chan V, Fearon DF, Merad M, et al. Understanding the tumor immune microenvironment (Time) for effective therapy. Nat Med (2018) 24(5):541–50. doi: 10.1038/s41591-018-0014-x

PubMed Abstract | CrossRef Full Text | Google Scholar

31. Zhou SL, Zhou ZJ, Hu ZQ, Huang XW, Wang Z, Chen EB, et al. Tumor-associated neutrophils recruit macrophages and T-regulatory cells to promote progression of hepatocellular carcinoma and resistance to sorafenib. Gastroenterology (2016) 150(7):1646–58.e17. doi: 10.1053/j.gastro.2016.02.040

PubMed Abstract | CrossRef Full Text | Google Scholar

32. Lin Z, Xu Q, Miao D, Yu F. An inflammatory response-related gene signature can impact the immune status and predict the prognosis of hepatocellular carcinoma. Front Oncol (2021) 11:644416. doi: 10.3389/fonc.2021.644416

PubMed Abstract | CrossRef Full Text | Google Scholar

33. Chiu DK, Tse AP, Xu IM, Di Cui J, Lai RK, Li LL, et al. Hypoxia inducible factor hif-1 promotes myeloid-derived suppressor cells accumulation through Entpd2/Cd39l1 in hepatocellular carcinoma. Nat Commun (2017) 8(1):517. doi: 10.1038/s41467-017-00530-7

PubMed Abstract | CrossRef Full Text | Google Scholar

34. Lin Y, Yang X, Liu W, Li B, Yin W, Shi Y, et al. Chemerin has a protective role in hepatocellular carcinoma by inhibiting the expression of il-6 and gm-csf and mdsc accumulation. Oncogene (2017) 36(25):3599–608. doi: 10.1038/onc.2016.516

PubMed Abstract | CrossRef Full Text | Google Scholar

35. Kubo N, Araki K, Kuwano H, Shirabe K. Cancer-associated fibroblasts in hepatocellular carcinoma. World J Gastroenterol (2016) 22(30):6841–50. doi: 10.3748/wjg.v22.i30.6841

PubMed Abstract | CrossRef Full Text | Google Scholar

36. Chakravarthy A, Khan L, Bensler NP, Bose P, De Carvalho DD. Tgf-Beta-Associated extracellular matrix genes link cancer-associated fibroblasts to immune evasion and immunotherapy failure. Nat Commun (2018) 9(1):4692. doi: 10.1038/s41467-018-06654-8

PubMed Abstract | CrossRef Full Text | Google Scholar

37. Liu T, Han C, Wang S, Fang P, Ma Z, Xu L, et al. Cancer-associated fibroblasts: An emerging target of anti-cancer immunotherapy. J Hematol Oncol (2019) 12(1):86. doi: 10.1186/s13045-019-0770-1

PubMed Abstract | CrossRef Full Text | Google Scholar

38. Garcia-Diaz A, Shin DS, Moreno BH, Saco J, Escuin-Ordinas H, Rodriguez GA, et al. Interferon receptor signaling pathways regulating pd-L1 and pd-L2 expression. Cell Rep (2019) 29(11):3766. doi: 10.1016/j.celrep.2019.11.113

PubMed Abstract | CrossRef Full Text | Google Scholar

39. Finn RS, Qin S, Ikeda M, Galle PR, Ducreux M, Kim TY, et al. Atezolizumab plus bevacizumab in unresectable hepatocellular carcinoma. N Engl J Med (2020) 382(20):1894–905. doi: 10.1056/NEJMoa1915745

PubMed Abstract | CrossRef Full Text | Google Scholar

40. Casak SJ, Donoghue M, Fashoyin-Aje L, Jiang X, Rodriguez L, Shen YL, et al. Fda approval summary: Atezolizumab plus bevacizumab for the treatment of patients with advanced unresectable or metastatic hepatocellular carcinoma. Clin Cancer Res Off J Am Assoc Cancer Res (2021) 27(7):1836–41. doi: 10.1158/1078-0432.CCR-20-3407

CrossRef Full Text | Google Scholar

41. Zhu Y, Yang J, Xu D, Gao XM, Zhang Z, Hsu JL, et al. Disruption of tumour-associated macrophage trafficking by the osteopontin-induced colony-stimulating factor-1 signalling sensitises hepatocellular carcinoma to anti-Pd-L1 blockade. Gut (2019) 68(9):1653–66. doi: 10.1136/gutjnl-2019-318419

PubMed Abstract | CrossRef Full Text | Google Scholar

42. Amilca-Seba K, Sabbah M, Larsen AK, Denis JA. Osteopontin as a regulator of colorectal cancer progression and its clinical applications. Cancers (2021) 13(15):3793. doi: 10.3390/cancers13153793

PubMed Abstract | CrossRef Full Text | Google Scholar