High expression of bone morphogenetic protein (BMP) 6 and BMP7 are associated with higher immune cell infiltration and better survival in estrogen receptor‑positive breast cancer.

There are conflicted experimental results on the role of bone morphogenetic proteins (BMPs) in cancer. Some results suggest that BMPs act as tumor suppressors while other findings indicate that BMPs are oncogenic factors. In the present study, we aimed to investigate the association of BMP expression and the survival of breast cancer patients utilizing The Cancer Genome Atlas (TCGA). High expression levels of BMP1 (P<0.001), BMP3 (P=0.002), BMP5 (P=0.002), BMP7 (P<0.001) and BMPR1A (P<0.001) were associated with better overall survival (OS). On the other hand, high expression levels of BMP6 (P<0.001), BMP8A (P=0.031), BMP8B (P<0.001) and BMPR1B (P=0.005) were associated with worse OS. Most of the BMPs demonstrated differential expression levels between normal and tumor tissues. High expression of BMP7 was found to be significantly associated with better prognosis in both estrogen receptor (ER)‑positive (ER+) (P<0.001) and ER‑negative (ER‑) tumors (P<0.001), whereas high expression of BMP6 was associated with better prognosis only in ER+ tumors (P=0.004); it was associated with worse prognosis in ER‑ tumors (P=0.006). In the ER+ tumors, high expression of BMP7 as well as BMP6 were associated with higher infiltration of anticancer immune cells and cytolytic activity. These results suggest that high expression levels of BMP7 as well as BMP6 possess higher anticancer immunity which results in better prognosis in ER+ breast cancer. In conclusion, BMP expression profiles may be useful for prognostication. Intervention in this pathway may serve to improve outcomes, manage metastatic disease and assist in clinical decision making on optimal therapy based on the risk of recurrence or metastasis.

Introduction

Bone morphogenetic proteins (BMPs) are members of the TGF-β superfamily and are known for their diverse function in cancer progression (1). BMPs were originally believed to play a role solely in the development of bone and cartilage as they were isolated from bone extracts (2); however, the various functions of BMPs have since been discovered. BMPs display aberrant expression in a variety of malignancies (3). In breast cancer, BMP4 and BMP7 are reported to be the most frequently overexpressed BMPs (4). However, their biological functions in breast cancer remain controversial. Some studies have reported that BMPs suppress breast cancer growth, whereas others have reported that BMPs promote it. BMP4 was reported to suppress breast cancer metastasis by upregulating immune responses and inhibition of myeloid-derived suppressor cells (5), whereas it was also reported that high expression of BMP4 promoted cancer cell migration, invasion and metastasis (6). BMP7 was also shown to both facilitate and inhibit cancer growth, migration, invasiveness and metastasis in different breast cancer cell lines (7). BMP receptors (BMPRs) have also reported to play a role in the carcinogenesis and progression of breast cancer (8). As such, there are multiple conflicting reports concerning the role of BMPs in breast cancer.

TGF-β regulates immune cell infiltration in the tumor microenvironment. TGF-β suppresses anticancer immune cell infiltration such as cytotoxic T cells and natural killer cells, whereas it promotes the function of pro-cancer immune cells, such as regulatory T cells (Tregs) and M2 macrophages (9). BMP6 was also found to induce hormone therapy resistance through macrophage-derived cytokines in prostate cancer (10). Thus, we hypothesized that BMPs regulate immune cell infiltration, resulting in differential survival of breast cancer patients based on the expression levels of BMPs.

Here, we investigated the association of BMP expression and prognosis in breast cancer utilizing The Cancer Genome Atlas (TCGA) (11,12). Furthermore, we explored the difference in the association between BMP expression and survival by estrogen receptor (ER) status and sought possible underlying mechanisms using a bioinformatical approach.

Materials and methods

Data acquisition

Clinical data and gene expression from RNA sequences in TCGA breast cancer cohort were downloaded through cBioPortal (13,14) and UCSC Genome Browser (http://genome.ucsc.edu/) and processed, as previously described (15–20). As a validation cohort, we utilized GSE1456, in which there are 62 Luminal A/B patients with gene expression and survival information from Gene Enrichment Omnibus (21). TCGA and GSE1456 are de-identified and publicly available cohorts, thus institutional review board approval was waived.

Gene set enrichment analysis (GSEA)

GSEA was conducted comparing BMP7 high and low as well as BMP6 high and low in both ER-positive (ER+) and ER-negative (ER−) cohorts with the same cutoffs of survival analyses among 50 hallmark sets (23) using software provided by the Broad Institute (http://software.broadinstitute.org/gsea/index.jsp) as described previously (24).

CIBERSORT

Tumor infiltrating immune cell composition was estimated from the gene expression profiles by CIBERSORT algorism (25). Immune cell fraction data were downloaded through TCIA (https://tcia.at/home) (26).

Cytolytic activity

Cytolytic activity was estimated utilizing GZMA and PRF1 expression (27) as described previously (28).

Statistical analysis

Overall survival (OS) was estimated by Kaplan-Meier curve with log-rank test, comparing high and low expression of each gene of interest. The cutoffs were determined by analyzing multiple cutoff points within the range of values and the optimal cutoff points were chosen based on the statistical significance of the varying points as previously reported (22).

Statistical comparisons of gene expression were performed using the Student t-test, and the clinical demographics were compared using the Fisher exact test. Pearson correlations were calculated based on the expression levels of the genes of interest for correlation matrix (Fig. 2C). All statistical analyses were performed using R software (http://www.r-project.org/) and Bioconductor (http://bioconductor.org/).

Figure 2.

Comparison of the expression of BMPs with ER status. (A) Tumor-suppressive BMPs in the whole cohort. (B) Oncogenic BMPs in the whole cohort. (C) Correlation of BMPs and BMPRs in the ER+ and ER− tumors. BMP, bone morphogenetic protein; BMPR, BMP receptor; ER, estrogen receptor; ER+, ER-positive; ER−, ER-negative.

Results

Various BMPs and BMPRs are associated with better OS, whereas others are associated with worse OS in breast cancer

There have been conflicting reports on the association between the expression levels of BMPs and patient survival (29). Thus, we analyzed the association between the expression of each BMP and OS of the breast cancer patients. High expression of BMP1 (P<0.001), BMP3 (P=0.002), BMP5 (P=0.002), BMP7 (P<0.001) and BMPR1A (P<0.001) was indicative of a significantly better OS of the breast cancer patients (Fig. 1A). Conversely, high expression of BMP6 (P<0.001), BMP8A (P=0.031), BMP8B (P<0.001) and BMPR1B (P=0.004) was indicative of a worse OS (Fig. 1C). These findings indicate that the former may act as tumor suppressors, and the latter as oncogenic factors. On the other hand, the expression levels of BMP2 (P=0.083), BMP4 (P=0.071), BMP10 (P=0.441), BMP13 (P=0.112) and BMP15 (P=0.050) were not associated with OS of breast cancer patients (Fig. S1). Of interest, cutoff points of BMP expression with the highest impact on OS were extremely low in some tumor-suppressive BMPs, and hence a poorer prognosis was shown in <10% of patients for BMP1, BMP7 and BMPR1A (Fig. 1A). Similar trends were shown in some oncogenic BMPs, BMP6 and BMPR1B, in which cutoff points with the highest impact of OS were extremely high and <10% of patients were classified as high expression with worse prognosis (Fig. 1C). In summary, we found that extremely low expression of tumor-suppressive BMPs and extremely high expression of oncogenic BMPs were associated with worse OS in breast cancer.

Figure 1.

Association of BMP expression with OS and expression comparison of tumor and adjacent normal tissues in TCGA breast cancer patients. (A) OS in regards to expression of tumor-suppressive BMPs. (B) Gene expression comparison of tumor-suppressive BMPs between tumor and adjacent normal tissues in matched pairs of TCGA breast cancer. (C) OS in regards to expression of oncogenic BMPs. (D) Gene expression comparison of oncogenic BMPs between tumor and adjacent normal tissues in matched pairs of TCGA breast cancer. BMP, bone morphogenetic protein; TCGA, The Cancer Genome Atlas; OS, overall survival.

Majority of oncogenic BMPs are expressed at a high level, and tumor-suppressive BMPs are expressed at a low in the tumors compared to the normal tissues

Given the associations between the expression levels of the BMPs and worse or better survival, we hypothesized that oncogenic BMPs are expressed at a higher level and tumor-suppressive BMPs are expressed at a lower level in tumors compared to adjacent normal tissues. A total of 114 breast cancers with matched tumor and adjacent normal tissues were present in TCGA, where the expression of BMPs and BMPRs was compared. All of the BMPs and BMPRs, which were associated with survival, were expressed at significantly different levels between the normal and tumor tissues. All of the tumor-suppressive BMPs except BMP1; i.e. BMP3, BMP5, BMP7 and BMPR1A, were expressed at a significantly lower level in the tumors, whereas all of the oncogenic BMPs except BMP6; i.e. BMP8A, BMP8B and BMPR1B, were expressed at a higher level in the tumors compared to the normal tissues (Fig. 1B and D). The expression levels of BMP10 and BMP15, which showed no association with OS, were not detectable in the vast majority of samples (106 normal and 108 tumor tissues, and 104 normal and 100 tumor tissues, respectively). BMP2 and BMP4, that did not show an association with OS, were also expressed at a lower level in the tumors (Fig. S1). These results support the notion that oncogenic BMPs are expressed higher, and tumor-suppressive BMPs are expressed lower in the tumors, except BMP1 and BMP6.

Expression of BMPs differs according to ER status

It is well known that certain signaling, such as HER2, crosstalks with ER signaling and the targeted drug sensitivity is enhanced by ER blockade (30). Signaling of BMPs is also known to crosstalk with that of ER (31,32). Therefore, it was of interest to investigate the association of ER positivity and BMP expression. We found that among the tumor-suppressive BMPs, BMP1 (P=0.002), BMP7 (P<0.001) and BMPR1A (P=0.010) were expressed significantly higher in the ER− tumors compared to the ER+ tumors, whereas there was no significant difference in the expression levels of BMP3 (P=0.247) and BMP5 (P=0.123) (Fig. 2A). On the other hand, among the oncogenic BMPs, BMP6 expression was significantly higher in the ER− tumors (P<0.001), whereas BMPR1B expression was significantly higher in the ER+ tumors (P<0.001) (Fig. 2B). Interactions of BMPs and BMPRs were also different between ER+ and ER− tumors (Fig. 2C). In the ER+ tumors, BMP6 was in the same cluster as BMP1, BMP2 and BMP4, where BMP6 exhibited a higher correlation with BMP1 and BMP2, and BMP7 was in the same cluster of BMPRs (BMPR1A and BMPR1B) (Fig. 2C). Whereas in the ER− tumors, BMP6 and BMP7 were in the same cluster with BMP2 and BMPR1A, showing a very weak correlation with each other (Fig. 2C).

Higher BMP7 expression is associated with better OS in both ER+ and ER− tumors

Since BMP7 was reported to have roles in both promotion and inhibition of cancer growth (4), and showed significant differential expression as well as differential interaction with other BMPs and BMPRs between ER+ and ER− tumors, we focused on BMP7 among the tumor-suppressive BMPs. Higher expression of BMP7 demonstrated better survival in both ER+ (P<0.001) and ER− (P<0.001) cohorts (Fig. 3A and B). Although higher expression of BMP7 was found to be associated with better prognosis in the both ER+ and ER− tumors, GSEA demonstrated that different gene sets were enriched between the ER+ and ER− tumors (Tables SI–SIII). In the ER+ tumors, immune response-related gene sets, such as inflammatory response (P<0.001), IL2/STAT5 signaling (P<0.001), TNFα signaling via NF-κB (P<0.001), IL6/JAK/STAT3 signaling (P<0.001) and allograft rejection (P=0.003), were enriched in the BMP7 high tumors (Fig. 3C, Table SI). As tumor aggressiveness feature-related gene sets, the epithelial-mesenchymal transition (EMT) gene set was also enriched (P<0.001), while cell proliferation, cell cycle-related gene sets were not enriched in the BMP7 high ER+ tumors (Table SI). On the other hand, no gene set was significantly enriched in the BMP7 low tumors. When we compared demographics of the patients with BMP7 high and low ER+ tumors, the patients with BMP7 high ER+ tumors were younger than those with BMP7 low ER+ tumors (Table I). In the ER− tumors, no immune reaction-related gene set was enriched neither in BMP7 high nor low tumors (Tables SII and SIII). These findings suggest that BMP7 plays a tumor-suppressive role in the ER+ tumor by evoking anticancer immunity, which has a more predominant effect than promotion of EMT, whereas it has different mechanisms in the ER− tumors.

Figure 3.

Association of BMP7 expression with OS of breast cancer patients by ER status. OS by BMP7 expression in (A) ER+ and (B) ER- tumors. (C) Gene sets enriched in the BMP7 high ER+ tumors. BMP, bone morphogenetic protein; OS, overall survival; ER, estrogen receptor; ER+, ER-positive; ER−, ER-negative.

Table I.

Patient demographics comparing BMP7 high and low as well as BMP6 high and low expression in the ER+ breast cancer patients.

	BMP7 expression			BMP6 expression
Characteristics	High (n=770) n (%)	Low (n=34) n (%)	P-value	High (n=608) n (%)	Low (n=196) n (%)	P-value
Age, years			0.001[a]			0.033[a]
<60	395 (51.3)	8 (23.5)		318 (52.7)	85 (43.4)
≥60	375 (48.7)	26 (76.5)		290 (47.7)	111 (56.6)
Race			0.658			0.167
Asian	36 (5.2)	2 (6.7)		26 (4.6)	12 (7.5)
African-American	106 (15.9)	3 (10.0)		81 (14.4)	28 (17.6)
Caucasian	548 (79.4)	25 (83.3)		454 (80.9)	119 (74.8)
Menopause status			0.056			<0.001[a]
Pre	160 (24.4)	3 (9.1)		142 (27.1)	21 (12.7)
Post	497 (75.6)	30 (90.9)		382 (72.9)	145 (87.3)
pT			0.343			0.663
T1/2	640 (83.1)	25 (75.8)		501 (82.4)	164 (84.1)
T2/3	130 (16.9)	8 (24.2)		107 (17.6)	31 (15.9)
pN			0.286			0.559
N0	337 (44.6)	18 (54.5)		266 (44.4)	89 (46.8)
N1/2/3	419 (55.4)	15 (45.5)		333 (55.6)	101 (53.2)
pM			0.474			0.084
M0	628 (97.7)	25 (96.2)		494 (98.2)	159 (95.8)
M1	15 (2.3)	1 (3.8)		9 (1.8)	7 (4.2)
Stage			0.842			0.572
Stage I/II	557 (73.9)	24 (81.8)		437 (73.2)	144 (75.4)
Stage III/IV	198 (26.1)	9 (18.2)		160 (26.8)	47 (24.6)
Histology			0.475			<0.001[a]
Infiltrating ductal carcinoma	523 (73.9)	18 (81.8)		381 (68.7)	160 (93.0)
Infiltrating lobular carcinoma	185 (26.1)	4 (18.2)		177 (31.7)	12 (7.0)

ER+, estrogen receptor-positive; BMP, bone morphogenetic protein.

P<0.05, significant difference.

Higher BMP6 expression is associated with better OS in ER+ tumors, whereas it is associated with worse OS in ER− tumors

Since high BMP6 expression was associated with significantly worse survival with lower expression in the tumor tissues, and different expression of BMP6 as well as different correlation with other BMPs and BMPRs between ER+ and ER− tumors were observed, we focused on BMP6 among the oncogenic BMPs. Despite the fact that BMP6 high tumors demonstrated worse OS in the whole cohort, that trend was only apparent in the ER− tumors (P=0.006) (Fig. 4B). Conversely, high expression of BMP6 was found to be associated with better survival in the ER+ tumors (P=0.004) (Fig. 4A). GSEA revealed that immune reaction-related gene sets, including inflammatory response (P<0.001), IL2/STAT5 signaling (P<0.001), TNFα signaling via NF-κB (P<0.001), IL6/JAK/STAT3 signaling (P<0.001) and allograft rejection (P<0.001) were enriched in the BMP6 high ER+ tumors (Fig. 4C and Table SIV). Similar to the BMP7 high ER+ tumors, the EMT gene set was also enriched (P<0.001) in the BMP6 high ER+ tumors (Table SIV). Interestingly, cell proliferation and cell cycle-related gene sets, such as E2F targets (P=0.008), G2/M checkpoint (P=0.019) and MYC targets v1 (33) (P=0.031), were enriched in the BMP6 low ER+ tumors, which suggest that a different mechanism may be in place than BMP7 (Fig. 4D and Table SV). The patients with BMP6 high ER+ tumors showed a higher proportion of younger age, pre-menopausal status, and infiltrating lobular carcinoma in histology compared to those with BMP6 low ER+ tumors (Table I). In the ER− tumors, MYC targets v2 (34) gene set was enriched in the BMP6 high tumors (P=0.041), whereas no immune reaction or cell proliferation gene set was enriched (Tables SVI and SVII). These findings suggest that the BMP6 high tumors have more anticancer immunity and the BMP6 low tumors have accelerated cell cycle; thus, high expression of BMP6 is associated with better survival in the ER+ tumors, whereas higher expression of MYC target genes may play a role in the BMP6 high ER− tumors, resulting in worse survival.

Figure 4.

Association of BMP6 expression with OS of breast cancer patients by ER status. OS by BMP6 expression in (A) ER+ and (B) ER− tumors. Gene sets enriched in (C) BMP6 high ER+ tumors, and (D) BMP6 low ER+ tumors. BMP, bone morphogenetic protein; OS, overall survival; ER, estrogen receptor; ER+, ER-positive; ER−, ER-negative.

BMP7 high as well as BMP6 high ER+ tumors have higher anticancer immune cell infiltration and cytolytic activity

The results that both BMP7 and BMP6 function as tumor suppressors in ER+ tumors were validated in another cohort, GSE1456 (21), in which the patients with high expression of BMP7 (P=0.005) as well as BMP6 (P=0.002) showed significantly better OS in Luminal A/B tumors (Fig. 5). Given these results and the BMP7 and BMP6 high expressing tumor enriched immune-related gene sets, we further hypothesized that the ER+ tumors with high BMP7 as well as high BMP6 expressions have higher anticancer immune cell infiltration. Utilizing CIBERSORT algorithm (25), we found higher anticancer immune cell infiltration, such as gamma-delta (γδ) T cells (P=0.029), M1 macrophages (P<0.001), and less pro-cancer immune cell infiltrations, such as M2 macrophages (P<0.001) and Tregs (P=0.020) in the BMP7 high ER+ tumors (Figs. 6A and S2). A similar trend was found in the BMP6 high ER+ tumors, where higher infiltration of anticancer immune cells, not only γδ T cells (P<0.001) and M1 macrophages (P<0.001), but also CD8+ T cells (P<0.001), CD4+ naïve T cells (P=0.008) and CD4+ activated memory T cells (P=0.003), and less infiltration of pro-cancer immune cells, M2 macrophages (P<0.001) and Tregs (P<0.001) in the BMP6 high ER+ tumors were observed (Figs. 6B and S3). Since higher anticancer immune cell infiltration was observed in both BMP7 and BMP6 high expressing ER+ tumors, we further analyzed the cytolytic activity of those tumors. As shown in Fig. 6C and D, high BMP7 as well as high BMP6 in the ER+ tumors demonstrated higher cytolytic activity compared to those with low expression tumors. These findings imply that ER+ tumors with high BMP7 as well as high BMP6 expression have not only higher immune cell infiltration, but also higher cytolytic activity, resulting in better prognosis.

Figure 5.

Validation of survival analyses utilizing Luminal A/B tumors in GSE1456. (A) OS by BMP7 expression. (B) OS by BMP6 expression. BMP, bone morphogenetic protein; OS, overall survival.

Figure 6.

Immune activity by BMP7 and BMP6 expression in the ER+ tumors. (A) Infiltrating immune cell comparison between the BMP7 high and low ER+ tumors. (B) Infiltrating immune cell comparison between the BMP6 high and low ER+ tumors. (C) Cytolytic activity, PD1 and PD-L1 expression between the BMP7 high and low ER+ tumors. (D) Cytolytic activity, PD1 and PD-L1 expression between the BMP6 high and low ER+ tumors. BMP, bone morphogenetic protein; ER+, ER-positive; PD1, programmed cell death protein 1; PD-L1, programmed death-ligand 1.

To ascertain whether these tumors are good candidates for immune check point inhibition, PD1 and PD-L1 expressions were compared between BMP7 high and low as well as BMP6 high and low expression in the ER+ tumors. PD1 and PD-L1 expression levels were significantly higher in the BMP7 high expressing as well as BMP6 high expressing ER+ tumors (Fig. 6C and D). These findings suggest that high expression of BMP7 and BMP6 in the ER+ tumors may be good candidates for immune check point inhibition.

Discussion

In the present study, we found that each BMP was associated differentially with OS of breast cancer patients. High expression levels of BMP1, BMP3, BMP5, BMP7 and BMPR1A were associated with better prognosis in the whole cohort of breast cancer patients, whereas those of BMP6, BMP8A, BMP8B and BMPR1B were associated with worse prognosis. BMP7 was associated with better prognosis in both the ER+ and ER− tumors. In contrast, higher BMP6 expression was found to be associated with better prognosis in the ER+ tumors, but was associated with worse prognosis in the ER− tumors. High expression of BMP7 as well as BMP6 was associated with higher immune response in the ER+ tumors, but not in the ER− tumors.

BMPs play complex roles in cancer progression, specifically in breast cancer (34). Our findings, that high expression levels of BMP7 and BMP6 are associated with better prognosis in ER+ breast cancer, is consistent with previous reports. Utilizing breast cancer cell lines, Takahashi et al demonstrated that BMP6 and BMP7 inhibit estrogen-induced proliferation (35). It was reported that higher BMP7 expression is associated with higher rates of bone metastasis, however, lower BMP7 expression is associated with local relapse in breast cancer (36). In our analysis BMP7 expression was found to be associated with a better overall survival, irrespective of ER status. In the same token, BMP6 expression has been demonstrated to suppress breast cancer metastasis in both ER+ and ER− tumors (37,38). It was also found that BMP6 silencing increased chemo-resistance to doxorubicin by upregulating multiple drug resistance glycoproteins in breast cancer (39). In our analysis, BMP6 expression was found to be associated with worse survival in the whole patient cohort. Notably, the reverse finding was observed in the ER+ tumors, where BMP6 plays a tumor-suppressive role. This suggests a complex interplay with the hormonal pathway.

There have been few reports that link BMP7 or BMP6 to immune response. During the wound healing process, BMP7 plays an anti-fibrotic role by causing inflammation (40). BMP6 expression was found to be elevated in the ovary of polycystic ovarian syndrome which is known to be associated with inflammation, when compared to normal ovarian tissues (41). Overall, the role of BMPs in breast cancer remains poorly understood and is associated with complex interactions with a multitude of pathways and effectors that may act to completely alter the effector function of BMP signaling.

Although novel findings were discovered, there are limitations to the present study. First, our results were confirmed in a few cohorts due to availability. Second, all of our findings were based on gene expression of the primary tumor, thus experimental approaches are needed to elucidate further molecular mechanisms.

In the context of ER expression, individualized patient BMP expression profiles may act as a prognostic tool to guide therapy or may even be used as a therapeutic target in itself. This revelation in findings will continue to evolve, and the correlation between the BMP activation in regards to different pathways and the effect on prognosis in the clinical realm remain to be clearly elucidated.