BMP2/BMPR1A is linked to tumour progression in dedifferentiated liposarcomas

Data acquisition

Protein expression data was downloaded from the Human Protein Atlas (Uhlén et al., 2010; Uhlén et al., 2015) (Cancer Atlas and the Cell Line Atlas; proteinatlas.org). Pan-cancer analysis was carried out using cBioPortal (cbioportal.org) (Cerami et al., 2012; Gao et al., 2013). With pathway data obtained from Cytoscape 3.2.1 (cytoscape.org) (Lopes et al., 2010). RNASeqV2 normalised gene results and clinical metadata were downloaded form the TCGA Data Portal (tcga-data.nci.nih.gov) (Sarcoma (provisional), Prostate (provisional), Breast (The Cancer Genome Atlas Network, 2012), Ovarian (The Cancer Genome Atlas Network, 2011), Lung (The Cancer Genome Atlas Network, 2014) and Head and Neck (The Cancer Genome Atlas Network, 2015)).

Data analysis

Normalised RNASeqV2 expression values and corresponding clinical metadata were manipulated in R Studio (Mac) 0.99.484 (R Studio: https://www.rstudio.com/). Combined data were analysed in Microsoft Excel (Mac 14.4.3) and R Studio with results plotted in GrpahPad Prism 6 (Mac). Comparisons between two groups were made using two-tailed t-tests or Mann–Whitney U-tests and more than two groups by one-way ANOVA or Kruskal–Wallis, multiple comparisons were corrected for by Bonferroni post hoc tests. Kaplan–Meier disease free survival graphs were constructed in GraphPad Prism 6 and presented as whole population versus experimental population. Comparisons between curves were made by Log-rank (Mantel-Cox) tests with the Mantel-Cox p value reported.

Differential expression analysis

For patient groups with observed differences in phenotype (pathway activation/disease free survival) differential expression analysis was performed on normalised RNASeqV2 gene counts. Data were combined in R studio and means across rows (genes) calculated. Multiple two tailed t-tests with Bonferroni correction were used to compare means. Using a false discovery rate of 0.05, differentially expressed genes were ranked based on log₁₀ fold-change. Ranked genes were analysed through the Gene Set Enrichment pre-ranked plugin (GSEAPreranked at genepattern.broadinstitute.org) with default parameters and 1,000 permutations (GenePattern version 3.9.4 (Hubble et al., 2009)). Gene sets c2.all.v5.0.symbols.gmt (Curated) and c6.all.v5.0symbols.gmt (Oncogenic signatures) were used for comparison. Gene sets in the top 10 from the dataset based on enrichments score were considered relevant.

Differential type 1 BMP receptor expression in normal tissues and sarcomas

Type I BMP receptors (BMPR1A and 1B) have well defined overlapping functions (Yoon et al., 2005) but may not be functionally redundant in all aspects of their biology. To begin to test this hypothesis, we first profiled expression of the type I and type II BMP receptors across normal tissues (Fig. S1). Whilst BMPR1A and BMPR1B showed similar expression across some tissues (for example heart, lung, skin, colon etc.) we observed differential expression in numerous tissues. For example, BMPR1B is highly expressed in all regions of the brain (cerebellum, hippocampus, cerebr BMP2/BMPR1A is linked to tumour progression in dedifferentiated liposarcomasal cortex and lateral ventricle) whereas BMPR1A has only moderate expression in the cerebral cortex, suggesting that BMPR1B signalling predominates in the central nervous system.

We next profiled expression patterns across cancer cell lines (Fig. S2 and Fig. 1A). In the three cell lines profiled by immunohistochemistry (IHC), BMPR1A exhibited both nuclear and cytoplasmic localisation, which matched BMPR2 expression patterns (Fig. S2). In contrast BMPR1B appeared predominantly localised to the nucleus, but not nucleolus. The osteosarcoma cell line U-2 OS showed highest expression of BMPR1A by IHC (Fig. 1A), which was echoed in RNAseq data from three sarcoma cell lines (Fig. 1B), where BMPR1A expression far exceeded that of BMPR1B. These data suggest that BMPR1A and BMPR1B may not be functionally redundant in all aspects of their biology.

As the largest differences in receptor expression patterns were observed in sarcoma cell lines and normal soft tissues, we elected to focus our analysis on BMP2 in the soft tissue sarcoma TCGA dataset. The observations from Figs. 1A and 1B (together with Figs. S1 and S2) that BMPR1A and BMPR1B have differential expression patterns led us to focus our analysis further on the effects of BMPR1A versus BMPR1B biased BMP2 signalling within this group of patients. As a proxy for this measure, we sought to identify patients with high BMP2 expression in the context of either high BMPR1A or BMPR1B expression.

By stratifying the 263 samples represented in the TCGA sarcoma dataset based on the ratio of BMPR1A to BMPR1B mRNA, we were able to define two clear groups (Fig. 1C). Although most patients expressed relatively similar levels of the two type I receptors there were clear groups who favoured either BMPR1A (red) or BMPR1B (blue) but who did not show differential expression of BMP2 (Fig. 1D). Based on the observed patterns of differential expression in normal tissues (Fig. S1) and cell lines (Fig. S2), and the conclusion that BMPR1A and BMPR1B may not be entirely functionally redundant, we reason that these two patient groups may show distinct outcomes to high BMP2 expression.

BMP2/BMPR1A expression is associated with extracellular matrix remodelling

To test the hypothesis that patients favouring BMPR1A would have differential biological responses to high levels of endogenous BMP2 to patients favouring BMPR1B we refined our patient groups from Fig. 1C. Referred to as BMPR1A-biased or BMPR1B-biased BMP2, we defined groups of patients in the upper quartile of the BMPR1A/BMPR1B ratio and in the upper quartile of BMP2 expression together with their reciprocals (high 1B/1A with high BMP2).

We first asked whether BMPR1A-biased BMP2 expression led to differential gene expression to BMPR1B-biased. In total we found 922 differentially expressed genes out of a possible 20,533 annotated transcripts, using a false discovery rate (FDR) of 0.05 (corrected p < 0.05). Of these, 252 were up-regulated and 669 down-regulated. Fig. 2A shows a volcano plot relating fold change (in patients with BMPR1A-biased versus BMPR1B-biased BMP2 expression) to p value (multiple two-tailed t-tests with Bonferroni post hoc). The insert panel in Fig. 2A highlights a selection of the most differentially expressed genes; the full list is available as Table S1. Amongst the genes negatively associated with BMPR1A-biased expression are WNT7B which has a role in osteoblast differentiation and mineralised nodule formation (Shao et al., 2016), whilst positively associated is SPON1 which inhibits excess bone formation (Palmer et al., 2014).

Although these data suggest that BMPR1A-biased expression is associated with reduced bone formation, which may be of importance to understanding rhBMP2 use in spinal surgeries, we sought to carry out a more unbiased analysis of our differential expression results in the form of a gene set enrichment analysis (GSEA). Figs. 2B and 2C show two of the most highly enriched gene sets, each associated with the matrisome, a computationally predicted collection of extracellular matrix (ECM) constituents and their regulatory networks (Naba et al., 2012).

Together, these data suggest that BMPR1A-biased BMP2 expression is associated with reduced bone formation and an increased turnover of ECM components. The interaction of tumour cells with ECM proteins can have profound effects on their biology, regulating gene expression, signal transduction and migration/invasion (Cassidy, 2014; Lu, Weaver & Werb, 2012). Moreover, BMP2 signalling has previously been linked to increased ECM production and degradation (Davidson et al., 2007). Thus, we next sought to determine whether BMPR1A-biased BMP2 expression had any clinical phenotype.

BMPR1A-biased BMP2 expression is prognostic of poor clinical outcome

Patients in the TCGA sarcoma dataset have a median disease free survival of 37 months (Fig. 3A). Those patients with BMPR1B biased expression (n = 44) have no difference in their prognosis when compared to the dataset as a whole. However, those with BMPR1A-biased BMP2 expression (n = 32) have a significantly worse outcome when disease free survival is used as the primary endpoint (p = 0.0087). To discount prognostic influences of BMPR1A, BMPR1B or BMP2 expression individually, we constructed Kaplan-Myer curves for the entire sarcoma dataset and compared patients with the highest and lowest 10% expression of BMPR1A, BMPR1B, BMPR2, BMP2 and BMP7. Despite several interesting trends, no significant associations were seen (Fig. S3). Thus, in addition to ECM remodelling and reduced bone formation, soft tissue sarcoma patients with BMPR1A-biased BMP2 expression are more likely to show tumour progression in the presence of high levels of endogenous BMP2.

Soft tissue sarcomas are classified based on their origin in the soft but not connective tissues of the body. They can be separated based on their cell of origin and are generally considered a heterogeneous mixture of diseases (Cormier & Pollock, 2004). Figure 3B shows the TCGA dataset broken down into subtypes of soft tissue sarcomas. The two best-represented subtypes are Leiomyosarcoma (LMS; smooth muscle sarcomas) and Dedifferentiated liposarcoma (DLS; adipocyte sarcomas), though the majority (60%) of the dataset is unclassified. Disease free survival is shown in Fig. 3C; LMS tends to have a better prognosis than DLS though data for the three other subtypes is limited. To determine whether BMPR1A biased expression was localised to a specific subtype, we profiled expression of BMPR1A, BMPR1B, BMPR2 and BMP2 across the three best-represented subtypes (UPS, LMS and DLS). No significant differences were noted, though UPS did trend toward lower overall expression of both ligand and receptors (Figs. S4). The two largest groups, LMS and DLS had indistinguishable expression values of BMP signalling nodules. Interestingly, BMPR1A-biased BMP2 expression conferred a worse prognosis on patients than ether LMS or DLS. Thus, we next sought to determine whether the prognostic benefits of BMPR1A-biased signalling were split evenly across these two major subtypes.

BMPR1A prognostic signalling is localised to dedifferentiated liposarcomas

We next profiled disease free survival in patients with LMS or DLS and investigated the effects of either BMPR1A-biased BMP2 expression or high BMPR1A expression alone. BMPR1A and BMPR1A-biased signalling appeared to have little effect in terms of disease free survival in patients with LMS (Fig. 4A), though high expression and R1A-biased BMP2 was highly detrimental in patients with DLS (Fig. 4B).

To investigate the mechanism for this effect, we compared global gene expression in DLS patients with BMPR1A-biased signalling versus those with BMPR1B-biased signalling. Despite working with a much reduced dataset of 10 versus 10 patients, we found 1,375 differentially expressed genes with a FDR of 0.05, the vast majority (1,061) of which were down-regulated (Fig. 4C); the full list is available as Table S2. Such a high number of differentially expressed transcripts compared to our early analysis in the whole sarcoma dataset (Fig. 2A) suggests patient groups which are more biologically robust. Amongst these differentially expressed genes, we noted a positive association between BMPR1A-biased BMP2 expression and GPR176, an orphan G-protein coupled receptor thought to be involved in atherosclerosis (Karagiannis et al., 2013). We performed GSEA on these differentially expressed genes and found a positive enrichment for a number of extracellular matrix and matrisome gene sets (Fig. 4D). These results further substantiate the observation that BMPR1A-biased BMP2 expression is associated with reduced bone formation and an increased turnover of ECM components as well as tumour progression.

Finally, we sought to determine whether the detrimental effects of BMP2-BMPR1A expression was a pan-cancer phenomenon. Several large studies performed by The Cancer Genome Atlas (TCGA) show copy number changes and mutational events in BMP2 interacting partners as defining features of that malignancy (Fig. S5). For example, two separate studies in Ovarian cancer show copy number gains as frequent events, which we later determined to be due to amplifications in the BMP2 and BMP7 loci (Fig. S6). Interestingly, mRNA up-regulation of both ligands and receptors were frequent events in each of the cancers profiled. Deep deletion of BMPR1A was a frequent event in prostate cancers, this has previously been shown to be related to the loss of the PTEN locus, which is associated with invasive forms of prostate cancer (Yoshimoto et al., 2012). In order to test whether our identified BMPR1A-biased BMP2 expression biomarker was useful across tumour types, we investigated disease free survival in Ovarian, Lung, Pancreatic and Prostate cancers, but found no significant alterations (Figs. S5). Thus, it appears that the detrimental effects of BMPR1A biased BMP2 expression are confined to the context of soft tissue sarcomas and in particular dedifferentiated liposarcomas.