Angiogenesis-Related Molecular Subtypes and a Novel Prognostic Signature in Clear Cell Renal Cell Carcinoma Patients

Introduction

Renal cell carcinoma (RCC) accounts for 90% of adult renal malignancies, and is regarded as the most lethal urological cancer.^1,2 Approximately 70–80% of cases with RCC was clear cell RCC (ccRCC).^2,3 Tumor heterogeneity leads to different survival outcomes in ccRCC patients.⁴ At present, there are no consistent molecular subtypes and prognostic molecular characteristics of ccRCC.⁵ It is urgent to explore effective biomarkers to establish risk stratification method and molecular subtypes to further investigate the potential molecular mechanisms of ccRCC.

It has been reported that tumor angiogenesis plays a vital role in the proliferation and metastasis of ccRCC, as well as antiangiogenic agents.⁶ The expression of vascular endothelial growth factor receptors (VEGFRs) and hypoxia inducible factor (HIF) was significantly increased in ccRCC cells, which promotes tumor angiogenesis.⁷ In recent years, the American Food and Drug Administration (FDA) approved tyrosine kinase inhibitors (TKI), including sorafenib, sunitinib, pazopanib etc., as vital anti-cancer targeted drugs by inhibiting VEGFRs.⁸ TKI targeted therapy might promote the expression of major histocompatibility complex class I (MHC-I) and the infiltration of diverse immune cells, which might contribute to the therapeutic efficacy of immune checkpoint inhibitors.⁹ Motzer et al¹⁰ conducted an integrated evaluation of 823 advanced RCC patients and identified seven molecular subsets with distinct angiogenesis, immune, cell-cycle, metabolism, and stromal programs. These results would contribute to molecularly stratifying RCC patients, explaining improved prognosis of sarcomatoid tumors to immunotherapy versus antiangiogenics alone, and developing personalized therapies. Hence, there might be a very promising value to develop effective risk stratification tools and molecular subtypes for ccRCC patients by using angiogenesis-related genes (ARGs). However, to our knowledge, there are no published studies on developing angiogenesis related molecular subtypes and prognostic signatures of ccRCC from the perspective of ARGs.

In this study, we utilized consensus clustering analysis to identify two angiogenesis molecular subtypes and explored whether the molecular subtypes were associated with antitumor immunity. Besides, we developed a signature based on ARGs for predicting prognosis of subtype 1 ccRCC by using univariate and multivariate Cox regression analyses. Internal and external validation were then performed, and the correlations of the ARG-based signature with tumor microenvironment, immune cell infiltration and immunotherapy were investigated in particular.

Materials and Methods

Data Acquisition and Processing

We downloaded the transcriptome, mutation, copy number variation and clinical data of ccRCC from TCGA database (https://portal.gdc.cancer.gov). There were 539 ccRCC patients with unabridged transcriptome data, 537 ccRCC patients with complete clinical data in TCGA database. The clinical data included age, gender, grade, T stage, N stage, M stage, and AJCC stage. There was a total of 475 ccRCC patients with complete clinical data and 91 ccRCC cases with unabridged transcriptome data in International Cancer Genome Consortium (ICGC) database (https://dcc.icgc.org/). After we merged transcriptome and survival information, there was a total of 90 ccRCC cases with unabridged transcriptome and survival data simultaneously in ICGC database, and 530 ccRCC cases with unabridged transcriptome and survival data simultaneously in TCGA database. All these data were combined into a matrix file using Perl language (http://www.perl.org/). Moreover, gene names of Ensemble IDs were converted into gene symbols utilizing the Ensemble database (http://asia.ensembl.org/signature.html).

The Molecular Signatures Database (MSigDB) (https://www.gsea-msigdb.org/gsea/msigdb) is regarded as a collection of annotated gene sets utilizing GSEA software. We obtained a total of 36 ARGs from HALLMARK ANGIOGENESIS gene set (H: hallmark gene sets; systematic name: M5944) in MSigDB database. All these ARGs were up-regulated during formation of blood vessels (angiogenesis). The detailed list of these ARGs was presented in Supplemental Document.

Identification of Angiogenesis-Related DEGs and Functional Enrichment

The mRNA expression matrix of 36 ARGs in TCGA cohort was extracted. We utilized the “limma” R package and Wilcoxon test to filter angiogenesis-related differentially expressed genes (DEARGs) between the normal samples and ccRCC samples in TCGA cohort. The cut-off value was false discovery rate (FDR) < 0.05 and log2 |fold change (FC)| > 1. The Metascape database (http://metascape.org/) was utilized for functional enrichment analysis of DEARGs. The mutation status and copy number variation of these 36 ARGs were also explored using TCGA data.

Identifying Angiogenesis Molecular Subtypes Using Consensus Clustering Analysis

Firstly, univariable Cox regression analysis was used to screen prognostic DEARGs that were associated with OS of ccRCC. We then performed consensus clustering analysis for identifying angiogenesis-related molecular subtypes of ccRCC based on DEARGs by using R package “ConsensusClusterPlus”. Next, we used R package “survival” and “survminer” to investigate the association between angiogenesis molecular subtypes and OS. R package “pheatmap” was used to show the association between angiogenesis molecular subtypes and clinicopathologic features (including T stage, N stage, M stage, AJCC stage, grade, gender, and age).

Development of a Novel Angiogenesis-Based Prognostic Signature for Subtype 1 ccRCC

All ccRCC patients with subtype 1 were randomly divided into training cohort and test cohort. We then performed univariate and multivariate Cox regression analyses in training cohort to establish an angiogenesis-based signature for predicting OS in subtype 1 ccRCC. Patients were divided into two groups, high-risk group and low-risk group, according to median value of the risk score. Survival analysis and time-dependent receiver operating characteristic (ROC) curve were performed to validate the performance of this angiogenesis-based signature. Univariate and multivariate independent prognostic analyses were used to investigate whether the angiogenesis-based prognostic signature was an independent indicator of OS of ccRCC.

Internal and External Validation of This Angiogenesis-Based Prognostic Signature

To perform internal and external validation for this novel angiogenesis-based prognostic signature, we calculated risk score for each case by risk formula in test cohort, subtype 2 cohort, and ICGC cohort. Patients in test cohort, subtype 2 cohort, and ICGC cohort were also divided into two groups, high-risk group and low-risk group, according to median value of the risk score.

Survival analysis and time-dependent ROC curve were also performed to validate the performance of angiogenesis-based signature. Moreover, the performance of this signature was also compared with previous signatures using ROC curve. The expression heatmap, the distribution of risk score and survival time of training cohort, test cohort, whole subtype 1 cohort, and ICGC cohort were presented using “pheatmap” R package.

Exploration of the Associations of this Angiogenesis-Related Signature with Tumor Microenvironment, Immune Cell Infiltration and Immune Functions

The ESTIMATE algorithm was applied to access tumor microenvironment (TME) scores while the CIBERSORT method was utilized to evaluate the score of 22 types of immune infiltrating cells.^11,12 Then, we calculated the infiltrating score of 22 immune cells using CIBERSORT method and tumor microenvironment (TME) scores using ESTIMATE algorithm. The associations of angiogenesis molecular subtypes with immune cell infiltration and tumor microenvironment were explored. Generally, the expression levels of immune checkpoint gene are closely related to the sensitivity of immunotherapy. We investigated the association of the angiogenesis molecular subtypes with expression level of PD-L1 gene to predict immunotherapy sensitivity. P value < 0.05 was considered statistically significant.

Besides, we calculated the infiltrating score of 16 immune cells and the activity of 13 immune-related pathways of TCGA cohort using single-sample gene set enrichment analysis (ssGSEA) in the “gsva” R package. Then, the associations of the angiogenesis-based signature with immune cells infiltration and immune-related pathways activity were explored.

Validation of Hub Biomarkers in Multi-Database

Survival analysis based on GEPIA database (http://gepia.cancer-pku.cn/index.html) was performed to evaluate the prognostic value of these three hub biomarkers. UALCAN database (http://ualcan.path.uab.edu/) is a portal for facilitating tumor subgroup gene expression and survival analyses. Expression levels and promoter methylation levels of hub genes were revealed using UALCAN database.

Results

Angiogenesis-Related DEGs and Functional Enrichment

The clinicopathologic data of TCGA cohort and ICGC cohort were demonstrated in Table 1. There was a total of 36 ARGs (Supplemental Document). The mutation status and copy number variation of these 36 ARGs were presented in Figure 1A and B. A total of 17 DEARGs were finally identified, including 3 downregulated genes and 14 upregulated genes (Table 2). The results of functional enrichment analysis of these DEARGs were presented in Figure 1C and D. These 17 DEARGs were mainly enriched in extracellular structure organization, response to wounding, response to growth factor, NABA PROTEOGLYCANS, integrin-mediated signaling pathway, post-translational protein phosphorylation, platelet degranulation, regulated exocytosis, positive regulation of cytokine production, cartilage development, HTLV-I infection, transmembrane receptor protein serine/threonine kinase signaling pathway, and embryonic morphogenesis.

Table 1 Clinicopathologic Data of TCGA Cohort and ICGC Cohort

Table 2 Angiogenesis-Related Differentially Expressed Genes

Figure 1 The mutation status (A) and copy number variation (B) of these 36 angiogenesis-related genes (ARGs). Functional enrichment analysis for angiogenesis-related differentially expressed genes (DEGs) (C). P-value of each gene in the network (D).

Identifying Two Angiogenesis Molecular Subtypes of ccRCC

Univariable Cox regression analysis revealed that there was a total of nine angiogenesis-related DEGs associated with OS, including COL5A2, JAG2, LPL, LUM, MSX1, NRP1, PF4, PRG2, TIMP1. We identified a total of two angiogenesis molecular subtypes of ccRCC, including 213 cases of subtype 1 and 317 cases of subtype 2 (Figure 2A and B). The OS of subtype 1 ccRCC was significantly decreased compared with that of subtype 2 ccRCC (P=0.001, Figure 2C). The correlation heatmap between the angiogenesis molecular subtypes and clinicopathologic features was presented in Figure 2D. As indicated by heatmap, these two molecular subtypes have significantly different AJCC stage and T stage.

Figure 2 Consensus clustering analysis for identifying two angiogenesis molecular subtypes of ccRCC (A and B). Comparison of OS between these two subtypes (C). The correlation heatmap between this angiogenesis molecular subtypes and clinicopathologic features (D). The expression level of PD-L1 between these two subtypes (E). The expression correlation of PD-L1 with nine angiogenesis-related prognostic DEGs (F). **P<0.01; *P<0.05.

The expression level of PD-L1 was significantly decreased in subtype 1 group in comparison with that in subtype 2 group (P <0.01, Figure 2E). The expression correlation of PD-L1 with these nine angiogenesis-related prognostic DEGs (including COL5A2, JAG2, LPL, LUM, MSX1, NRP1, PF4, PRG2, TIMP1) was demonstrated in Figure 2F. The results showed that the expression level of PD-L1 was significantly associated with the expression level of JAG2 and NRP. The immune score, stromal score, and ESTIMATE score in subtype 1 group was significantly increased compared with that in subtype 2 group while the tumor purity score in subtype 1 group was significantly decreased compared with that in subtype 2 group (P<0.05, Figure 3A–D). The infiltrating proportions of naive B cells, regulatory T cells, M0 macrophages, M2 macrophages in subtype 1 group were significantly increased compared with that in subtype 2 group while the infiltrating proportions of CD8 T cells, T follicular helper cells, monocytes, M1 macrophages, resting dendritic cells in subtype 1 group were significantly decreased compared with that in subtype 2 group (P<0.05, Figure 3E).

Figure 3 Relationship between angiogenesis molecular subtypes and tumor immune microenvironment. Correlation of angiogenesis molecular subtypes with ESTIMATE score (A), tumor purity (B), immune score (C), stromal score (D). Correlation of angiogenesis molecular subtypes with tumor-infiltrating immune cells (E).

Development and Validation of a Novel Angiogenesis-Based Signature for Prognosis Prediction of Subtype 1 ccRCC

There was a total of 213 cases of type 1 ccRCC, including 108 in the training cohort and 105 in test cohort. A total of 90 ccRCC cases in ICGC database were used as an external validation cohort. We performed multivariate Cox regression analysis to develop a novel angiogenesis-based signature using three ARGs in training cohort for prognosis prediction of subtype 1 ccRCC. The calculation formula of risk score is shown as follows: risk score = 0.202457231804962 * MSX1 + 0.000522859217901573 * TIMP1 - 0.109529744114114 * JAG2 (Table 3).

Table 3 Univariate and Multivariate Cox Regression Analysis to Develop Angiogenesis-Based Prognostic Signature of Subtype 1 ccRCC

The difference in OS between high-risk and low-risk group was statistically significant in training cohort (P = 0.009), test cohort (P = 0.024), whole type 1 cohort (P < 0.001), and ICGC cohort (P = 0.041), respectively. High-risk score was associated with significantly poor OS in comparison with low-risk score in training cohort, test cohort, whole type 1 cohort and ICGC cohort, respectively. The area under ROC curve (AUC) for one year OS prediction was 0.732, 0.710, 0.725 and 0.645 in training cohort, test cohort, whole type 1 cohort and ICGC cohort, respectively, suggesting the potential power of this prognostic signature for prognosis prediction of subtype 1 ccRCC (Figure 4). The expression heatmap, the distribution of risk score and survival time of training cohort, test cohort, whole subtype 1 cohort and ICGC cohort were demonstrated in Figure 5.

Figure 4 Internal and external validation of a novel angiogenesis-based prognostic signature in subtype 1 ccRCC. The survival analysis between high and low risk group and corresponding area under ROC curve in training cohort (A and B), test cohort (C and D), the whole type 1 cohort (E and F) and validation cohort (G and H).

Figure 5 The expression heatmap, the distribution of risk score and survival time of training cohort (A), test cohort (B), the whole type 1 cohort (C) and ICGC cohort (D).

Independent Prognostic Analysis and Comparison with Previous Models

Univariate and multivariate independent prognostic analyses showed that this signature was an independent predictor for subtype 1 ccRCC, indicating the great performance of this signature (P < 0.05, Table 4). Besides, the difference in OS between high-risk and low-risk group of subtype 2 ccRCC in TCGA cohort was also statistically significant (Figure 6A). The area under ROC curve (AUC) of our signature for OS prediction was higher in comparison with YangSu’s signature, YueheGuo’s signature, JingminZhou’s signature, YueWu’s signature, suggesting that the power of this prognostic signature was superior to other signatures reported in previous studies (Figure 6B). Stage III–IV, grade 3–4 and T3-4 were significantly associated with higher risk score compared with stage I–II, grade 1–2 and T1-2, respectively (P < 0.05, Figure 7A–D). The correlation heatmap among this signature, clinicopathologic features, angiogenesis molecular subtypes and immune score was presented in Figure 7E.

Table 4 Univariate and Multivariate Independent Prognostic Analysis

Figure 6 External validation in subtype 2 ccRCC (A). Comparison of this signature with previous models (B).

Figure 7 Comparison of this signature with T stage (A), AJCC stage (B), grade (C) and PD-L1 expression (D). The correlation heatmap among this signature, clinicopathologic features, angiogenesis molecular subtypes and immune score (E).

Association of the Angiogenesis-Related Signature with Immune Cell Infiltration and Immune Function

We quantified the immune cell infiltration and immune functions using ssGSEA to further explore the association of this angiogenesis-based signature with immune status. Interestingly, the infiltrating proportions of aDCs, macrophages, Tfh, Th1 cells and Th2 cells in high-risk group were significantly increased compared with that in low-risk group, while the infiltrating proportions of mast cells in high-risk group were significantly decreased compared with that in low-risk group (P<0.05, Figure 8A). Moreover, the results also showed that patients with high-risk score have higher scores of immune function or pathways than those with low risk, including APC co-inhibition, CCR, para-inflammation, T cell co-inhibition, T cell co-stimulation. However, patients with high risk had a lower score of Type II IFN response than those with low risk (P<0.05, Figure 8B).

Figure 8 Association of this signature with immune cell infiltration (A) and immune functions (B). ***P<0.001; **P<0.01; *P<0.05.

Validation of Hub Biomarkers in Multi-Database

It was also found that expression levels of MSX1, TIMP1 and JAG2 were significantly increased in ccRCC tissues compared with that in normal tissues. Using ULCAN database, we indicated that MSX1 has significantly higher levels of promoter methylation in ccRCC tissues in comparison with that in normal tissue while the promoter methylation level of JAG2 was significantly lower in ccRCC tissues compared with that in normal tissue. However, the promoter methylation level of TIMP1 between ccRCC and normal tissue was not statistically significant (Figure 9). The survival analysis revealed that high expression of TIMP1 was significantly associated with poor OS and disease-free survival (DFS) of ccRCC patients, while high expression of JAG2 was significantly associated with higher OS of ccRCC patients. The expression level of MSX1 was not associated with OS and DFS of ccRCC patients. The expression level of JAG2 was not associated with DFS of RCC patients (Figure 10).

Figure 9 Validation of mRNA expression levels and promoter methylation of MSX1 (A and B), TIMP1 (C and D) and JAG2 (E and F) in ULCAN database.

Figure 10 Correlation of the expression levels of MSX1 (A and B), TIMP1 (C and D) and JAG2 (E and F) with OS and disease-free survival in GEPIA data.

Discussion

It has been confirmed that paracrine and autocrine factors promoting angiogenesis is vital in balancing cell proliferation, orientation, migration, and elongation, which might result in micro-vessels’ formation and establishment.¹³ Increasing evidence indicates that aberrant angiogenesis in tumor microenvironment is of great significance in processes of carcinogenesis, progression, and metastasis of ccRCC.¹⁴ However, there is no published study focusing on the correlation of angiogenesis with prognosis evaluation and molecular subtypes of ccRCC. In this study, we identified a total of two angiogenesis molecular subtypes of ccRCC using consensus clustering analysis. The OS of subtype 1 ccRCC was significantly decreased compared with that of subtype 2 ccRCC. Further investigation showed that these two molecular subtypes had significantly varying tumor microenvironments. The expression level of PD-L1 was significantly decreased in subtype 1 group in comparison with that in subtype 2 group.

Most importantly, we developed a novel prognostic signature related to angiogenesis for subtype 1 ccRCC. Univariate and multivariate Cox regression analyses revealed that this novel prognostic signature was a vital independent predictor for OS of ccRCC. Internal and external validation indicated the potential role of this angiogenesis-related signature for subtype 1 ccRCC. There have been several prognostic models reported previously. Chu et al¹⁵ constructed a novel seven protein signature for predicting prognosis. Xing et al utilized glycolysis-related genes to successfully establish and externally verify a novel signature for predicting OS of ccRCC.¹⁶ Wang et al¹⁷ successfully constructed a prognostic model for ccRCC based on apoptosis-related genes, which would help clinicians make the prognostic assessment. However, there has been no study developing prognostic signature of ccRCC from the point of angiogenesis and molecular subtypes. In this study, we developed a novel angiogenesis-related prognostic signature for subtype 1 ccRCC for the first time. Most importantly, both internal validation and external validation suggested the satisfactory predictive power of this signature. Furthermore, we interestingly found that the area under ROC curve (AUC) of our signature for OS prediction was higher in comparison with YangSu’s signature,¹⁸ YuheGuo’s signature,¹⁹ JingminZhou’s signature,²⁰ YueWu’s signature,²¹ suggesting that the power of this prognostic signature might be superior to other signatures reported in previous studies. Furthermore, we demonstrated that stage III–IV, grade 3–4 and T3-4 were significantly associated with higher risk score compared with stage I–II, grade 1–2 and T1-2, respectively, and that high and low-risk score ccRCC patients have significantly varying immune cell infiltration and immune function.

These three vital ARGs (including MSX1, TIMP1 and JAG2), which were used to construct angiogenesis-related prognostic signature, could be considered as candidate markers of ccRCC. Dias et al²² found that extracellular vesicle-derived TIMP1 in plasma could be used as an adverse potential prognostic biomarker of ccRCC. Armstrong et al²³ detected 23 types of angiokines in ccRCC patients in a randomized open-label Phase II trial, and found that the expression level of TIMP1 increased with progression on everolimus or sunitinib treatment. He et al²⁴ reported that RUNX3 played a vital role in inhibiting ccRCC proliferation and metastasis of cRCC cells by regulating TIMP1. Yang et al²⁵ considered that IL-8 could recruit tumor associated neutrophils and induce the expression of JAG2 in tumor associated neutrophils, which play an important role in regulating immune evasion microenvironment in epithelial ovarian cancer. Hatano et al revealed that JAG2 overexpression has potential prognostic value in oral squamous cell carcinoma.²⁶ He et al²⁷ demonstrated that interaction of JAG2 and PRAF2 facilitates colorectal carcinoma cells’ migration and invasion through non-canonical Notch and non-EMT-dependent pathways. MSX1 could inhibit angiogenesis to suppress tumorigenesis.²⁸ Son et al²⁸ showed that the interaction of PIASy and MSX1 mediated the inhibition of tumor angiogenesis. Yue et al²⁹ reported that msh homeobox 1 (MSX1) could suppress breast tumorigenesis through arresting G1/S cell cycle and apoptosis, and the methylation of MSX1 might be regarded as an epigenetic biomarker for early diagnosis of breast cancer. However, there has been no study focusing on the role of MSX1 and JAG2 in ccRCC. In this study, we found that expression levels of MSX1, TIMP1 and JAG2 were significantly increased in RCC tissues compared with normal tissues. High expression of TIMP1 was significantly associated with poor OS and DFS of ccRCC patients while high expression of JAG2 was significantly associated with higher OS of ccRCC patients. Besides, we revealed that MSX1 has significantly higher levels of promoter methylation in ccRCC tissues in comparison with normal tissue, while the promoter methylation level of JAG2 was significantly lower in ccRCC tissues compared with normal tissue.

However, we have to acknowledge that several limitations in our study were inescapable. First of all, retrospective design with limited sample size would result in selection bias. Larger prospective studies should be performed for generalizability. Secondly, all of the work in this study was purely based on bioinformatics analysis and did not achieve any mechanistic explanation. Further in vitro and in vivo experiments into the underlying mechanism are warranted.

In this study, we identified two angiogenesis-related molecular subtypes for ccRCC. These two molecular subtypes have significantly different prognosis, tumor microenvironment, and immune checkpoint inhibitor sensitivities. Besides, we developed and successfully validated a novel angiogenesis-related gene signature for predicting prognosis of subtype 1 ccRCC, which was superior to several previous signatures. The angiogenesis-based molecular subtypes and prognostic signature would contribute to further understanding of molecular mechanism in ccRCC.