Integrated RNA-seq and DNase-seq analyses identify phenotype-specific BMP4 signaling in breast cancer

RESEARCH ARTICLE Open Access Integrated RNA-seq and DNase-seq analyses identify phenotype-specific BMP4 signaling in breast cancer M. Ampuja 1,2*† , T. Rantapero 1† , A. Rodriguez-Martinez 1,2† , M. Palmroth 1 , E. L. Alarmo 1 , M. Nykter 1 and A. Kallioniemi 1,2 Abstract Background: Bone morphogenetic protein 4 (BMP4) plays an important role in cancer pathogenesis. In breast cancer, it reduces proliferation and increases migration in a cell line-dependent manner. To characterize the transcriptional mediators of these phenotypes, we performed RNA-seq and DNase-seq analyses after BMP4 treatment in MDA-MB-231 and T-47D breast cancer cells that respond to BMP4 with enhanced migration and decreased cell growth, respectively. Results: The RNA-seq data revealed gene expression changes that were consistent with the in vitro phenotypes of the cell lines, particularly in MDA-MB-231, where migration-related processes were enriched. These results were confirmed when enrichment of BMP4-induced open chromatin regions was analyzed. Interestingly, the chromatin in transcription start sites of differentially expressed genes was already open in unstimulated cells, thus enabling rapid recruitment of transcription factors to the promoters as a response to stimulation. Further analysis and functional validation identified MBD2, CBFB, and HIF1A as downstream regulators of BMP4 signaling. Silencing of these transcription factors revealed that MBD2 was a consistent activator of target genes in both cell lines, CBFB an activator in cells with reduced proliferation phenotype, and HIF1A a repressor in cells with induced migration phenotype. Conclusions: Integrating RNA-seq and DNase-seq data showed that the phenotypic responses to BMP4 in breast cancer cell lines are reflected in transcriptomic and chromatin levels. We identified and experimentally validated downstream regulators of BMP4 signaling that relate to the different in vitro phenotypes and thus demonstrate that the downstream BMP4 response is regulated in a cell type-specific manner. Keywords: Bone morphogenetic protein, Breast cancer, NGS, RNA-seq, DNase-seq, Transcription factor Background Despite many advances in diagnostics and therapeutics, breast cancer remains the leading cause of cancer death in women [1]. Bone morphogenetic proteins (BMPs) are a group of growth factors that are important players during development [2, 3] but also contribute to cancer formation and progression [4–6]. As a subfamily of the transforming growth factor β(TGF-β) protein superfamily, BMPs are extracellular ligands that bind as dimers to their specific transmembrane receptors and activate the intracellular SMAD signaling pathway leading to phosphorylation of receptor-regulated SMADs (SMAD1/5/9). The activated SMADs bind to SMAD4 and the complex translocates to the nucleus where it regulates the expression of BMP target genes [7, 8]. Alternatively, BMP signals are also mediated through the activation of ERK, JNK and p38 mitogen-activated protein kinase pathways [7, 8]. The functional consequences of BMP signaling depend on the BMP ligand and tissue type. We and others have shown that BMP4 reduces the proliferation of breast cancer cell lines, while simultaneously inducing migration and invasion in a subset of cell lines [9–11]. Similar * Correspondence: [email protected] † Equal contributors 1 BioMediTech, University of Tampere, Tampere, Finland 2 Fimlab Laboratories, Tampere, Finland © The Author(s). 2017 Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated. Ampuja et al. BMC Genomics (2017) 18:68 DOI 10.1186/s12864-016-3428-1 dualistic effects upon BMP4 stimulation have also been reported in other tumor types [12]. Concordantly, data from breast cancer patient samples point to a correlation between elevated BMP4 levels and reduced proliferation as well as an increased risk of recurrence [13]. These BMP4-related effects that seem either detrimental (reduced cell growth) or beneficial (increased mobility) for the cancer cells are likely to be mediated by specific BMP4 target genes. The identification of such target genes is thus important since it may allow generation of effective cancer therapies targeting each phenotype independently. We have previously searched for BMP4 target genes in a set of breast cancer cell lines that predominantly respond to BMP4 treatment by reduction of proliferation [14]. Here, we used next-generation sequencing (NGS) technologies (RNA-seq and DNase-seq) to uncover BMP4-mediated transcriptional events with a specific focus on comparing cells in which BMP4 has opposing effects, namely antiproliferative and promigratory. Out of the nine breast cancer cell lines we have previously studied, T-47D shows one of the most prominent growth reductions and MDA-MB-231 cells display the most overt induction of migration [9, 10], and were thus selected for this study. RNA-seqmethodquantifiesthelevelofgeneexpression across the genome [15] while DNase-seq allows identification of open chromatin regions that are sensitive to digestion by the DNase I endonuclease [16]. Open chromatin regions are considered as sites where transcriptional regulation can take place since they are accessible for regulatory molecules to bind and exert their function. By combining data from RNA-seq and DNase-seq, and using additional data analysis tools, it was possible to identify candidate transcription factors involved in the observed transcriptional responses. This approach thus provides the means to better understand the transcriptional events that link BMP4 signaling and its resulting phenotypes. Results We performed RNA-seq and DNase-seq analyses in two breast cancer cell lines, T-47D and MDA-MB-231. The cell lines were treated with BMP4 and vehicle control for 3 h, thus allowing us to specifically focus on early response events. Both vehicleand BMP4-treated cell lines were sequenced (see methods). BMP4-elicited transcriptional regulation is highly divergent in the two breast cancer cell lines with different functional responses to BMP4 Sequencing reads from RNA-seq and DNase-seq were aligned to the human genome and further analyzed as described in the methods. To confirm that the two datasets were consistent, we compared the chromatin openness as determined by DNase-seq signal at the transcription start site (TSS) to the expression level of the gene as determined by RNA-seq. As expected, we found that the increased openness of TSS globally correlated with increased gene expression (Additional file 1: Figure S1, Panels A and B). However, the variance is high, indicating that the differences in the chromatin state only partly explain gene expression patterns. Next we compared the expression levels from RNA-seq between the vehicleand BMP4-treated cells. This analysis identified 91 differentially expressed genes (DEGs) in MDA-MB-231, of which 58 were upregulated and 33 downregulated (Additional file 2: Table S1). In T-47D, there were 203 DEGs, of which 160 were upregulated and 43 were downregulated (Additional file 3: Table S2). In total, 10 DEGs (ATOH8,BDKRB2,BMF,GS1-124 K5.4, ID1,ID2,ID3,SKIL,SMAD6,andSMAD9)wereshared by the two cell lines and all of them were upregulated except GS1-124 K5.4 which was downregulated in both cell lines. To illustrate that BMP4 induces markedly divergent transcriptional responses in these two cell lines, we generated a heatmap to show the expression levels of the protein-coding DEGs (Fig. 1a). Using the DNase-seq data, we examined the chromatin status at the transcription start sites (TSSs) of these protein-coding DEGs. For the majority of the cases the chromatin was open at the TSS before BMP4 stimulation (approximately 86% of all DEGs in both cell lines) (Additional file 1: Figure S1, Panels C and D). For the remaining DEGs, we observed either opening or closing of the TSS after stimulation or no change in the closed chromatin status (Fig. 1a). These data indicate that, at this early time point, the BMP4-induced differential expression mainly involves genes whose transcription does not require changes in the chromatin status at TSS. The DEG lists included a number of genes involved in the canonical BMP pathway. As expected, ID1, ID2 and ID3, known BMP4 target genes, were upregulated in both cell lines (Fig. 1b). Similarly, the receptorregulated SMAD9 was upregulated in both cell lines whereas no significant difference in the other receptorregulated SMADs or SMAD4 expression was observed. Among the inhibitory SMADs, SMAD6wasupregulated in both cell lines and SMAD7 in T-47D. In addition, the BMP type I receptor BMPR1A and negative regulators of BMP signaling, NOG and BAMBI, were upregulated in T-47D while in MDA-MB-231 their expression was not significantly changed (Fig. 1b). Thus BMP4 stimulation leads to expression changes having characteristics of both feedback and feedforward loops. We then evaluated whether the differentially expressed genes participate in specific biological processes and Ampuja et al. BMC Genomics (2017) 18:68 Page 2 of 15 A B Fig. 1 The RNA-seq and DNAse-seq data reveal cell line-specific responses to BMP4. aGene expression levels of differentially expressed protein-coding genes converted to log2 scale are shown for both cell lines and treatments, upregulated genes on the left and downregulated genes on the right. The status column denotes the cell line in which the gene is differentially expressed. The rightmost columns indicate the status of the chromatin at transcription start sites (TSS) of the DEGs as measured by DNAse-seq. bIllustration of the differentially expressed components of the BMP signaling pathway upon BMP4 treatment Ampuja et al. BMC Genomics (2017) 18:68 Page 3 of 15 especially assessed whether the non-common DEGs have differing functions. To this end we used DAVID to search for GO terms enriched in the sets of noncommon protein-coding DEGs. In MDA-MB-231, most of the enriched terms were related to cell migration whereas organ development and morphogenesis as well as intracellular signaling were the most significant GO terms in T-47D (Table 1). These findings imply that the transcriptional changes are indeed likely to explain the dissimilarities in the phenotypic responses of these two cell lines to BMP4 treatment. Thereafter, we also wanted to investigate whether the expression levels of DEGs could be linked with survival in breast cancer patients. For this purpose, we used the data publicly available in the TCGA database. The results showed that 20 DEGs in the MDA-MB-231 and 46 DEGs in the T-47D cells associated with either good or poor prognosis (Additional file 4: Tables S3 and S4). Of the nine shared protein-coding DEGs, four (ATOH8, ID3,SMAD6 and SMAD9) were correlated with survival, all being associated with poor prognosis. To validate the results of the RNA-seq analysis and to extend the scope of the study beyond the 3 h time point in two cell lines, qRT-PCR was used to study the expression levels of 15 selected DEGs in MDA-MB-231 and T-47D cells as well as in five additional breast cancer cell lines (BT-474, HCC-1954, MCF-7, MDAMB-361, and MDA-MB-436) and one normal breast epithelial cell line (MCF-10A) treated with BMP4 and vehicle for 3, 6 and 24 h. The genes were selected based on their expression levels and reported cancer association in the literature, and five of these were upregulated according to the RNA-seq in both MDA-MB-231 and T-47D. The expression patterns of the majority of the genes showed similarities across the cell line panel and time points with the clear exception of MDA-MB436, in which the expression changes were very limited (Fig. 2). Particularly the five shared genes (ATOH8,ID2, SKIL,SMAD6 and SMAD9)aswellasDLX3 were consistently upregulated upon BMP4 treatment throughout the time series thus confirming that they represent common BMP4 target genes. The remaining genes showed more variability with altered expression typically in only two to three cell lines, suggesting that their expression is likely to be influenced by factors that are cell line-specific. Chromatin landscape and dynamics following BMP4 treatment To gain more insight into the changes of chromatin structure during BMP4 treatment, we performed peak detection in a genome-wide manner to identify the areas of open chromatin. The peak detection approach was benchmarked by comparison to publicly available DNase-seq data of unstimulated T-47D cell line from ENCODE (see methods), showing that most of the peaks identified in our data are present also in ENCODE samples (Additional file 5: Table S5). After filtering procedures (see methods), the numbers of identified DNase hypersensitive sites (DHSs) in the MDA-MB-231 cell line were 89,830 and 97,349 in vehicleand BMP4-treated samples, respectively. In T47D, the corresponding numbers were 68,000 and 73,881. To obtain a unified set of peaks for both conditions, the overlapping DHSs were merged resulting in a total of 106,154 DHSs in MDA-MB-231 and 110,028 in T-47D. After the merging, the fraction of shared DHSs between BMP4 and vehicle control in MDA-MB-231 samples was 75% while the fraction of unique DHSs in the vehicle was 9% and correspondingly in the BMP4 sample 16% (Additional file 6: Figure S2). In the T-47D cell line, Table 1 Gene ontology analysis Cell line GO accession GO term Number of genes Adjusted p-value MDA-MB-231 GO:0030334 regulation of cell migration 5 2.0 × 10 −2 GO:0030335 positive regulation of cell migration 4 2.3 × 10 −2 GO:2000145 regulation of cell motility 5 2.4 × 10 −2 GO:2000147 positive regulation of cell motility 4 2.5 × 10 −2 GO:0051272 positive regulation of cellular component movement 4 2.7 × 10 −2 T-47D GO:0048513 animal organ development 45 2.6 × 10 −8 GO:0035556 intracellular signal transduction 41 4.5 × 10 −8 GO:0009887 organ morphogenesis 22 4.0 × 10 −7 GO:0009966 regulation of signal transduction 36 4.5 × 10 −6 GO:0007166 cell surface receptor signaling pathway 34 9.5 × 10 −5 The DAVID Functional Annotation Tools was used to reveal significantly enriched GO categories among the differentially expressed protein-coding genes. The analysis was done independently for each cell line and shared differentially expressed genes were omitted. The top five biological function GO terms are shown Ampuja et al. BMC Genomics (2017) 18:68 Page 4 of 15 the fraction of shared DHSs between the two conditionswas27%whereasthefractionofuniqueDHSsin the vehicle was 34% and in the BMP4 sample 39% (Additional file 6: Figure S2). Annotation of the merged DHSs to genomic features revealed a similar distribution in the two cell lines in the vehicle-treated condition, with the largest fraction (>30%) of DHSs locating in introns (Fig. 3a). When comparing the distributions of the BMP4-induced DHSs between the cell lines apparent resemblances were also observed. In both cell lines, the proportion of DHSs associated with intronic and intergenic regions increased after BMP4 stimulation with a corresponding decrease at other genomic locations, including the promoter regions (Fig. 3b). Fig. 2 Expression levels of selected BMP4 target genes by qRT-PCR in a breast cancer cell line panel. The expression levels of 15 DEGs were measured after 3, 6 and 24 h of BMP4 treatment in the indicated cell lines. The color code illustrates the relative expression levels in the BMP4-treated sample as compared to the corresponding vehicle control. FC = Fold change, n.a. = mRNA level too low to allow reliable measurement A B Fig. 3 Distribution of open chromatin regions. Annotation of open chromatin regions in MDA-MB-231 and T-47D after (a) vehicle treatment (basal openness) and (b) BMP4 treatment (consisting only of the chromatin that opened after BMP4 treatment) Ampuja et al. BMC Genomics (2017) 18:68 Page 5 of 15 To assess the functional impact of the BMP4-induced global changes in the chromatin structure we conducted an enrichment analysis using GREAT [17] which maps the DHSs to putative regulatory regions of genes and conducts a gene ontology enrichment analysis. The results highlighted e.g. cell motility and organ morphogenesis as enriched biological functions for MDA-MB-231 and T-47D, respectively (Additional file 7: Tables S6 and S7). These results are consistent with those obtained by enrichment analysis of the differentially expressed genes from RNA-seq (Table 1) and thereby suggest that, together with specific target genes, BMP4-induced changes at chromatin level may contribute to the emergence of the different BMP4-mediated phenotypes. Transcription factor binding site enrichment analysis in open chromatin regions of promoters reveals transcription factors involved in BMP4 signaling regulation Based on our TSS openness analysis (Fig. 1a), a dominant feature of our data is that the chromatin of the putative BMP4 target genes (identified by RNA-seq) is open already in vehicle-treated cells. This is further supported by our genome-wide peak analysis, where the promoter regions were not overrepresented after the treatment (Fig. 3b). Therefore, the alterations in the chromatin state only partially explain gene expression differences induced by the BMP4 treatment. However, differential transcription factor binding to open promoters may explain the different responses in the cell lines. Therefore we performed transcription factor (TF) motif binding analysis. To assess which TFs might be regulators of the BMP4 response, the sequences of open chromatin sites in the proximal promoters of upregulated genes were analyzed with a total of 426 position weight matrixes (PWMs), representing 401 individual TFs or TF-complexes (see methods). For each TF we calculated an enrichment score (see methods) for the number of binding sites in either MDA-MB-231 or T-47D cells. This analysis led to the identification of candidate regulator TFs, including multiple members of the SMAD family of TFs, as expected, as well as a number of shared common regulator TFs. To focus on biologically relevant candidates, we filtered out those TFs that were not expressed based on our RNA-seq data. In addition, we included only those TFs whose binding sites (TFBSs) in open chromatin regions of the promoters of DEGs were enriched in one and depleted in the other cell line. The top 15 TFs that are expressed in both cell lines but have a high enrichment score only in one of the cell lines are listed in Tables 2 and 3. Examples of target gene promoters with binding motifs for predicted TFs are shown in Fig. 4a. For more in-depth functional analysis we selected particular TFs from the top enriched candidates using the following criteria: 1) a binding motif with a quality category of A-C in the HOCOMOCO database, 2) relevance in the context of our model based on literature, 3) not a highly common regulator or part of a large TF family, and 4) high expression level of the TF (>1000 reads) in at least one cell line and differential expression between cell lines according to the RNA-seq. The Table 2 Top 15 transcription factors enriched in MDA-MB-231 cells TF name Motif Selection by: TF binding sites Ref. sites Expected sites in ref. Ratio of enrichment Mean read count MYBL2 MYBB_f1 2, 3, 4 12 2930 6.2 1.92 2197 BACH1 BACH1_si 1, 2, 3 15 3904 8.3 1.81 531 MYC MYC_f1 1, 2, 4 10 2698 5.7 1.74 3044 MAFK MAFK_si 2, 3 16 4428 9.4 1.70 688 RELA TF65_f2 1, 2, 4 19 5467 11.6 1.63 1398 PPARA PPARA_f1 1, 2, 3 9 2747 5.8 1.54 185 NFIA/B/C/X a 1, 2, 3 15 4669 9.9 1.51 b NFIL3 NFIL3_si 1, 2, 3 11 3494 7.4 1.48 474 FOXA2 FOXA2_f1 1, 2, 3 36 11477 24.4 1.47 434 REL REL_do 1, 2, 3 17 5422 11.5 1.47 69 ZFHX3 ZFHX3_f1 2, 3 46 14683 31.2 1.47 66 RXRB RXRB_f1 1, 2, 4 20 6414 13.6 1.47 1015 SMARCC1 SMRC1_f1 1, 4 20 6443 13.7 1.46 1478 ETV5 ETV5_f1 2, 3 16 5199 11.1 1.45 641 NR3C1 GCR_si 1, 2, 4 15 4910 10.4 1.44 1087 The ratio of enrichment is the result of dividing the number of TF binding sites by the number of expected sites. Motifs are derived from the HOCOMOCO database. a NFIA + NFIB + NFIC + NFIX_f2, b Read count range (51, 148, 748, 444, respectively). Ref. reference Ampuja et al. BMC Genomics (2017) 18:68 Page 6 of 15 last criteria was used to ensure methodological success in subsequent functional assays. With the criteria described above CBFB, HIF1A, and MBD2 were selected for further study. Of these, MBD2 had a large number of binding sites in the promoters of our DEGs while binding sites of the other two TFs were less widespread. Inaddition,SMAD4wasusedasapositivecontrol. As SMAD4 is a known regulator of BMP signaling, we performed co-occurrence analysis of the binding sites between our three candidate TFs and the SMAD motifs. We found that the MBD2 motif was significantly co-localized with the GC-rich SMAD4 consensus motifs CGCC (P= 1.1e-9), GCCGnCGC (P=1.3e-14), and GGCGCC (P= 2e-10). As binding sites for CBFB or HIF1A were less frequent across DEGs, statistical significance for co-localization with SMAD motifs could not be reliably evaluated. However, we did find several promoters where SMAD binding sites co-localized with these factors. Silencing of selected TFs (SMAD4, CBFB, HIF1A, and MBD2) was then used to further evaluate their impact on BMP4 signaling. After 48 h of silencing, the cells were treated with BMP4 for 24 h and the mRNA levels of the validated DEGs were measured to assess whether the silencing influences BMP4 target gene expression (Fig. 4b and Additional file 8: Figure S3). Downregulation of SMAD4 was able to reverse the BMP4-mediated change in the expression of all the tested target genes in both MDA-MB-231 and T-47D cells (Fig. 4c) indicating that these expression changes are indeed transmitted via the canonical BMP pathway. For most of the target genes, MBD2 silencing led to abrogation of the BMP4mediated induction in gene expression in both cell lines. In T-47D cells, similar data was also obtained for most of the genes upon CBFB (9/10) and HIF1A depletion (6/10). However in MDA-MB-231, silencing of HIF1A resulted exclusively in upregulation of the target genes and both enhanced and diminished expression was seen after CBFB downregulation. Of note, silencing of all of the TFs in T-47D cells led to the enhanced expression of the DLL1 gene, which was consistent with it being downregulated upon BMP4 treatment. These data imply that the TFs may function as either repressors or enhancers of BMP4 target gene expression in a context-dependent manner. Discussion We have previously characterized transcriptional responses of breast cancer cell lines to BMP4 by using microarray technology[14].However,inthatstudywefocusedonly on cells that respond to BMP4 by reduced proliferation. Efforts by others to examine BMP signaling target genes have concentrated exclusively on non-cancerous cells [18–20]. Here we set out to uncover the transcriptional responses of breast cancer cell lines with different phenotypes by using one cell line that responds to BMP4 by reduced proliferation (T-47D) and another that reacts with increased migration (MDA-MB-231). Being able to uncover the mechanisms of these two different responses is essential for the understanding of the role of BMP4 in breast cancer pathogenesis. To this end, we used a substantially new approach of combining DNaseseq, RNA-seq and functional experiments. In order to find the early mediators of BMP4 response, we treated the cells with BMP4 or vehicle control for Table 3 Top 15 transcription factors enriched in T-47D cells TF name Motif Selection by: TF binding sites Ref. sites Expected sites in ref. Ratio of enrichment Mean read count MBD2 MBD2_si 1, 2, 3, 4 101 6664 39.6 2.55 571 TFAP2A AP2A_f2 1, 2, 3 115 10363 61.6 1.87 941 E4F1 E4F1_f1 2, 3 18 1750 10.4 1.73 310 SP1 SP1_f1 1, 2 392 41453 246.3 1.59 838 CUX1 CUX1_f1 1, 2, 3 13 1462 8.7 1.50 141 E2F2 E2F2_f1 1, 2 17 1941 11.5 1.47 215 AHR AHR_si 1, 2, 3 9 1030 6.1 1.47 791 SP2 SP2_si 1, 2 140 16512 98.1 1.43 672 CREB1 CREB1_f1 1, 2, 3 23 2720 16.2 1.42 177 CBFB PEBB_f1 1, 2, 3, 4 46 5461 32.4 1.42 457 ZIC2 ZIC2_f1 1, 2, 3 46 5487 32.6 1.41 118 ZFX ZFX_f1 1, 2, 3 127 15650 93.0 1.37 287 HIF1A HIF1A_si 1, 2, 3, 4 15 1890 11.2 1.34 1847 E2F3 E2F3_si 1, 2, 3 16 2019 12.0 1.33 322 XBP1 XBP1_f1 1, 3, 4 12 1545 9.2 1.31 22744 The ratio of enrichment is the result of dividing the number of TF binding sites by the number of expected sites. Motifs are derived from the HOCOMOCO database. Ref. reference Ampuja et al. BMC Genomics (2017) 18:68 Page 7 of 15 3 h. At this time point, the canonical BMP pathway through SMAD1/5/9 is already activated [9]. The results of RNA-seq revealed that the cell lines responded to BMP4 by upregulating or downregulating a set of genes that were mostly cell line-specific, with only ten common DEGs identified. Consistent with the sequencing data, validation with qRT-PCR across multiple time points (3, 6, and 24 h) and five additional cell lines further confirmed in a wider context the existence of common BMP4 target genes as well as cell line-specific expression patterns. Of the ten shared DEGs, three were known BMP4 target genes (ID1-3) and two members of the BMP A BC Fig. 4 Examples of predicted TFBSs and the impact of transcription factors on BMP4 target gene expression. aThe predicted binding sites of transcription factors MBD2, HIF1A and CBFB are depicted at the promoters of NOG,SMAD7 and ID1 genes, respectively. In addition, known BMP-response elements (BRE) located near the binding sites are illustrated. bThe TFs were silenced and the cells were treated with BMP4 or vehicle control followed by measurement of target gene expression by qRT-PCR. Examples of relative expression levels of SKIL after HIF1A silencing in MDA-MB-231 cells (top panel) and NOG expression after MBD2 silencing in T-47D cells (bottom panel) are shown. cGraphical summary of the TF silencing experiments. The order of the genes is identical to that in Fig. 2. Blue color (decreased target gene expression) denotes TFs that were essential for target gene expression and red color those whose silencing led to enhanced target gene expression. Not applicable indicates cases where BMP4 did not alter the baseline gene expression. Data on the DLL1 gene, which is downregulated in T-47D upon BMP4 treatment, are highlighted with a bold line Ampuja et al. BMC Genomics (2017) 18:68 Page 8 of 15 signaling pathway (SMAD6,SMAD9) [3]. The activation of the inhibitory SMAD6 indicates a negative feedback loop, which in T-47D is reinforced by the upregulation of BMP antagonist NOG and the pseudoreceptor BAMBI. On the other hand, activation of the receptor-regulated SMAD9 seems to point to a positive feedback loop, as alongside other R-SMADs, SMAD9 has been found to enhance BMP signaling [21, 22]. However, one study indicated that SMAD9 may have an inhibitory role in BMP signaling [23]. In any case, upregulation of SMAD9 due to BMP4 treatment has also been recently reported in various cell types, for example in primary fibroblasts, hepatocellular carcinoma and melanoma cells [24]. To understand the function of the cell line-specific DEGs, we used GO analysis to segregate the DEGs into biological process categories and discovered that the results reflected the response of the cell lines to BMP4. Processes related to migration were enriched in the MDA-MB-231 cells, whereas more diverse responses were found in T-47D, including categories comprised of signaling, development and morphogenesis. These findings were corroborated by the DNase-seq data, where we found that BMP4-induced global open chromatin sites were enriched with the same biological categories that were found with RNA-seq data. While categories associated with signaling were observed in both cell lines, in MDA-MB-231 those related to migration were enriched. These data extend our previous results showing enrichment of differentially expressed genes in GO categories that were associated with the BMP4-induced decrease in proliferation [14]. Taken together, the different responses of the cell lines to BMP4 are reflected both at the transcriptional and chromatin levels. In the analysis of TSS chromatin state we could observe changes in only a few of the genes that were differentially regulated by BMP4. This might be due to the fact that the 3 h stimulation of BMP4 is too short for most of the TSSs to change their chromatin status. Moreover, we could observe that in many cases the chromatin was already open at the TSS, in which case further changes are not needed to enhance the transcriptional activity. Together with the observation that there is a large variation between the chromatin status and gene expression when we extend the analysis to the whole set of protein-coding genes, it can be concluded that the chromatin state of TSSs explains the observed expression patterns only to a small extent. This result was not unexpected, as gene expression is also commonly regulated from regions located far from the TSS, such as enhancers [25, 26]. With genome-wide detection of open chromatin areas we noticed that BMP4 stimulation induces opening of the chromatin mostly in the intronic and intergenic regions. This is consistent with the fact that changes in the TSS and promoter regions were observed with only a few of the differentially expressed genes. Opening of the intronic sequences may indicate increased level of RNA polymerase activity at gene bodies. Chromatin opening at intergenic regions might suggest that additional regulatory control is being attained in large extent through distal regulatory elements such as enhancers and silencers. Thus, already at the early 3-h time point we are able to observe conformational changes that cells may utilize in more detailed regulation of the BMP response. Unfortunately, based on this analysis we were not able to define a specific transcription factor chromatin signature that could be used to define BMPspecific regulatory sequences. Hence detailed analysis of the putative enhancer regions would require more specific measurement data about the chromatin interactions in these cells. To further characterize the regulation of BMP4 target genes, we analyzed transcription factor binding sites (TFBSs) in the open chromatin regions located on gene promoters. Among the top 15 enriched TFs, there were a few which had previously been linked to BMP target gene regulation. For example, XBP1 and RELA have been shown to be repressors of BMP target genes Xvent-2 and Id1, respectively [27, 28]. Using enrichment of the TFBSs between cell lines as well as other criteria, we selected three TFs (CBFB, HIF1A, and MBD2) for functional characterization and silenced them in the two cell lines. In addition, we used SMAD4, a key component of the canonical BMP pathway, as our positive control and indeed SMAD4 was required for transcriptional regulation of all the BMP4 target genes in the assay. Although BMPs can signal through alternative pathways [7, 8], this result points to regulation through the canonical pathway. In contrast, the response to other transcription factors was more variable and cell line-specific. MBD2 is a methyl-CpG-binding transcription factor that plays a role in development [29, 30]. Several studies have shown that MBD2 acts as a transcriptional repressor by recruiting co-repressor complexes to promoters, which in turn leads to formation of repressive chromatin through chromatin remodelling [31, 32]. However, there is also evidence that MBD2 can activate transcription by removing methylation from CpG islands located in promoters [33]. In both cell lines, MBD2 seemed to act mainly as an activator of transcription, although its role was more prominent in MDA-MB-231. In our analysis, MBD2 had a large number of binding sites across DEGs and it was highly expressed in both cell lines, consistent with the observed behavior in the silencing experiment. The key role of MBD2 in controlling the BMP4 response suggests that DNA methylation may be involved in BMP4 signaling. Ampuja et al. BMC Genomics (2017) 18:68 Page 9 of 15