Vascular endothelial growth factor receptor-2 and its association with tumor immune regulatory gene expression in hepatocellular carcinoma

Ze-Long Liu; Ling-Ling Zhu; Jing-Hua Liu; Zhang-Ya Pu; Zhi-Ping Ruan; Jiang Chen

doi:10.18632/aging.104119

Research Paper Volume 12, Issue 24 pp 25172—25188

Vascular endothelial growth factor receptor-2 and its association with tumor immune regulatory gene expression in hepatocellular carcinoma

Ze-Long Liu^{1,5,
,} , Ling-Ling Zhu^{2,5,
,} , Jing-Hua Liu^{3,
*,} , Zhang-Ya Pu^5, , Zhi-Ping Ruan^5, , Jiang Chen^4,5, ,

¹ Division of Interventional Ultrasound, The First Affiliated Hospital of Sun Yat-sen University, Guangzhou, Guangdong Province, China
² Lung Cancer Center, West China Hospital of Sichuan University, Chengdu, Sichuan Province, China
³ Department of Hepatobiliary Surgery and Professor Cai’s Laboratory, Linyi People's Hospital, Linyi, Shandong Province, China
⁴ Department of General Surgery, Sir Run Run Shaw Hospital, Zhejiang University, Hangzhou, Zhejiang Province, China
⁵ Department of Radiation Oncology, Massachusetts General Hospital, Harvard Medical School, Boston, MA 02115, USA

* Equal contribution

Received: July 29, 2020 Accepted: September 9, 2020 Published: November 20, 2020

https://doi.org/10.18632/aging.104119
How to Cite

Copyright: © 2020 Liu et al. This is an open access article distributed under the terms of the Creative Commons Attribution License (CC BY 3.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.

Abstract

Anti-vascular endothelial growth factor (anti-VEGF) drugs have long been the only first-line treatment for advanced or unresectable hepatocellular carcinoma (HCC). Recently, the combination of bevacizumab (an anti-VEGF drug) and atezolizumab (an immune checkpoint blockade, ICB) has been proven to have superior efficacy over sorafenib. However, the complex association between VEGF signaling pathway and tumor immune microenvironment is still largely unknown. Here, we analyzed the RNA sequencing and clinical data of 365 HCC patients obtained from The Cancer Genome Atlas to investigate the potential correlation between VEGF signaling pathway and tumor immune microenvironment, including immune cell infiltration, 66 immune markers, genomic instability, and immune-related pathways. Our study revealed that VEGF signaling pathway score was positively correlated with immune cell infiltration and the expression profile of 66 immune markers. Enrichment analysis indicated that genes differentially expressed between two VEGF score subtypes were enriched in many immune-related Gene Ontology terms. Most importantly, both VEGF signaling pathway and activated CD8+ T cells were positively correlated with prognosis. Our findings suggest the co-activation of VEGF signaling pathway and tumor immune microenvironment in HCC patients, indicating the underlining mechanism of combination therapy including anti-VEGF drugs and ICBs.

Introduction

Hepatocellular carcinoma (HCC) is the fourth leading cause of cancer-related death worldwide [1]. Surgical resection has been adopted as the standard treatment for most patients with early-stage HCC [2, 3]. However, up to 70% of patients experience tumor recurrence within five years after surgery [2, 3]. Therefore, the biggest challenge in the clinical management of HCC is to prevent tumor recurrence and to improve the prognosis of advanced or unresectable disease.

Angiogenesis, the formation of new blood vessels, is crucial in tumor development, growth, and metastasis because blood vessels carry nutrients and facilitate the spread of tumor cells. Vascular endothelial growth factor (VEGF) signaling is important in angiogenesis and has been identified as a therapeutic target in diverse cancer types, including HCC [4]. Among the VEGF receptors, VEGF-receptor 2 (VEGFR2) is now considered the main receptor, and VEGFR2 activation promotes endothelial cell mitogenesis and vascular permeability [4]. Anti-VEGF monotherapy has been the only standard treatment for advanced or unresectable HCC over the past decade after sorafenib was approved by the US Food and Drug Administration (FDA) [5]. After sorafenib, lenvatinib was approved as a first-line treatment, while several other anti-VEGF drugs, including regorafenib, cabozantinib, and ramucirumab, were approved as second-line treatments [6–9]. Sorafenib provided longer overall survival (OS) compared with a placebo and lenvatinib was non-inferior to sorafenib [5, 6]. Despite this, patients with advanced HCC continue to have unsatisfactory clinical outcomes, with an OS of approximately one year [5–9]. To improve patient prognosis, there is an urgent medical need for novel treatments or improved clinical benefits of anti-VEGF drugs.

Over the last decade, cancer immunotherapy such as immune checkpoint blockade (ICB) therapy, has made tremendous breakthroughs in a variety of malignancies. However, only a small subset of HCC patients respond to ICB monotherapy [10, 11]. The major hurdles in achieving satisfactory benefits are high inter- and intra-tumor heterogeneity, various immune microenvironment, lack of predictive markers, and insufficient efficacy of monotherapy. The IMbrave 150 study was a phase 3 clinical trial that compared the efficacy of atezolizumab (ICB) plus bevacizumab (anti-VEGF drug) versus sorafenib as a first-line treatment for patients with unresectable HCC [12]. Recently, the results of this study demonstrated the impressive clinical benefits of the combination therapy of atezolizumab and bevacizumab, which resulted in a significantly prolonged OS [12]. Based on this, the US FDA has approved atezolizumab plus bevacizumab as a novel first-line treatment for advanced or unresectable HCC. In another phase 3 trial, LEAP-002, early clinical data has also indicated promising efficacy for the combination of lenvatinib and pembrolizumab (ICB) [13]. However, the underlying mechanisms of the combination of anti-VEGF drugs and ICBs have not yet been elucidated. Dual PD-1 and VEGFR2 blockade exhibited a synergistic antitumor effect in mouse models of HCC by promoting vascular normalization and enhancing antitumor immune responses [14]. Although the VEGF signaling pathway contributes to immunosuppressive tumor microenvironment, the association between VEGF signaling pathway and immune-related genes is poorly understood [15]. Recently, the expression profile of immune checkpoint genes was shown to be related to the prognosis of HCC patients [16]. Of note, NRP1, an immune checkpoint gene that is also involved in VEGF-VEGFR2 signaling, was found to be an important prognostic factor in HCC [16].

Therefore, we aimed to investigate the association between VEGF signaling pathway, tumor immune microenvironment, and the prognosis of HCC patients using data obtained from The Cancer Genome Atlas (TCGA; n = 365) in this study. It is important to gain a more comprehensive understanding of the field to improve the clinical management of HCC patients, especially for decision making regarding anti-VEGF treatments, ICB therapies, and their combination in patients with advanced HCC.

Results

Immune cell infiltration patterns associated with VEGF score subtypes

RNA sequencing and clinical data of 371 HCC patients were downloaded from the TCGA database, among which six were excluded due to incomplete follow-up information. Based on calculated enrichment scores of the VEGF signaling pathway, patients were divided into two VEGF subtypes: the high score subtype with scores in the top one-third (n = 122) and the low score subtype with scores in the bottom two-thirds (n = 243).

Based on the published data from the TCGA group, we compared the tumor purity and the stromal fraction between the high and the low VEGF score subtype. This revealed a significantly lower tumor purity and higher stromal fraction in the high VEGF score subtype (P < 0.0001) (Figure 1A). More importantly, the high VEGF score subtype exhibited a significantly higher leukocyte fraction (P < 0.0001) (Figure 1A). Based on the ESTIMATE results, we found that the high VEGF score subtype had a significantly higher immune score, stromal score, and Estimate score (all with P < 0.0001) (Figure 1B). Infiltration of 28 immune cell types was estimated for each by using ssGSEA scores. scores. Compared with the low VEGF score subtype, the high VEGF score subtype had relatively higher immune cell infiltration, including cells contributing to anti-tumor immunity (e.g., activated CD8+ T cells, type 1 T helper cells, and natural killer cells) and cells contributing to pro-tumor suppression (e.g., regulatory T cells, immature dendritic cells, and myeloid-derived suppressor cells, MDSCs) (Figure 1C). For each VEGF score subtype, Spearman’s correlation test suggested a positive association between these two categories of immune cells in the local tumor microenvironment (Figure 1D). The high VEGF score group exhibited both higher anti-tumor immunity and pro-tumor suppression (Figure 1E).

$Immune cell infiltration patterns of the two VEGF score subtypes. The high VEGF score subtype showed lower tumor purity, higher stroma fraction, and higher leukocyte fraction (A). Higher immune score, stromal score, and Estimate score based on the ESTIMATE algorithm were observed in the high VEGF score subtype (B). The high VEGF score subtype had relatively higher immune cell infiltration, including cells contributing to both anti-tumor immunity and pro-tumor suppression (C). A positive association between these two categories of immune cells in the local tumor microenvironment was observed in both VEGF score subtypes (D). The high VEGF score subtype showed both higher anti-tumor immunity and pro-tumor suppression (E).$

Figure 1. Immune cell infiltration patterns of the two VEGF score subtypes. The high VEGF score subtype showed lower tumor purity, higher stroma fraction, and higher leukocyte fraction (A). Higher immune score, stromal score, and Estimate score based on the ESTIMATE algorithm were observed in the high VEGF score subtype (B). The high VEGF score subtype had relatively higher immune cell infiltration, including cells contributing to both anti-tumor immunity and pro-tumor suppression (C). A positive association between these two categories of immune cells in the local tumor microenvironment was observed in both VEGF score subtypes (D). The high VEGF score subtype showed both higher anti-tumor immunity and pro-tumor suppression (E).

For specific immune cell types, the high VEGF score subtype had significantly higher infiltration of 26 immune cell types (except for memory B cells and type17 T helper cells) (Figure 2). Spearman’s correlation tests also indicated positive correlations between VEGF signaling pathway ssGSEA scores and the ssGSEA scores of 28 immune cell types (Supplementary Figure 1). Although the absolute levels of immune cells varied between the two VEGF score subtypes, the relative cell fractions of most immune cells estimated by the CIBERSORT algorithm appeared to be similar (Supplementary Figure 2A, 2B). Notably, higher fractions of M1 macrophages, M0 macrophages, and gamma delta T cells were observed in the high VEGF score subtype (Supplementary Figure 2A, 2C).

Comparison of the infiltration of 28 immune cell types between the two VEGF score subtypes. The high VEGF score subtype exhibited significantly higher infiltration of 26 immune cell types (except for memory B cells and type17 T helper cells).

Figure 2. Comparison of the infiltration of 28 immune cell types between the two VEGF score subtypes. The high VEGF score subtype exhibited significantly higher infiltration of 26 immune cell types (except for memory B cells and type17 T helper cells).

Immune-regulatory gene expression profiles associated with VEGF score subtypes

The high VEGF score subtype exhibited higher expression of 66 immune markers associated with immune stimulation or suppression (Figure 3). The correlations between the 66 immune marker expression scores in all patients, patients of the high VEGF score subtype, and patients of the low VEGF score subtypes were investigated by Spearman’s correlation tests and visualized in Supplementary Figures 3–5, respectively. We then focused on several important checkpoint molecules including PD-1, PD-L1, CTLA-4, IDO1, and LAG3, and found that the expression levels of these molecules were all significantly higher in the high VEGF score subtype (Figure 4A–4E).

Figure 3. The high VEGF score subtype showed a higher expression of 66 immune markers associated with immune stimulation or suppression.

Expression levels of important checkpoint molecules including PD-1 (A), PD-L1 (B), CTLA-4 (C), IDO1 (D), and LAG3 (E) were all significantly higher in the high VEGF score subtype. Significantly lower TMB (F), SNV neoantigen number (G), and indel neoantigen number (H) were observed in the high VEGF score subtype.

Figure 4. Expression levels of important checkpoint molecules including PD-1 (A), PD-L1 (B), CTLA-4 (C), IDO1 (D), and LAG3 (E) were all significantly higher in the high VEGF score subtype. Significantly lower TMB (F), SNV neoantigen number (G), and indel neoantigen number (H) were observed in the high VEGF score subtype.

To investigate the genomic changes associated with each VEGF score subtype, we compared the TMB, single-nucleotide variant (SNV) neoantigen number, and indel neoantigen number between the two VEGF score subtypes. We found significantly lower TMB, SNV neoantigen number, and indel neoantigen number in the high VEGF score subtype than in the low VEGF score subtype (Figure 4F–4H), which indicated relatively lower genome instability in the high VEGF score subtype.

DEGs between VEGF score subtypes were enriched in immune-related GO terms

In all, 495 DEGs between the two VEGF score subtypes were identified, among which 456 were upregulated, and 39 were downregulated in the high VEGF score subtype. GO pathway enrichment analysis revealed the following top GO terms: adaptive immune response, complement activation (classical pathway), humoral immune response mediated by circulating immunoglobulin in biological process (Figure 5A); immunoglobulin complex, immunoglobulin complex (circulating), external side of plasma membrane in cellular components (Figure 5B); antigen binding, immunoglobulin receptor binding, extracellular matrix structural constituent in molecular functions (Figure 5C). Moreover, we constructed a network of the top 20 GO summary terms via the Metascape website, which revealed the relationships between immune-related GO terms, angiogenesis-related GO terms, and stroma related GO terms (Figure 5D).

DEGs between the two VEGF score subtypes enriched in immune-related GO terms. GO pathway enrichment analysis revealed that immune-related GO terms ranked top in biological process (A), cellular components (B), and molecular functions (C). The network of the top 20 GO summary terms revealed the relationships between immune-related GO terms, angiogenesis-related GO terms, and stroma related GO terms (D).

Figure 5. DEGs between the two VEGF score subtypes enriched in immune-related GO terms. GO pathway enrichment analysis revealed that immune-related GO terms ranked top in biological process (A), cellular components (B), and molecular functions (C). The network of the top 20 GO summary terms revealed the relationships between immune-related GO terms, angiogenesis-related GO terms, and stroma related GO terms (D).

Among the 495 DGEs, 172 were immune-related genes according to the data obtained from ImmPort. To further investigate these immune-related genes, a PPI network consisting of 44 nodes and 112 edges was constructed using the STRING database and Cytoscape v3.6.1 (Figure 6A). CCL19 exhibited the highest connectivity degree (15) in the PPI network. Analysis via the MCODE algorithm revealed two important modules with MCODE scores of 8.444 (10 nodes and 38 edges) and 4 (4 nodes and 6 edges) (Figure 6B, 6C). These two top modules included 10 genes (including CCL11, CCL19, CCL21, CCL22, CCR4, CCR8, CXCL5, CXCL14, MMP9, and SSTR5) and four genes (including CD1E, CD19, CD79A, and CR2), respectively (Figure 6B, 6C).

PPI network of immune-related genes. The PPI network consisted of 44 nodes and 112 edges (A). Two important modules were observed via the MCODE algorithm (B, C).

Figure 6. PPI network of immune-related genes. The PPI network consisted of 44 nodes and 112 edges (A). Two important modules were observed via the MCODE algorithm (B, C).

High VEGF score subtype and high levels of activated CD8+ T cells predict favorable prognosis

Given the prognostic value of CD8+ T cells in a variety of solid tumors, we then divided the 365 HCC patients into a high activated CD8+ T cell group (CD8Hi, n = 122) and a low activated CD8+ T cell group (CD8Lo, n = 243) based on ssGSEA scores. Compared with the low VEGF score subtype, the high VEGF score subtype had a significantly longer recurrence free survival (RFS, P = 0.0191) and the tendency of a longer OS (P = 0.0685) (Figure 7A, 7B). We also observed significantly better RFS and OS in the CD8Hi group when compared with the CD8Lo group (P < 0.0001 and P = 0.0034, respectively) (Figure 7C, 7D). We then categorized the 365 patients into four groups (VEGFHiCD8Hi, VEGFHiCD8Lo, VEGFLoCD8Hi, and VEGFLoCD8Lo) based on the ssGSEA scores of the VEGF signaling pathway and activated CD8+ T cells. Upon comparing RFS and OS between the groups, both suggested better prognosis in the VEGFHiCD8Hi group compared with the VEGFHiCD8Lo group and the VEGFLoCD8Lo group, but not the VEGFLoCD8Hi group (Figure 7E, 7F).

Survival analysis of 365 HCC patients. High VEGF score subtype (A, B) and high activated CD8+ T cells (C, D) predicted favorable prognosis. Improved RFS and OS were observed in the VEGFHiCD8Hi group compared to the VEGFHiCD8Lo group and the VEGFLoCD8Lo group, but not the VEGFLoCD8Hi group (E, F).

Figure 7. Survival analysis of 365 HCC patients. High VEGF score subtype (A, B) and high activated CD8+ T cells (C, D) predicted favorable prognosis. Improved RFS and OS were observed in the VEGFHiCD8Hi group compared to the VEGFHiCD8Lo group and the VEGFLoCD8Lo group, but not the VEGFLoCD8Hi group (E, F).

Discussion

Cancer immunotherapies, particularly ICBs targeting the PD-1/PD-L1 and CTLA-4 immune checkpoints, have made great success in a variety of malignancies over the past decade and have revolutionized treatment strategies for many cancer patients. However, for HCC, the response rate and prognostic benefit of ICBs alone appear to be extremely limited [10, 11]. Recent evidence has impressively demonstrated that combination therapy of anti-VEGF drugs and ICBs provides promising efficacy in HCC patients [12, 13]. The success of such combination therapy requires a better understanding of the complex associations between the VEGF signaling pathway and the tumor immune microenvironment.

Pathologic angiogenesis plays a critical role in tumor development and progression. There is now much evidence demonstrating that the VEGF/VEGFR2 signaling pathway is one of the most important pathways regulating angiogenesis. Activation of the VEGF signaling pathway promotes endothelial cell proliferation and vascular permeability, which can facilitate immune cell infiltration [8]. Additionally, by altering the adhesion molecules expressed on immune cells (e.g., integrin ligands intercellular adhesion molecule 1) and endothelial cells (e.g., vascular cell adhesion protein 1), VEGF can control the trafficking of immune cells to tumors [15]. However, VEGF has also been demonstrated to induce immunosuppression through several mechanisms, including inhibiting the trafficking, proliferation, and effector function of cytotoxic T cells and the maturation and antigen presentation of dendritic cells, promoting the recruitment and proliferation of immunosuppressive cells, and leading to a hypoxic and low-pH immunosuppressive microenvironment [15]. In our study, using RNA sequencing data from TCGA database, we found that there were more infiltrated immune cells, including both immune-active cells and immune-suppressive cells, in the high VEGF score subtype, which supports previous theories. Moreover, positive correlations between VEGF scores and immune cell scores were observed for most of the 28 immune cell types. These findings implied that both immune-active cells and immune-suppressive cells might enter the tumor microenvironment when VEGF-related molecules are upregulated. This also indicates the potential for feedback from the recruitment and infiltration of immune-active cells to facilitate immune suppression. Notably, we found a lower fraction of M1 macrophages and a higher fraction of M2 macrophages in the low VEGF score subtype. High infiltration of M2 macrophages in the tumor stroma has been reported to be associated with a down-regulating antitumor immune response [30]. In contrast, the high VEGF score subtype exhibited a higher fraction of M1 macrophages, indicating a potentially more responsive tumor status to ICBs. This result was inconsistent with the previous findings that VEGF could contribute to macrophage polarization from M1 to M2 [15]. Given that the total amount of macrophages was also higher in the high VEGF score subtype and the spatial organization of macrophages was unclear, the relationship between VEGF and macrophages requires further investigation. Additionally, the fraction of gamma delta T cells was relatively higher in the high VEGF score subtype. Gamma delta T cells have been found to possess cytotoxic antitumor activity and insufficient levels of functional gamma delta T cells is associated with tumor recurrence and poor prognosis [31–33]. Given that the role of gamma delta T cells in HCC is poorly understood, further studies are needed to validate our finding.

Our data revealed that the expression levels of 66 immune markers, both immune stimulation-related genes and immune suppression-related genes, were higher in the high VEGF score subtype. Notably, immune checkpoint genes, including PD-1, PD-L1, CTLA-4, IDO1, and LAG3, also exhibited significantly higher expression levels in the high VEGF score subtype. This indicated that the tumor was more likely to be responsive to ICB treatments in high VEGF score subtype. Additionally, the high VEGF score subtype had both lower TMB and neoantigen numbers compared to the low VEGF score subtype, indicating lower genetic instability. High immune infiltration is thought to lower genetic instability by eliminating tumor subclones and inhibiting tumor evolution [34]. On the contrary, insufficient immune infiltration might allow immune escape, leading to tumor evolution and increased genomic instability [34]. Together, these findings suggest that the VEGF high score subtype is more immune inflamed and may be more likely to benefit from ICB treatments. In agreement with previous findings, numerous DEGs between the two VEGF score subtypes were immune-related genes (172 of 495). To comprehensively explore the differences in gene expression profiles, we conducted GO pathway enrichment analysis, which indicated that many immune-related GO terms were ranked top in biological process, cellular components, and molecular functions. Other top GO terms were primarily associated with stromal components (e.g., extracellular matrix, extracellular matrix structural constituent), which was also consistent with the high stromal score in the high VEGF score subtype. The network of the top 20 summary GO terms further indicated the interplay between the immune microenvironment, the development of blood vessels, and the stroma. Additional PPI analysis focused on the immune-related genes revealed two important modules that included 10 genes (including CCL11, CCL19, CCL21, CCL22, CCR4, CCR8, CXCL5, CXCL14, MMP9, and SSTR5) and four genes (including CD1E, CD19, CD79A, and CR2), respectively. The cytokines in the first module primarily display chemotactic activities for miscellaneous immune cells, which explain the higher immune infiltration in the high VEGF score subtype. Genes in the second module were associated with encoding important components of B cells, which are important in humoral immunity.

Although both immune active factors and immune suppressive factors were present in the local tumor microenvironment, the former appeared to be the dominating components. Our survival analysis also supported the findings mentioned above and indicated a better prognosis (especially RFS) for the high VEGF score subtype. Another possible explanation is that HCC patients of the high VEGF score subtype might respond better to widely used anti-VEGF therapies. Given the prognostic value of CD8+ T cells in various malignancies, we next focused on them. Activated CD8+ T cells, in agreement with previous literature, predicted favorable clinical outcomes in HCC patients [35–38]. We also observed that the survival curves of the VEGFLoCD8Lo group and the VEGFHiCD8Lo group were similar and lower than the other two groups (the VEGFHiCD8Hi group and the VEGFLoCD8Higroup), suggesting an important role for immune infiltration accompanying angiogenesis in eliminating the tumor.

One of the biggest limitations of our study is that there is no information regarding the spatial organization of immune cells in the tumors. As previously demonstrated, center tumor immune cells, but not invasive margin immune cells, primarily contribute to tumor elimination and favorable prognosis [39, 40]. Additionally, although we considered the whole VEGF/VEGFR2 pathway rather than single molecules to better represent angiogenesis, vessel development in tumors is indeed very complicated and could be affected by many other factors. Finally, we demonstrated the correlation between VEGF signaling and tumor immune microenvironment, but not a causal relationship. Given that our findings are based on bioinformatics analysis, further laboratory experiments are needed to validate these results.

Conclusions

In conclusion, our study found the potential co-activation of VEGF signaling pathway and tumor immune microenvironment, indicating a potential benefit for combination therapy including anti-VEGF and ICBs in HCC patients.

Materials and Methods

Data sources

The level 3 RNA sequencing and clinical data of HCC patients were retrieved from TCGA data portal (https://portal.gdc.cancer.gov).

Gene signatures and single-sample gene set enrichment analysis (ssGSEA) scores

The enrichment score of the VEGF signaling pathway (the KEGG_VEGF_SIGNALING_PATHWAY, Supplementary Figure 6), was calculated using ssGSEA. The KEGG_VEGF_SIGNALING_PATHWAY only focused on the VEGF/VEGFR2 axis. The gene set representing the KEGG_VEGF_SIGNALING_PATHWAY was downloaded from the GSEA Molecular Signatures Database (MSigDB) v7.1 (https://www.gsea-msigdb.org/gsea/downloads.jsp) and includes a total of 76 genes [17, 18]. Data on tumor purity, stromal fraction, and leukocyte fraction were obtained from a previously published study from the TCGA group [19]. For the quantification of infiltrated immune cells, we utilized a previously published calculation method based on the expression profiles of 782 genes from 28 immune cell types [20]. The degree of immune infiltration was estimated by the ssGSEA scores through the Gene Set Variation Analysis (GSVA) package in R v3.6.2 [21]. We obtained immune scores, stromal scores, and Estimate scores, which were calculated by using the ESTIMATE algorithm, for TCGA cohorts from the official data portal (https://bioinformatics.mdanderson.org/estimate/) [22]. The CIBERSORT algorithm (https://cibersort.stanford.edu) was used to determine the relative abundances of 22 immune cell types in the tumor microenvironment [23].

Immune signature and immune markers

To investigate the expression of immune markers, we calculated expression scores using a panel of 66 immune markers that are associated with immune stimulation or suppression in the tumor microenvironment [24]. To determine if VEGF signaling pathway subtypes are associated with genomic changes, we downloaded tumor mutation burden (TMB) and neoantigen data from two previously published studies from the TCGA group [19, 25].

Identification of differentially expressed genes (DEGs) and Gene Ontology (GO) pathway enrichment analyses

RNA sequencing data was analyzed using the limma package in R v3.6.2 to identify DEGs with cutoff values of |log2FC| ≥ 1.5 and FDR < 0.05 [26]. GO pathway enrichment analyses of DEGs were performed using the Metascape website (http://metascape.org) [27].

Functional enrichment of protein-protein interaction (PPI) network and module analysis of immune-related DEGs

Immune-related genes were extracted from the identified DEGs, based on the immune-related gene list downloaded from the Immunology Database and Analysis Portal (ImmPort) (https://immport.niaid.nih.gov). We utilized the STRING database to construct and download PPIs among the immune-related DEGs [28]. Cytoscape v3.6.1 was used to reconstruct the PPI network and calculate the connectivity degree of each node [29]. Hub clusters in the network were detected and visualized using the molecular complex detection (MCODE) algorithm.

Statistical analysis

Continuous parameters were compared using the Student’s t test. The correlations between the VEGF signaling pathway ssGSEA scores and the ssGSEA scores across 28 immune cell types were determined by Spearman’s correlation test. The correlation between the immune marker expression scores was determined by Spearman’s correlation test and was visualized using the corrplot package in R v3.6.2. To visualize and compare immune cell infiltration patterns and immune signatures across different VEGF signaling pathway subtypes, we generated heatmaps using the pheatmap package in R v3.6.2. Bubble plots were plotted using the ggplot2 package in R v3.6.2. For survival analysis, Kaplan–Meier curves were plotted and compared using the Cox proportional-hazards regression model. In all analyses, a P value of a two-sided test less than 0.05 was considered to be statistically significant. All statistical analyses were performed with GraphPad Prism v8.3.0 and R v3.6.2.

Supplementary Materials

Supplementary Figures

Author Contributions

ZLL, JHL, LLZ, ZYP, and JC contributed to conception and design of the work. ZLL, JHL, LLZ, and ZYP contributed to the acquisition, analysis, and interpretation of the data for the work. ZLL, JHL, and LLZ contributed to manuscript writing. JC contributed to manuscript review and revision.

Acknowledgments

We thank the patients and investigators who participated in TCGA for providing data.

Conflicts of Interest

The authors declare no potential conflicts of interest.

Funding

This work was supported by the National Natural Science Foundation of China (grant number: 81702981), the China Postdoctoral Science Foundation Funded Project (grant number: 2020T130584), and the China Scholarship Council (grant number: 201806325017).

References

1. Bray F, Ferlay J, Soerjomataram I, Siegel RL, Torre LA, Jemal A. Global cancer statistics 2018: GLOBOCAN estimates of incidence and mortality worldwide for 36 cancers in 185 countries. CA Cancer J Clin. 2018; 68:394–424. https://doi.org/10.3322/caac.21492 [PubMed]
2. European Association for the Study of the Liver. EASL clinical practice guidelines: management of hepatocellular carcinoma. J Hepatol. 2018; 69:182–236. https://doi.org/10.1016/j.jhep.2018.03.019 [PubMed]
3. Pinyol R, Montal R, Bassaganyas L, Sia D, Takayama T, Chau GY, Mazzaferro V, Roayaie S, Lee HC, Kokudo N, Zhang Z, Torrecilla S, Moeini A, et al. Molecular predictors of prevention of recurrence in HCC with sorafenib as adjuvant treatment and prognostic factors in the phase 3 STORM trial. Gut. 2019; 68:1065–75. https://doi.org/10.1136/gutjnl-2018-316408 [PubMed]
4. Apte RS, Chen DS, Ferrara N. VEGF in signaling and disease: beyond discovery and development. Cell. 2019; 176:1248–64. https://doi.org/10.1016/j.cell.2019.01.021 [PubMed]
5. Llovet JM, Ricci S, Mazzaferro V, Hilgard P, Gane E, Blanc JF, de Oliveira AC, Santoro A, Raoul JL, Forner A, Schwartz M, Porta C, Zeuzem S, et al, and SHARP Investigators Study Group. Sorafenib in advanced hepatocellular carcinoma. N Engl J Med. 2008; 359:378–90. https://doi.org/10.1056/NEJMoa0708857 [PubMed]
6. Kudo M, Finn RS, Qin S, Han KH, Ikeda K, Piscaglia F, Baron A, Park JW, Han G, Jassem J, Blanc JF, Vogel A, Komov D, et al. Lenvatinib versus sorafenib in first-line treatment of patients with unresectable hepatocellular carcinoma: a randomised phase 3 non-inferiority trial. Lancet. 2018; 391:1163–73. https://doi.org/10.1016/S0140-6736(18)30207-1 [PubMed]
7. Bruix J, Qin S, Merle P, Granito A, Huang YH, Bodoky G, Pracht M, Yokosuka O, Rosmorduc O, Breder V, Gerolami R, Masi G, Ross PJ, et al, and RESORCE Investigators. Regorafenib for patients with hepatocellular carcinoma who progressed on sorafenib treatment (RESORCE): a randomised, double-blind, placebo-controlled, phase 3 trial. Lancet. 2017; 389:56–66. https://doi.org/10.1016/S0140-6736(16)32453-9 [PubMed]
8. Abou-Alfa GK, Meyer T, Cheng AL, El-Khoueiry AB, Rimassa L, Ryoo BY, Cicin I, Merle P, Chen Y, Park JW, Blanc JF, Bolondi L, Klümpen HJ, et al. Cabozantinib in patients with advanced and progressing hepatocellular carcinoma. N Engl J Med. 2018; 379:54–63. https://doi.org/10.1056/NEJMoa1717002 [PubMed]
9. Zhu AX, Kang YK, Yen CJ, Finn RS, Galle PR, Llovet JM, Assenat E, Brandi G, Pracht M, Lim HY, Rau KM, Motomura K, Ohno I, et al, and REACH-2 study investigators. Ramucirumab after sorafenib in patients with advanced hepatocellular carcinoma and increased α-fetoprotein concentrations (REACH-2): a randomised, double-blind, placebo-controlled, phase 3 trial. Lancet Oncol. 2019; 20:282–96. https://doi.org/10.1016/S1470-2045(18)30937-9 [PubMed]
10. El-Khoueiry AB, Sangro B, Yau T, Crocenzi TS, Kudo M, Hsu C, Kim TY, Choo SP, Trojan J, Welling TH 3rd, Meyer T, Kang YK, Yeo W, et al. Nivolumab in patients with advanced hepatocellular carcinoma (CheckMate 040): an open-label, non-comparative, phase 1/2 dose escalation and expansion trial. Lancet. 2017; 389:2492–502. https://doi.org/10.1016/S0140-6736(17)31046-2 [PubMed]
11. Zhu AX, Finn RS, Edeline J, Cattan S, Ogasawara S, Palmer D, Verslype C, Zagonel V, Fartoux L, Vogel A, Sarker D, Verset G, Chan SL, et al, and KEYNOTE-224 investigators. Pembrolizumab in patients with advanced hepatocellular carcinoma previously treated with sorafenib (KEYNOTE-224): a non-randomised, open-label phase 2 trial. Lancet Oncol. 2018; 19:940–52. https://doi.org/10.1016/S1470-2045(18)30351-6 [PubMed]
12. Cheng AL, Qin S, Ikeda M, Galle P, Ducreux M, Zhu A, Kim TY, Kudo M, Breder V, Merle P, Kaseb A, Li D, Verret W, et al. LBA3IMbrave150: Efficacy and safety results from a ph III study evaluating atezolizumab (atezo) + bevacizumab (bev) vs sorafenib (Sor) as first treatment (tx) for patients (pts) with unresectable hepatocellular carcinoma (HCC). Annals of Oncology. 2019; 30:ix186–87. https://doi.org/10.1093/annonc/mdz446.002
13. Llovet JM, Kudo M, Cheng AL, Finn R, Galle P, Kaneko S, Meyer T, Qin S, Dutcus C, Chen E, Dubrovsky L, Zhu A. Lenvatinib (len) plus pembrolizumab (pembro) for the first-line treatment of patients (pts) with advanced hepatocellular carcinoma (HCC): Phase 3 LEAP-002 study. Journal of Clinical Oncology. 2019; 37: TPS4152-TPS. https://doi.org/10.1200/JCO.2019.37.15_suppl.TPS4152
14. Shigeta K, Datta M, Hato T, Kitahara S, Chen IX, Matsui A, Kikuchi H, Mamessier E, Aoki S, Ramjiawan RR, Ochiai H, Bardeesy N, Huang P, et al. Dual programmed death receptor-1 and vascular endothelial growth factor receptor-2 blockade promotes vascular normalization and enhances antitumor immune responses in hepatocellular carcinoma. Hepatology. 2020; 71:1247–61. https://doi.org/10.1002/hep.30889 [PubMed]
15. Fukumura D, Kloepper J, Amoozgar Z, Duda DG, Jain RK. Enhancing cancer immunotherapy using antiangiogenics: opportunities and challenges. Nat Rev Clin Oncol. 2018; 15:325–40. https://doi.org/10.1038/nrclinonc.2018.29 [PubMed]
16. Xu D, Liu X, Wang Y, Zhou K, Wu J, Chen JC, Chen C, Chen L, Zheng J. Identification of immune subtypes and prognosis of hepatocellular carcinoma based on immune checkpoint gene expression profile. Biomed Pharmacother. 2020; 126:109903. https://doi.org/10.1016/j.biopha.2020.109903 [PubMed]
17. Subramanian A, Tamayo P, Mootha VK, Mukherjee S, Ebert BL, Gillette MA, Paulovich A, Pomeroy SL, Golub TR, Lander ES, Mesirov JP. Gene set enrichment analysis: a knowledge-based approach for interpreting genome-wide expression profiles. Proc Natl Acad Sci USA. 2005; 102:15545–50. https://doi.org/10.1073/pnas.0506580102 [PubMed]
18. Mootha VK, Lindgren CM, Eriksson KF, Subramanian A, Sihag S, Lehar J, Puigserver P, Carlsson E, Ridderstråle M, Laurila E, Houstis N, Daly MJ, Patterson N, et al. PGC-1alpha-responsive genes involved in oxidative phosphorylation are coordinately downregulated in human diabetes. Nat Genet. 2003; 34:267–73. https://doi.org/10.1038/ng1180 [PubMed]
19. Thorsson V, Gibbs DL, Brown SD, Wolf D, Bortone DS, Ou Yang TH, Porta-Pardo E, Gao GF, Plaisier CL, Eddy JA, Ziv E, Culhane AC, Paull EO, et al, and Cancer Genome Atlas Research Network. The immune landscape of cancer. Immunity. 2018; 48:812–30.e14. https://doi.org/10.1016/j.immuni.2018.03.023 [PubMed]
20. Charoentong P, Finotello F, Angelova M, Mayer C, Efremova M, Rieder D, Hackl H, Trajanoski Z. Pan-cancer immunogenomic analyses reveal genotype-immunophenotype relationships and predictors of response to checkpoint blockade. Cell Rep. 2017; 18:248–62. https://doi.org/10.1016/j.celrep.2016.12.019 [PubMed]
21. Hänzelmann S, Castelo R, Guinney J. GSVA: gene set variation analysis for microarray and RNA-seq data. BMC Bioinformatics. 2013; 14:7. https://doi.org/10.1186/1471-2105-14-7 [PubMed]
22. Yoshihara K, Shahmoradgoli M, Martínez E, Vegesna R, Kim H, Torres-Garcia W, Treviño V, Shen H, Laird PW, Levine DA, Carter SL, Getz G, Stemke-Hale K, et al. Inferring tumour purity and stromal and immune cell admixture from expression data. Nat Commun. 2013; 4:2612. https://doi.org/10.1038/ncomms3612 [PubMed]
23. Chen B, Khodadoust MS, Liu CL, Newman AM, Alizadeh AA. Profiling tumor infiltrating immune cells with CIBERSORT. Methods Mol Biol. 2018; 1711:243–59. https://doi.org/10.1007/978-1-4939-7493-1_12 [PubMed]
24. Cancer Genome Atlas Research Network. Comprehensive and integrative genomic characterization of hepatocellular carcinoma. Cell. 2017; 169:1327–41.e23. https://doi.org/10.1016/j.cell.2017.05.046 [PubMed]
25. Hoadley KA, Yau C, Hinoue T, Wolf DM, Lazar AJ, Drill E, Shen R, Taylor AM, Cherniack AD, Thorsson V, Akbani R, Bowlby R, Wong CK, et al, and Cancer Genome Atlas Network. Cell-of-origin patterns dominate the molecular classification of 10,000 tumors from 33 types of cancer. Cell. 2018; 173:291–304.e6. https://doi.org/10.1016/j.cell.2018.03.022 [PubMed]
26. Ritchie ME, Phipson B, Wu D, Hu Y, Law CW, Shi W, Smyth GK. Limma powers differential expression analyses for RNA-sequencing and microarray studies. Nucleic Acids Res. 2015; 43:e47. https://doi.org/10.1093/nar/gkv007 [PubMed]
27. Zhou Y, Zhou B, Pache L, Chang M, Khodabakhshi AH, Tanaseichuk O, Benner C, Chanda SK. Metascape provides a biologist-oriented resource for the analysis of systems-level datasets. Nat Commun. 2019; 10:1523. https://doi.org/10.1038/s41467-019-09234-6 [PubMed]
28. Szklarczyk D, Morris JH, Cook H, Kuhn M, Wyder S, Simonovic M, Santos A, Doncheva NT, Roth A, Bork P, Jensen LJ, von Mering C. The STRING database in 2017: quality-controlled protein-protein association networks, made broadly accessible. Nucleic Acids Res. 2017; 45:D362–68. https://doi.org/10.1093/nar/gkw937 [PubMed]
29. Shannon P, Markiel A, Ozier O, Baliga NS, Wang JT, Ramage D, Amin N, Schwikowski B, Ideker T. Cytoscape: a software environment for integrated models of biomolecular interaction networks. Genome Res. 2003; 13:2498–504. https://doi.org/10.1101/gr.1239303 [PubMed]
30. Tan B, Shi X, Zhang J, Qin J, Zhang N, Ren H, Qian M, Siwko S, Carmon K, Liu Q, Han H, Du B, Liu M. Inhibition of rspo-Lgr4 facilitates checkpoint blockade therapy by switching macrophage polarization. Cancer Res. 2018; 78:4929–42. https://doi.org/10.1158/0008-5472.CAN-18-0152 [PubMed]
31. Nielsen MM, Witherden DA, Havran WL. Γδ T cells in homeostasis and host defence of epithelial barrier tissues. Nat Rev Immunol. 2017; 17:733–45. https://doi.org/10.1038/nri.2017.101 [PubMed]
32. Cai XY, Wang JX, Yi Y, He HW, Ni XC, Zhou J, Cheng YF, Jin JJ, Fan J, Qiu SJ. Low counts of γδ T cells in peritumoral liver tissue are related to more frequent recurrence in patients with hepatocellular carcinoma after curative resection. Asian Pac J Cancer Prev. 2014; 15:775–80. https://doi.org/10.7314/apjcp.2014.15.2.775 [PubMed]
33. Yi Y, He HW, Wang JX, Cai XY, Li YW, Zhou J, Cheng YF, Jin JJ, Fan J, Qiu SJ. The functional impairment of HCC-infiltrating γδ T cells, partially mediated by regulatory T cells in a TGFβ- and IL-10-dependent manner. J Hepatol. 2013; 58:977–83. https://doi.org/10.1016/j.jhep.2012.12.015 [PubMed]
34. Karn T, Jiang T, Hatzis C, Sänger N, El-Balat A, Rody A, Holtrich U, Becker S, Bianchini G, Pusztai L. Association between genomic metrics and immune infiltration in triple-negative breast cancer. JAMA Oncol. 2017; 3:1707–11. https://doi.org/10.1001/jamaoncol.2017.2140 [PubMed]
35. Gabrielson A, Wu Y, Wang H, Jiang J, Kallakury B, Gatalica Z, Reddy S, Kleiner D, Fishbein T, Johnson L, Island E, Satoskar R, Banovac F, et al. Intratumoral CD3 and CD8 t-cell densities associated with relapse-free survival in HCC. Cancer Immunol Res. 2016; 4:419–30. https://doi.org/10.1158/2326-6066.CIR-15-0110 [PubMed]
36. Sideras K, Biermann K, Verheij J, Takkenberg BR, Mancham S, Hansen BE, Schutz HM, de Man RA, Sprengers D, Buschow SI, Verseput MC, Boor PP, Pan Q, et al. PD-L1, galectin-9 and CD8⁺ tumor-infiltrating lymphocytes are associated with survival in hepatocellular carcinoma. Oncoimmunology. 2017; 6:e1273309. https://doi.org/10.1080/2162402X.2016.1273309 [PubMed]
37. Xu X, Tan Y, Qian Y, Xue W, Wang Y, Du J, Jin L, Ding W. Clinicopathologic and prognostic significance of tumor-infiltrating CD8+ T cells in patients with hepatocellular carcinoma: a meta-analysis. Medicine (Baltimore). 2019; 98:e13923. https://doi.org/10.1097/MD.0000000000013923 [PubMed]
38. Ye L, Li Y, Tang H, Liu W, Chen Y, Dai T, Liang R, Shi M, Yi S, Chen G, Yang Y. CD8+CXCR5+T cells infiltrating hepatocellular carcinomas are activated and predictive of a better prognosis. Aging (Albany NY). 2019; 11:8879–91. https://doi.org/10.18632/aging.102308 [PubMed]
39. Yao Q, Bao X, Xue R, Liu H, Liu H, Li J, Dong J, Duan Z, Ren M, Zhao J, Song Q, Yu H, Zhu Y, et al. Prognostic value of immunoscore to identify mortality outcomes in adults with HBV-related primary hepatocellular carcinoma. Medicine (Baltimore). 2017; 96:e6735. https://doi.org/10.1097/MD.0000000000006735 [PubMed]
40. Sun C, Xu J, Song J, Liu C, Wang J, Weng C, Sun H, Wei H, Xiao W, Sun R, Tian Z. The predictive value of centre tumour CD8⁺ T cells in patients with hepatocellular carcinoma: comparison with immunoscore. Oncotarget. 2015; 6:35602–15. https://doi.org/10.18632/oncotarget.5801 [PubMed]

Research Paper Volume 12, Issue 24 pp 25172—25188