Activation of the TGFβ1-induced signaling pathway in MCF10A and MCF7 cells
MCF10A nontumorigenic breast epithelial cells (isolated from the mammary gland of a patient with fibrocystic disease) and MCF7 breast adenocarcinoma cells (estrogen receptor-positive, ER+) were selected for the study of TGFβ1 (a key member of the TGFβ superfamily) signaling. Both cell lines are epithelial and are characterized by CDH1 expression (Figure S1A). TGFβ1 treatment induced morphological changes (shape elongation) that were noticeable around the third day of treatment in both cell lines, with MCF7 cells additionally undergoing changes leading to cell death (PARP1 cleavage was observed, indicating caspase activation; Figure S1B, C). TGFβ acts through specific receptors to activate intracellular pathways leading to the phosphorylation of SMAD2/3 proteins. Therefore, to monitor the activation of the pathways triggered by TGFβ1, we assessed the level of SMAD3 phosphorylation at Ser423/425 (Figure S1B, C). Pilot experiments showed its strong phosphorylation and nuclear accumulation after just one hour of TGFβ1 treatment. Only in MCF10A cells was this accompanied by the significant nuclear accumulation of the SNAI1 and SNAI2 transcription factors (which are known to promote CDH1 repression; [32], [33], [34]). In addition, experiments showed that in MCF10A cells, the response to TGFβ1 treatment manifested by SMAD3 phosphorylation was extinguished over the next few hours. However, in MCF7 cells, SMAD3 phosphorylation remained constant for at least 24 hours (until a new dose of TGFβ1 was administered) (Figure S1B, C). These initial experiments showed that noncancerous and cancerous breast epithelial cells (MCF10A and MCF7, respectively) responded differently to TGFβ1 stimulation. Therefore, we performed gene expression profiling to further clarify the differences between the two cell lines. We treated the cells for six consecutive days, with daily administration of TGFβ1 for two hours. The TGFβ1-containing medium was changed to fresh medium, which was collected after another 22 hours (so-called conditioned medium, CM) and could be administered to other cells. This enabled us to assess the ability of extracellularly secreted mediators to induce a response in naïve cells (bystander effect). RNA-seq analyses were performed on cells either directly treated with TGFβ1 or treated with conditioned media for 1, 2, 3, 4, 5, or 6 days (Figure S2).
Comparison of transcriptional responses to direct and indirect stimulation with TGFβ1
To compare the effects of direct TGFβ1 treatment and exposure to CM, we assessed changes in gene expression relative to that in untreated cells (day 0, Ctr) (Supplementary data 1). Initially, the number of affected genes (change at least 2-fold, padj < 0.05) was significantly greater (several-fold) in MCF10A cells than in MCF7 cells, but after six days of treatment, the effect was stronger in MCF7 cells (Figs. 1A, B and S3A). Only part of this gene pool was common in both cell lines (Figs. 1B and S3A). The response to CM (number of affected genes) was delayed compared to that to direct TGFβ1 treatment in both cell lines. However, it reached the same level by day three in MCF10A cells, while it was stabilized at a low level in MCF7 cells (Fig. 1A). The poor transcriptional response to conditioned media in MCF7 cells putatively resulted from cell death induced after TGFβ1 treatment, which reduced the number of mediator-secreting cells (as suggested by PARP1 cleavage; see Figures S1B, C, and 5B). Nevertheless, the majority of genes affected in cells stimulated by CM exhibited similar changes in expression upon direct TGFβ1 treatment, especially in MCF7 cells (Figs. 1B and S3B).
Hierarchical clustering analysis based on expression trends identified six clusters of genes affected by TGFβ1 in MCF10A and MCF7 cells (Fig. 1C) (Supplementary data 2). Genes activated in both cell lines were grouped into clusters 4 and 5. These gene sets were putatively regulated by, among others, the SMAD3 and SMAD4 transcription factors (prediction based on ChEA3 transcription factor enrichment analysis [30]) (Figure S4A), which were also identified as major upstream regulators in the analysis of all genes activated after treatment (Figure S4B). A similar pattern was visible in cluster 6 but only in MCF10A cells. Assignment to terms from the hallmark collection showed that these clusters were enriched in genes primarily associated with the epithelial mesenchymal transition (M5930) (Fig. 1C). On the other hand, cluster 1 contained genes repressed after treatment in both cell lines (and predicted to be regulated by MYC and E2F family of transcription factors; Figure S4A). This cluster was enriched in genes related to e2f targets (M5925) and the g2m checkpoint (M5901), which both reflect cell proliferation ability, as well as myc targets v1 (M5926) and myc targets v2 (M5928). Clusters 2 and 3 showed the opposite pattern in both cell lines, i.e. repression of expression in one cell line was accompanied by activation in the other. Cluster 3 mainly showed inhibition of the expression of genes related to the estrogen response early (M5906) in the MCF7 line (Fig. 1C). Moreover, estrogen receptor alpha (ESR1/ERα) was found to be an upstream regulator of this set of genes (Figure S4A), as well as the entire set of genes downregulated after TGFβ1 treatment in MCF7 cells (Figure S4B).
More detailed gene set enrichment analyses (GSEAs) based on the whole transcript dataset revealed enrichment in the tgf beta signaling (M5896) and, partially overlapping, epithelial mesenchymal transition (M5930) pathways from the hallmark collection in both cell lines after both TGFβ1 (TGF) and conditioned media (CM) treatment (mainly upregulation) (Fig. 1D; Supplementary data 3). The lack of significant differences between TGF- or CM-treated MCF10A cells (TGF vs CM) indicates a similar response to both treatments. Minor differences between both treatments in MCF7 cells likely result from weaker responses in CM-treated cells. Comparisons of responses to TGFβ1 in MCF10A and MCF7 cells (MCF10A vs MCF7) revealed differences even in signaling pathways enriched in both cell lines, for example, the tgf beta signaling (M5896), epithelial mesenchymal transition (M5930), p53 pathway (M5939), and estrogen response early (M5906) (Fig. 1D). However, GSEA revealed no significant differences in terms reflecting cell proliferation ability, i.e., e2f targets (M5925) and, partially overlapping, the g2m checkpoint (M5901) (with apparent inhibition of gene expression). Nonetheless, genes regulated by MYC (especially a subgroup associated with myc targets v2) were inhibited more effectively after treatment in MCF10A cells (Fig. 1D). Moreover, greater enrichment in TGFβ1-treated MCF7 cells was observed in cholesterol homeostasis (M5892) and angiogenesis (5944), while protein secretion (M5910) was more enriched in MCF10A cells. The comparison of TGFβ1-induced changes after 6 days of treatment shown in the scatterplots (Fig. 2A) further illustrates the differences and similarities between cell lines in the response of individual genes associated with the identified signaling pathways. Notably, however, even in pathways similarly enriched in both cell lines, individual genes responded differently (e.g., CDKN1A, and HMGB3 from the e2f targets pathway or SMAD3, SLC7A5, and DMD from the g2m checkpoint). On the other hand, even in pathways differentially enriched in both cell lines, individual genes responded similarly (e.g., MCM2, MCM4, MCM5, MCM6 from myc targets pathways, encoding components of the minichromosome maintenance protein complex) (Fig. 2A).
Interestingly, we observed large differences between both cell lines at baseline (i.e., without TGFβ1 treatment), mainly in the expression of genes associated with epithelial mesenchymal transition, p53 pathway, and estrogen response early (Fig. 1D). In particular, many genes associated with the first two pathways were expressed at higher levels in untreated MCF10A cells (e.g., those encoding collagens COL16A1, COL4A1, COL4A2, or the matrix metalloproteinase MMP14), while those associated with the estrogen response early were expressed at higher levels in untreated MCF7 cells (Fig. 2B), consistent with their estrogen receptor status. Consequently, SMAD3 and ESR1 were identified as upstream regulators of differentially expressed genes in MCF10A and MCF7 cells, respectively (Figure S4B). Additionally, genes associated with the tgf beta signaling pathway were differentially expressed in both cell types (Fig. 1D), as shown in detail on the map from the KEGG database (Figure S5A).
Differences in the TGFβ signaling pathway between MCF10A and MCF7 cells
Detailed analysis of the expression of genes involved in EMT and TGFβ signaling revealed that many genes (e.g., TGFB1, TGFBR2, SMAD2, NOG, SNAI2, ZEB2, CDH2, FN1, and VIM) were initially (Ctr, day 0) expressed at higher levels in MCF10A cells than in MCF7 cells (Fig. 3). Based on this finding, we assumed that TGFβ1 signaling through TGFBR2 (TGFβ receptor type II), SMAD2 (receptor-regulated SMAD), SNAI2, and ZEB2 may be more important for maintaining the fibrocystic phenotype (manifested by higher expression of CDH2, FN1, and VIM) of the MCF10A cell line than the cancerous phenotype of MCF7 cells. Genes involved in EMT and TGFβ signaling also responded differently to TGFβ1 stimulation in both cell lines (Figs. 2A, S5B). TGFB1 and SMAD3 transcript levels significantly increased (reaching a maximum on days 2–3 of treatment) only in MCF7 cells while they did not change substantially in MCF10A cells (Fig. 3A). Therefore, we postulate that in MCF7 cells, TGFβ1 signaling is more likely to be prolonged through SMAD3 than SMAD2. Consistent with this finding, phosphorylated/active SMAD3 was more efficiently stabilized in MCF7 cells (Figures S1B, C, and 5B). In parallel with the activation of R-SMADs, factors that deactivate the pathway can be upregulated after TGFβ1 treatment. These are SMAD6 and SMAD7, i.e., inhibitory SMADs, or NOG, a potent inhibitor of BMP (bone morphogenetic protein) cytokines, which are known to trigger signaling through SMADs (similar to TGFβ but via SMAD1, SMAD5, or SMAD8). Indeed, we observed SMAD6, SMAD7, and NOG upregulation and only SMAD6 transcript levels in MCF7 cells (initially high) did not change. However, it is noteworthy that SMAD6 and SMAD7 (but not NOG) transcript levels were considerably lower in untreated MCF10A cells (Fig. 3A). This presumably facilitates the maintenance of higher signaling activity through R-SMADs in these cells, resulting in the enrichment of the epithelial mesenchymal transition pathway (110 of 200 genes were expressed at higher levels in MCF10A cells than in MCF7 cells; Fig. 2B). The TGFβ pathway can also be inactivated by SMAD-specific E3 ubiquitin protein ligases (SMURF1 and SMURF2, whose transcript levels increased more strongly after TGFβ1 treatment in MCF10A cells) and latent-transforming growth factor beta-binding proteins (LTBP1, LTBP2, LTBP3), furin (FURIN), and decorin (DCN), whose transcript levels were differentially regulated after TGFβ1 treatment in both cell lines (Figure S5B, Supplementary data 1).
Similar changes in the expression of genes involved in TGFβ signaling were observed in cells directly stimulated with TGFβ1 and those treated with CM (although the response was frequently delayed or weaker in the second scenario), indicating the secretion of TGFβ pathway mediators, including TGFβ1 itself, into the media. Genes encoding other potential regulators of the TGFβ signaling pathway, including TGFB2, BMP2, BMP4, BMP6, BMP7 (BMP cytokines), INHBA, INHBB, INHBE (activins and inhibins), or GDF11, GDF15 (growth differentiation factors) and the receptors BMPR1B, BMPR2, ACVR1, and ACVR1C were upregulated in at least one cell line (Figure S5B and Supplementary data 1).
Among the transcription factors involved in the later stages of TGFβ and EMT signaling, we observed a several-fold increase in SNAI1 (as well as ZEB1 and ZEB2) transcript levels starting from the first day of treatment only in MCF10A cells (stimulated by either TGFβ1 or CM), while SNAI2 was more effectively upregulated in MCF7 cells (but its transcript was still at a lower level than that in MCF10A cells) (Fig. 3B). Surprisingly, increased expression of SNAI1 and SNAI2 was accompanied by an increase (albeit small) rather than the expected decrease in CDH1 transcript levels in MCF10A cells, while the decrease in CDH1 levels was not significant in MCF7 cells (Fig. 3C). However, the expression of mesenchymal markers, CDH2, VIM, and FN1, was induced to some extent in both cell lines. Nevertheless, western blot analysis showed that the CDH1 protein levels decreased in both cell lines after TGFβ1 treatment and confirmed the accumulation of VIM in MCF10A cells, reflecting its expression at the mRNA level (Fig. 5C). In conclusion, although the TGFβ and EMT signaling pathways were enriched in both cell lines starting from the first day of treatment with TGFβ1 or conditioned media, there were considerable differences between the cell lines (partially resulting from differences existing before treatments).
TGFβ1 treatment affects the expression of genes associated with cell cycle progression
In addition to inducing EMT, TGFβ is known to inhibit cell cycle progression by blocking the G1 phase resulting from decreased expression of the MYC proto-oncogene [35] and IDs (inhibitors of DNA binding) [36]. Accordingly, MYC and ID1 expression was inhibited, but only in MCF10A cells (Fig. 4A). Consequently, MYC targets (M5926 and M5928 from the hallmark collection) were also largely downregulated mainly in MCF10A cells (Figs. 1D and 2A). In contrast, the expression of ID3 was upregulated (more clearly in MCF7 cells; Figure S7A). Nevertheless, the downregulation of genes involved in cell cycle progression, both during the G1/S transition (e2f targets, M5925) and during the G2/M transition (g2m checkpoint, M5901), was also observed in MCF7 cells (Fig. 1D; see also Figure S6A). However, our data revealed that cell cycle inhibition could be achieved by different mechanisms, such as MYC inhibition in MCF10A cells and a p53-dependent mechanism in MCF7 cells (Figs. 1D, 2A, and S6B). The strong upregulation of TP63 (according to ChEA3, it may be responsible for the upregulation of ~ 9% of genes in MCF7 cells; Figure S4B) and cyclin-dependent kinase (CDK) inhibitors: CDKN1A (encoding p21) and CDKN2B (encoding p15) was observed specifically in MCF7 cells (Fig. 4A, D; although some other CDKN genes, e.g., CDKN2C, were downregulated; Figure S7B).
Therefore, growth arrest and DNA damage (GADD45A and GADD45B) genes were upregulated (Fig. 5D), while genes encoding activators of cell cycle progression were downregulated (Figure S6A) after TGFβ1 treatment in both cell lines. This phenomenon applies to several family members of the E2F transcription factors, namely, E2F1, E2F2, E2F4, and E2F8 (Figs. 4A, S7A), their dimerization partner TFDP1, and the pocket protein RBL1, but not RB1 and RBL2 (Figure S7A), cyclins and CDKs, i.e., G1-specific CCNE2, CCNA2 and its targets: CDK2 during S phase, and CDK1 during the transition from G2 to M, G2/M-specific: CCNB1, CCNB2, CCNF, but not cyclins D (Figs. 4B, C, S7B), cell division cycle proteins: CDC6, CDC7, CDC20, CDC25A, CDC25C, CDC45 (Figs. 4E, S7C), and cell division cycle associated proteins, such as CDCA2, CDCA3, CDCA5, CDCA7, CDCA7L, and CDCA8 (Figs. 4F, S7D). Furthermore, some genes encoding DNA polymerases, e.g., POLA1, POLA2, POLD1, POLD2, POLE, and POLE2 (Figs. 4G, S7E), and proteins involved in chromosomal replication, e.g., origin recognition complex subunits, ORC1 and ORC6, as well as all minichromosome maintenance complex components, MCM2-7, centromere proteins, CENPA, CENPE, CENPF, CENPH, CENPM, and others (Figs. 4H, S7F), were downregulated in both cell lines. Interestingly, some genes encoding proteins involved in DNA repair (BRCA1, BRCA2, DCLRE1A, DCLRE1B, ERCC6L, and MRE11) were more effectively inhibited in MCF7 cells (Figs. 4I, S7G). Consequently, the marker of proliferation MKI67 and the proliferating cell nuclear antigen PCNA were downregulated in both cell lines (Fig. 4J).
Inhibition of estrogen signaling and induction of cell death after TGFβ1 treatment in MCF7 cells
Treatment with TGFβ1 affected genes related to estrogen signaling, which is essential for maintaining the phenotype of MCF7 cells. Specifically, many members of the estrogen response early (hallmark M5906) were downregulated in the MCF7 cell line (Fig. 1D). This may be due to the downregulation of the estrogen receptor ESR1, both at the mRNA and protein levels (Fig. 5A, B), which was correlated with decreased expression of several direct ESR1 targets (known to be stimulated by estrogen; [37]), e.g., RET, TFF1, and IGFBP4 (Fig. 5C). It is noteworthy that transcription factor enrichment analysis (using ChEA3) identified ESR1 as a regulator of ~ 10% of all TGFβ1-downregulated genes (98 out of 975) in MCF7 cells (while downregulated E2F family members were identified as transcriptional regulators in both cell lines; Figure S4B). This indicates that reducing ESR1 expression may have a large impact on the expression profile and may contribute to the inhibition of proliferation in TGFβ1-treated MCF7 cells. TGFβ1 treatment also inhibited the expression of PGR (progesterone receptor), while GPER1 (G-protein coupled estrogen receptor) was upregulated (Fig. 5A). TGFβ1 can also deregulate various signaling pathways (NOTCH, WNT/β-catenin, PI3K/AKT, EGFR, RAS/RAF/MAPK) known to stimulate cell growth, survival and differentiation in other molecular subtypes of breast cancer (Figure S8; specifically, the activation of epidermal growth factor receptor, EGFR, expression in MCF7 cells is of interest; Fig. 5A).
The induction of cell death in MCF7 cells was the most important functional difference between MCF7 and MCF10A cells in response to TGFβ1. This could be executed by caspases 7 and 9 (caspase 3 is not expressed in MCF7 cells due to a mutation leading to a frameshift in the CASP3 gene). The accumulation of active forms caspases correlated with PARP1 cleavage (Fig. 5B) and indicated the induction of caspase-dependent apoptosis. Analysis of the apoptosis network (KEGG ID: hsa04210; Figure S9) revealed upregulation of some prosurvival genes in one cell line or the other (e.g. BCL2A1/A1, GADD45A, and GADD45B; Fig. 5D). On the other hand, the expression of proapoptotic genes: BAK1 (BAK), BAX (BAX), BBC3 (PUMA), BCL2L11 (BIM), FAS, HRK, TP53AIP1 (p53AIP1), TNFRSF10A and TNFRSF10B (TRAILR2) was preferentially induced in MCF7 cells (only BMF was similarly induced in both cell lines; Figs. 5E and S9). We assumed that apoptotic MCF7 cells had reduced secretory functions (hallmark protein secretion collection was enriched only in MCF10A cells; Fig. 1D), which resulted in a weak propagation of TGFβ1 signaling through conditioned media (particularly in terms related to cell cycle progression; Fig. 4).
The observed proapoptotic and antiestrogenic effects of TGFβ1 treatment in MCF7 cells may suggest the prognostic value of TGFB1 expression in breast cancer patients. Interestingly, the Kaplan Meier analysis of a large cohort of breast cancer patients (n = 2976) integrated by Győrffy [38] (available at https://kmplot.com/analysis/index.php?p=service&cancer=breast_rnaseq_gse96058) revealed that high TGFB1 levels are associated with a better prognosis, but only in luminal A patients (Figure S10). These tumors are characterized by the presence of ER and/or PR, the absence of HER2, and low expression of Ki-67.