Primary prostate tumor VEGF/Sema alterations
To examine why prostate cancer is less susceptible to VEGF inhibition, we first compared gene expression between primary prostate tumors and a cancer type that typically does respond to VEGF inhibitors, renal cell carcinoma [25]. Using the TCGA RNA-Seq dataset, we found that while two key pro-angiogenic ligands, VEGFA and PGF, were up-regulated in renal cell carcinoma, they were down-regulated in primary prostate adenocarcinoma (Fig. 1a, b for prostate, Fig. 1e, f for kidney). This could indicate a lack of VEGF signaling for VEGF inhibitors to target, making attempts at targeting the VEGF pathway in prostate cancer futile. This pattern was observed for other VEGF ligands: these genes were down-regulated or unchanged in prostate cancer, whereas they were up-regulated in renal cell carcinoma (Additional file 1: Figures S2-S5). Similarly, the receptors VEGFR2 (KDR) and NRP1 were decreased or relatively unchanged in prostate primary tumors while they were up-regulated in renal cell carcinoma (Fig. 1c, d for prostate, Fig. 1g, h for kidney). This pattern was consistent across other VEGF receptors (Additional file 1: Figures S2-S5). The up-regulation of both ligands and receptors in renal cell carcinoma clearly suggests the importance of VEGF signaling in that case. However, the up- or down-regulation of VEGF ligands and receptors in cancer types may not fully explain the response to VEGF inhibitors in prostate cancer. Subsets of tumors may have pro-angiogenic gene expression despite the anti-angiogenic pattern in the tumors as an overall group. Additionally, other pathways that interact with the VEGF pathway may affect signaling.
We considered an additional set of genes that have been demonstrated to affect VEGF signaling: the semaphorins and their plexin receptors, bringing the total number of VEGF/Sema-related genes under consideration to 39. We expanded our comparison of normal prostate tissue and primary prostate tumors to include multiple microarray datasets in addition to the TCGA dataset (Additional file 1: Table S1, note that only GSE6919, GSE21034, TCGA, and GSE35988 contain both normal prostate tissue and primary prostate tumors). For primary tumor versus normal tissue comparisons, we measured differences in gene expression using a two-tailed t-test with multiple testing correction using the Benjamini-Hochberg procedure. We found that in addition to VEGF ligands, many semaphorins were also down-regulated in prostate cancer (Fig. 1i). In particular, class 3 semaphorins, which have potential anti-angiogenic effects due to their ability to compete with VEGF for neuropilin binding, were down-regulated. Not all gene expression changes were consistent across all datasets; ones that were significant (defined as q-value less than 0.05) in two or more datasets were marked with blue or red boxes to the right. The subset of genes consistently down-regulated in primary tumors included three out of five VEGF ligands (VEGFA, VEGFB, and VEGFC) and five out of seven class 3 semaphorins (SEMA3A, SEMA3B, SEMA3C, SEMA3D, and SEMA3E). These results made the overall impact on angiogenesis unclear as both pro-angiogenic VEGF signals and anti-angiogenic semaphorin signals were reduced in primary tumors.
A pro-lymphangiogenic gene expression signature is associated with aggressive primary tumors
To assess whether some subsets of primary prostate tumors may have differing potential benefits from VEGF signaling inhibitors, we used partial least squares discriminant analysis (PLS-DA), a multivariate algorithm that allowed us to simultaneously consider both the pattern of VEGF/Sema gene expression and effects on an output variable. As an output, we used a binary variable indicating whether biochemical recurrence (BCR) eventually occurred. Data for BCR/follow-up times were available in the TCGA and GSE21034 datasets, with 10 and 27 BCR events, respectively (Additional file 1: Table S3). Since BCR-negative samples with short follow-up times may eventually undergo recurrence, we used only the ten and 27 BCR-negative samples with the longest follow-up times in the TCGA and GSE21034 datasets, respectively. These groups of indolent samples had follow-up times greater than 3.3 years in the TCGA dataset and greater than 5.2 years in the GSE21034 dataset. These minimum indolent follow-up times were greater than all but one of the times at which BCR occurred in the TCGA dataset and all but three in the GSE21034 dataset. The TCGA PLS-DA model was effective in differentiating aggressive and indolent tumors, with a 92 ± 3 % training accuracy (Fig. 2a). The GSE21034 PLS-DA model, on the other hand, was less effective with only a 66 ± 1 % training accuracy (Fig. 2b). This discrepancy in performance between TCGA and GSE21034 datasets was seen for other biomarkers meant to distinguish aggressive from indolent tumors (Additional file 1: Figure S7). Notably, the VEGF/Sema PLS-DA model (and PLS-DA models based on the genes from other biomarkers) lost some of its predictive ability when correcting for Gleason score, but still was significantly prognostic (Additional file 1: Figure S8). ROC curves from leave-one-out cross-validation barely deviated from a 45-degree line (Additional file 1: Figure S6), indicating that these models would likely be ineffective in distinguishing aggressive and indolent tumors in other datasets. ROC curves for PLS-DA models trained with the 478 genes of the “angiome”, a set of angiogenesis-related genes [26], improved but were still fairly poor: the AUCs for TCGA and GSE21034 were 0.72 and 0.66, respectively (Additional file 1: Table S6). The 1,233 genes of the extended angiome [26] also led to poor PLS-DA models, with AUCs of 0.74 and 0.69 (Additional file 1: Table S7). Nonetheless, the PLS-DA models provided gene expression signatures with potential prognostic significance: when the PLS discriminant scores were used as predictive variables in Kaplan-Meier survival analysis, aggressive and indolent tumors had significantly different outcomes (Fig. 2c, d).
The PLS-DA models also yielded information regarding patterns of VEGF/Sema expression associated with aggressive prostate tumors. The gene-specific loadings vectors (arrows in Fig. 2a, b) are the projections of each gene onto the latent variable space, and show the extent to which different genes contribute to the variability in the response variable. The loadings revealed that the association of VEGF/Sema genes with aggressiveness was different between the two datasets. For the TCGA dataset, aggressive tumors were associated with high expression of FIGF (VEGFD), NRP2, PLXNA1, PLXNB1, PLXNB3, and SEMA5B and low expression of PLXNB2 and SEMA4D. The first two of these genes, FIGF (VEGFD) and NRP2, mediate pro-lymphangiogenesis signals. The process of lymphangiogenesis may provide tumors with a route by which to escape their tissue of origin into the bloodstream [27]. In the GSE21034 model, high expression of VEGFC, KDR, and NRP1 and low expression of VEGFA and SEMA4B were associated with aggressive tumors. Although this signature was different from the one found in the TCGA dataset, the high expression of VEGFC also suggested a potential for lymphangiogenic activity in aggressive primary tumors.
Metastatic prostate tumors are associated with a pro-angiogenic signature
While primary tumors are often treated with surgery and radiation, targeted therapeutics such as VEGF inhibitors are used more often in metastatic disease. Therefore, we next considered VEGF/Sema gene expression in datasets with metastatic tumors (Additional file 1: Table S1, note that only GSE6919, GSE21034, GSE32269, and GSE35988 contain both primary and metastatic tumors, while only GSE6919, GSE21034, GSE38241, and GSE35988 contain both normal prostate tissue and metastatic tumors). In contrast to the reduced expression of VEGF ligands in primary tumors, metastatic tumors tended to have higher expression of the major pro-angiogenic ligand, VEGFA (Fig. 3a, b). This ligand was up-regulated in metastases relative to primary tumors and normal tumors in the GSE6919, GSE35988, and GSE38241 datasets, but was actually down-regulated in GSE21034 and GSE32269. Notably, metastatic samples in the three datasets with up-regulated VEGFA were all obtained from warm autopsy programs where samples were processed rapidly upon the death of the patient. Metastatic samples in GSE32269 were from bone marrow biopsies of live patients, and no details were given regarding how metastatic samples were obtained in the GSE21034 dataset. The class 3 semaphorins were down-regulated in metastases relative to normal prostate tissue (Fig. 3b), suggesting that the loss of semaphorin expression in primary prostate tumors was maintained upon metastasis. SEMA3C was further down-regulated relative to primary tumors as well. Other expression alterations recurrent across datasets included up-regulation of NRP1, PLXNA1, and PLXNA3 relative to both normal tissue and primary tumors. These three genes participate in class 3 semaphorin signaling, while only NRP1 participates in VEGF signaling. KDR and NRP2 were recurrently down-regulated in metastases relative to normal tissue but not relative to primary tumors.
We used PLS-DA to show that the VEGF/Sema gene expression alterations present in metastatic prostate tumors differed significantly from both normal tissue and primary tumors. PLS-DA models comparing metastases with primary tumors in the GSE35988 dataset and with normal tissue in the GSE38241 dataset led to large separation between the PLS scores of the two classes (Fig. 3c, d and Additional file 1: Figure S9). The leave-one-out cross-validation accuracy was 100 % for both of these comparisons.
VEGF/Sema alterations are consistent across multiple metastases within individual patients
The metastasis samples in the GSE38241 dataset included three to four metastases per patient from five different patients, providing the opportunity to analyze VEGF/Sema gene expression in both inter- and intra-patient contexts. We observed a high degree of consistency between metastases from the same patient, with VEGFA consistently up-regulated and several class 3 semaphorins, KDR, and NRP2 consistently down-regulated (Fig. 4a; the numbers 16, 17, 21, 22, and 30 on the x-axis correspond to different patients). Some genes that were not significantly altered when comparing all metastases to normal samples did show patient-specific alterations: VEGFC and SEMA6A were up-regulated in some metastases and down-regulated in others. SEMA6A was consistent within patients, while VEGFC was not.
To assess the overall consistency of metastases from individual patients, we performed consensus K-means clustering of the metastatic samples and compared the clusters to the patients of origin. Repeating the clustering on random subsets of the data yielded a consensus matrix with clear separations of the groups when the number of clusters was five (Fig. 4b, c and Additional file 1: Figure S10). This did not perfectly separate the patients, but there was a significant association between clusters and patients (χ2 test p = 0.001). We also found that if we expanded our set of genes beyond VEGF/Sema to include other ligands and receptors important in angiogenesis (specifically, the EGF/ErbB family, HGF/Met, the Ang/Tie family, the PDGF/PDGFR family, and the IGF/IGFR family), we obtained consensus K-means clusters that resulted in perfect separation of the patients into clusters (Fig. 4d, e). This supports the hypothesis that, although signaling pathways that contribute to angiogenesis may vary from patient to patient, the multiple metastases that can arise from a single prostate tumor would likely all respond (or not respond) to the same or similar therapies. That this broader set of receptor tyrosine kinases and their associated ligands improves patient stratification is consistent with the need to understand a broader angiome beyond VEGF.
Computational modeling of the VEGF/Sema pathway stratifies patients
To further analyze the contributions of gene expression changes to angiogenesis-related signaling, we developed a computational model consisting of ordinary differential equations (ODEs) describing the ligand-receptor binding kinetics of the five VEGF ligands and seven class 3 Semaphorins. A detailed model is essential, given the high inter-individual variability in gene expression, and given that predicting the outcome on signaling of simultaneous changes in expression of multiple genes encoding competing ligands and receptors becomes difficult. For a more detailed description of the model development, see
Methods
and Additional file 2. A key component of the model is that it incorporates both tumor cells and tumor endothelial cells, and these cells can express both the ligands and receptors of the VEGF and Sema pathways. Instead of using an ‘average’ model, the tumor cell ligand and receptor production rates were varied based on the gene expression data in the GSE35988 dataset, which included normal prostate, primary prostate tumors, and metastatic tumors. This enabled us to run many simulations – one for each individual – and to predict the amounts of VEGF and Semaphorin signaling complexes in a patient-specific manner.
The results of simulating many individuals with tumors were divided into groups (benign, localized, metastatic) and the distribution of predicted signaling outputs across the population are calculated. Binding of the two major isoforms of VEGFA (VEGF165 and VEGF121) to VEGFR-2 and VEGFR-1 on endothelial cells was lowest in primary (localized) tumors, whereas the collective binding of the class 3 Semaphorins to endothelial plexins decreased in primary tumors relative to normal tissue and decreased further in metastatic tumors (Fig. 5a). The receptor binding trends closely followed ligand expression levels, with VEGFA expression being the predominant factor driving VEGFR-1/2 binding, while the total Sema3 expression accounted for most of the variation in Sema3-NRP-PlxnA binding (Fig. 5b). The receptor binding profiles on tumor cells (Additional file 1: Figure S12) were similar to those on endothelial cells. These results indicated that primary tumors are associated with both pro-angiogenic (decreased Sema3-NRP-PlxnA binding) and anti-angiogenic (decreased VEGFR-2 binding) alterations, making it difficult to predict whether primary tumors would benefit from therapies that inhibit VEGF signaling. On the other hand, the alterations in metastatic tumors were all pro-angiogenic: VEGFR-2 binding was higher and Sema3-NRP-PlxnA binding was lower. The range of VEGFR-2 binding was highest in metastatic samples, going beyond both the low and high ends of the range of normal samples. This suggested an important role for patient selection in the use of anti-angiogenic therapies, as only the patients with high baseline VEGFR-2 signaling would be expected to respond.
Aside from anti-angiogenic signaling initiated by Sema3 binding to neuropilins and plexins, Sema3s also may inhibit angiogenesis by competitively displacing VEGF from neuropilins. In our model, we found only a weak competitive effect: the least squares fits for the non binding ligand-receptor pairs (VEGFA effects on Sema3-NRP-PlxnA and Sema3 effects on VEGFA-VEGFR-2 in Fig. 5b) had slight negative slopes, but the least squares models explained a very small proportion of the overall variance. This was due to the fact that the amount of NRP1 and NRP2 present in ligand-containing signaling complexes was low relative to the total amount of neuropilins present on endothelial cells. If the amount of endothelial NRP1 and/or NRP2 were lower or the ligand secretion rates were higher, competition would be expected to have more impact. Thus, although direct VEGF-Sema competition for neuropilin is included in the model, the quantitative impact of the competition is predicted to be small.
Our simulation results suggest that determining patients with sensitivity to anti-angiogenic therapies would require biomarkers consisting of multiple predictor variables. Each individual tumor has a different predicted level of VEGFA-VEGFR-2 and low Sema3-NRP-PlxnA activity (Fig. 5c). The red and blue gradients, and the arrows on the axes, indicate direction of increasing pro-angiogenic signaling for VEGFR2 (accelerator ‘ON’) and PlxnA1 (brake ‘OFF’), respectively. Metastatic samples predominantly have high simulated VEGFA-VEGFR-2 and low Sema3-NRP-PlxnA (lower right quadrant of Fig. 5c: accelerator ‘ON’, brake ‘OFF’), and would be expected to benefit the most from VEGF-targeting therapies. However, this is not true of all metastatic tumors; a minority of patients with high VEGFA-VEGFR-2 also have high Sema3-NRP-PlxnA, possibly negating any clinical benefit of an anti-VEGFA agent. That combination is also predicted to be the most common for the benign tissues (Fig. 5c, green symbols). There are also a subset of metastatic tumors predicted to have low VEGFA-VEGFR2 signaling. This heterogeneity amongst the metastatic tumors reinforces the need to bring an individualized understanding to therapeutic selection. For those with localized disease (Fig. 5c, orange symbols) the characteristic signaling state is low VEGFA-VEGFR2 and high Sema-Plexin (upper left quadrant of Fig. 5c: accelerator ‘OFF’, brake ‘ON’), correlating with low angiogenesis potential and no metastases. Further analysis of this model and availability of both gene expression and clinical anti-VEGF outcome data will allow us to develop simulation-based biomarkers.