Estrogen receptor β target gene expression reveals novel repressive functions in aggressive breast cancer.

Introduction
Inflammatory breast cancer (IBC) is an aggressive form of breast cancer with a unique clinical presentation and poorly defined biology that accounts for a substantial proportion of breast cancer-related deaths in the United States each year1,2. As with other forms of aggressive breast cancer, the combined use of neoadjuvant chemotherapy with loco-regional surgery and radiotherapy has improved clinical outcomes. However, the 5-year survival rate remains around 30%, significantly lower than that of non-IBC2–4. The poor outcome has been linked to the intricate biology of IBC and the lack of effective targeted therapies4. With a low incidence of estrogen receptor αlpha (ERα)-positivity in IBC, most patients do not benefit from endocrine therapy, underscoring the need to develop new molecular targets2.
Efforts to identify novel molecular determinants of IBC aggressiveness through genomic analysis of human tumors have been unsuccessful5. Similarly, functional studies focusing on pathways involved in actin-based cell migration, lymphatic invasion and tumor microenvironment responses have not substantially advanced our understanding of the etiology of the aggressive IBC tumors5–7. Ongoing clinical trials are primarily exploring treatments used for other types of cancer, highlighting the need for deeper understanding of IBC biology to discover specific molecular targets for this form of breast cancer8,9.
Given the lack of ERα in approximately 70% of IBC tumors, we investigated whether the second ER subtype, ERβ, mediates estrogen effects in IBC tumors10. Consistent with the proposed tumor suppressor role of ERβ, we and others have initially reported its anti-metastatic effects in triple-negative breast cancer (TNBC) through the inhibition of epithelial-to-mesenchymal transition (EMT)11–14. More recently, we demonstrated similar anti-tumor activity of ERβ in IBC, associating its tumor expression with reduced metastasis in patients and xenografts, as well as with inhibition of actin-based cell migration and Rho GTPase signaling15,16. However, our previous targeted approach did not fully elucidate the underlying signaling involved in the anti-metastatic activity of ERβ, emphasizing the need for a comprehensive genomics approach to study our IBC models.
We applied combined chromatin immunoprecipitation followed by sequencing (ChIP-seq) and RNA-seq in various IBC cell lines with different ERβ levels, both treated and untreated with ERβ ligands. By comparing the ERβ transcriptome and cistrome in cancer cells that express endogenous and transfected ERβ, we demonstrated gene expression profiles supporting the anti-tumor function of endogenous ERβ. Our approach further identifies ERβ as an essential factor in tissue differentiation during development and a regulator of metabolism in cancer progression and metastasis. Finally, we identified novel target genes that may represent key mediators of estrogen signaling in aggressive cancers and potentially unique targets for therapeutic intervention.

Results

ERβ-DNA binding in inflammatory breast cancer cells
We previously showed that ERβ suppresses metastasis in IBC by inhibiting actin-based cell migration and its associated Rho GTPase signaling15. Our current study aims to elucidate the mechanisms of gene regulation employed by ERβ to control metastatic signals in aggressive forms of breast cancer, using IBC as a model due to its metastatic behavior and availability of suitable cell lines with endogenous ERβ (Fig. 1A)2,15.
To identify enriched ERβ binding sites, we upregulated the receptor in two IBC cell lines and activated it with ligand treatment, followed by chromatin immunoprecipitation and sequencing (ChIP-seq) (Fig. 1A). We conducted these experiments in ERβ knockout (ERβΚΟ) HER2-positive IBC KPL4 cells and in ERβ knockout cells re-expressing FLAG-tagged ERβ (ERβKO + ERβ), using an anti-FLAG antibody for immunoprecipitation as previously described (Fig. 1A)15. Additionally, we examined ERβ-DNA binding in control and ERβ-expressing SUM149 TNBC IBC cells (Fig. 1A). We also analyzed ERβ-bound chromatin in ERβKO + ERβ cells following hormone deprivation for three days and subsequent activation of the receptor by incubation with the selective agonist LY500307 for 2 h, optimizing ERβ-chromatin binding13,17. We carried out three biological replicates, and the sequencing generated between 25 M and 65 M reads per sample, except for replicate 3 in ERβ-expressing KPL4 cells, which only generated ∼5 M reads, below the minimum cutoff of 10 M reads needed for accurate transcription factor binding identification. Replicates with different ERβ levels and upon ligand treatment showed high peak reproducibility, except for replicate 3 in KPL4 and SUM149 cells with altered expression of the receptor, which were excluded from further analysis due to low reproducibility and insufficient read depth. Heat maps in Fig. 1B and C illustrate the enrichment of ERβ-bound DNA after comparing DNA binding in ERβKO and ERβKO + ERβ cells and treatment with vehicle and LY500307. High peak reproducibility was observed with 1988, 1544, and 1885 common sites among replicates in each cell line and ligand treatment (Fig. 1D, E, and F). After combining evidence for ERβ-bound DNA across replicates using the multiple sample peak calling (MSPC) approach, 3399, 45243, and 11583 significant sites were considered for ERβ-expressing KPL4 and SUM149 cells and ERβKO + ERβ cells after treatment with LY500307 agonist, respectively.

Distribution of ERβ-DNA binding and transcription factor motifs across the genome
To understand the transcriptional mechanism of ERβ action in IBC, we mapped the relative location of significant ERβ-binding sites in the genome of ERβKO + ERβ KPL4 and ERβ-expressing SUM149 cells, including after treatment with the ERβ agonist LY500307. We considered both binding sites proximal to the promoter (mapped to −1000 bp to +100 bp from the transcription start site, TSS) and cis-regulatory elements. As shown in Fig. 1G and I, both upon upregulation of ERβ and activation by agonist treatment in KPL4 cells, only a small percentage of ERβ-binding sites were within the promoter areas (6.2% and 17.5%, respectively) and exons (9.9% and 19.1%), while most were within introns (51.3% and 39.3%) or intergenic regions (32.7% and 24.1%), consistent with ERα- and ERβ-binding patterns in other breast cancer cell lines13,17,18. In contrast, SUM149 cells showed a different enrichment pattern with a similar number of sites within introns, exons and near the TSS (Fig. 1H). To identify partner transcription factors of ERβ, we scanned the ERβ-binding sites for individual motif occurrences (p < 0.0001). This analysis revealed the canonical estrogen response element (ERE) as the most abundant motif in KPL4 cells following upregulation of ERβ and activation by LY500307 (Fig. 1J, L). In SUM149 cells where we detected more significant peaks than KPL4 cells, ERE motifs appeared in ~50% of peaks (Fig. 1K and Supplementary Fig. 1). Additionally, KLF and TCF motifs were prominently enriched in both KPL4 and SUM149 cells, both upon upregulation of ERβ and activation by agonist treatment, suggesting a universal tethering mechanism in IBC involving networks with specific ERβ-interacting factors (Fig. 1J, K and L). We also observed an overrepresentation of the well-known interacting motifs such as Forkhead box (FOX) and AP-1, as well as zinc finger response elements (ZNF460, ZNF263), which were not previously reported as enriched motifs in ERβ-bound chromatin (Fig. 1J–L)17–19.
To confirm the association of ERβ with specific transcription factors, we identified enriched motif patterns using the Spamo algorithm. We obtained significant motif hits for ERβ (ERE, primary) and its partner transcription factors (secondary) within the provided binding sites and identified significant co-occurrences at specific distances from each other. Notably, SUM149 cells exhibited more significant secondary motif hits compared to KPL4 cells, attributed to a higher number of detected binding sites in SUM149 cells (Fig. 1E). In KPL4 cells, following ERβ upregulation or LY500307-mediated activation, motifs like KLF13, ZNF460, RPN4, MET31, RAV2, and DOF5.8 were found at a significant distance from ERE (Supplementary Fig. 2). Similarly, in SUM149 cells, motifs including ZNF460, FOXM, AP1, GAF1, DOF5.8, MEF2A, STAT2, and FOXO, showed similar enrichment patterns. Interestingly, both ligand-treated KPL4 and ERβ-expressing SUM149 cells shared common enriched motifs (RPN4, MET31, RAV2, ZNF460, and DOF5.8) positioned identically relative to ERE, indicating the formation of specific transcription factor complexes during the ERβ-DNA binding.

ERβ-DNA binding in inflammatory breast cancer
To further investigate ERβ-DNA binding in IBC, we examined common binding sites in KPL4 and SUM149 cells using ChIP-seq data. We also compared the DNA binding of ERβ in KPL4 cells following upregulation of the receptor in serum-containing media and activation by LY500307 in estrogen-depleted media. We identified 807 common sites between the two cell lines and 2348 sites across different conditions of ERβ upregulation and activation (Fig. 2A, B). In addition, we scanned the common binding sites for enriched DNA-binding elements. The search confirmed the overrepresentation of EREs and captured the previously reported motifs like KLF, AP1, FOX, and TCF among the top significant and frequently occurring motifs in ERβ binding sites within IBC cells (Fig. 2C, D)17,18. Pathway enrichment analysis of the genes located near these 807 and 2348 common binding sites revealed significant (p < 0.0001) biological processes associated with cell development, differentiation, and organ morphogenesis, critical in breast cancer development and progression (Fig. 2E, F). ERβ-associated functions observed in both IBC cell lines were similarly enriched following both receptor upregulation and ligand activation, emphasizing the specificity of ERβ binding in IBC cells (Fig. 2E, F).
Additionally, we analyzed the co-occurrence of the top four most frequent motifs in common binding sites located within 20,000 bp upstream and downstream of the TSS in both KPL4 and SUM149 cells. As shown in the Venn diagram in Fig. 2G, the majority of binding sites contained combinations such as ERE and KLF motifs (335), ERE and TCF (300), ERE and AP1 (258) and all four motifs together (214). In contrast, a very small subset of common binding sites contained exclusively ERE (2), KLF (5), TCF (1), or AP1 (1) sequences, suggesting that ERβ likely regulates gene expression in IBC cells primarily through interactions with other transcription factors, often involving direct binding of the receptor to the cognate ERE (Fig. 2G).
To examine whether specific functions of ERβ are mediated through direct DNA binding or interactions with other DNA motifs, we searched for enriched functions associated with genes located near binding sites containing ERE motifs and other than ERE motifs17,18. As shown in Fig. 2G, the number of common binding sites containing EREs (354) was significantly higher than those without (49). Our analysis links ERβ direct binding to functions related to development and differentiation, whereas its interaction with other transcription factors is associated with functions like cell migration, cytosolic, and endosomal transport (Fig. 2H).

ERβ-associated gene expression correlates with DNA-binding
In addition to identifying genes located near ERβ-binding sites, we used RNA-seq tο examine changes in gene expression in IBC cells with varying levels of ERβ. Based on our previous findings of endogenous ERβ expression in IBC cells15, we compared, for the first time, gene expression in cells with both endogenous and transfected ERβ, which has been associated with tumor repressive activity11,12,14,15,19–22. We also assessed mRNA changes in ERβ-expressing cells following treatment with the agonist LY500307.
As shown in the multi-dimensional scaling analysis of RNA-seq data (Supplementary Fig. 3A), replicates cluster together, while different conditions were distinct and clearly separated. The analysis identified significantly differentially expressed genes that exhibited an absolute log2FC ≥ 0.585 (≥1.5-fold difference and p < 0.05) due to changes in ERβ expression and activity, listing the top 1000 altered genes in Supplementary Data 3-5. The highest number of differentially expressed genes (3335) was detected in KPL4 cells after knocking out ERβ, followed by upregulation of the receptor in knockout cells (1007) and its activation in ERβ-expressing cells with LY500307 (125) (Fig. 3A).
Supplementary Fig. 3B shows a significant number of upregulated (171) and downregulated (397) genes overlapping in KPL4 cells with endogenous and transfected ERβ. As expected, heat map analysis of the top 20 most highly up- and down-regulated genes in KPL4 cells with different ERβ levels and activity identified similar gene expression profiles, with a significant (20%) proportion of common top altered genes (8/40) in two different heatmaps (Fig. 3B, C). These primarily include novel ERβ target genes that regulate differentiation (CRISP3, KLK5), tumor immunity (SPP1, ANXA1, VTCN1), metabolism (KLK5, CES1), as well as previously reported targets (ELMO1) that regulate actin-based cell migration and mediate the anti-metastatic activity of ERβ in IBC (Fig. 3B, C)15. These findings represent the first unequivocal proof of the presence and functionality of endogenous ERβ, which, similar to the transfected receptor, regulates the expression of genes that are linked to its tumor suppressor activity11,12,14,15,19–22. Common altered genes (CA12, SYTL5) were also found among the top 20 up- and down-regulated genes in IBC cells with altered ERβ expression and those treated with the ERβ agonist, further validating the specificity of our sequencing analysis and the quality of ERβ-expressing models (Fig. 3B–D).
To differentiate direct targets of ERβ from those indirectly regulated by the receptor, we correlated the ERβ binding sites from our MSPC analysis in KPL4 and SUM149 cells with the significantly altered genes (FDR < 0.05) in KPL4 cells with different ERβ levels and after their treatment with LY500307 and validated the integration of ChIP-seq peaks with differentially expressed genes using the Binding and Expression Target Analysis (BETA) method (Supplementary Fig. 4A). The heatmap in Supplementary Fig. 4B depicting significantly altered genes in KPL4 cells with different ERβ levels that contain ERβ binding sites in promoter, gene body, upstream of the promoter and downstream of the gene, shows similar trends of target gene regulation in cells with endogenous and transfected receptor. As shown in Fig. 3E, F, several common genes (TFF1, TRIM31) were among the top 10 up- and down-regulated genes with ERβ binding sites within and upstream of the promoter in cells with varying levels of ERβ and those following treatment with the ERβ agonist. Similarly, differentially expressed genes (BACE2) with ERβ-binding sites in KPL4 and SUM149 cells suggest a universal mechanism of ERβ transcriptional regulation in IBC (Fig. 3G).
Finally, we verified both the CHIP and RNA sequencing results by ChIP-qPCR and qPCR, respectively, using a number of selected significant binding sites and genes (Fig. 2H, I). qPCR confirmed not only the genes detected in RNA-seq analysis but also the mode of alteration in the presence of endogenous and transfected ERβ (Fig. 2I).

The transcriptional activity of ERβ dictates regulation of metastasis-associated biological functions
We previously demonstrated that ERβ elicits anti-metastatic effects in IBC15. To fully define the signaling pathways and biological processes that are linked to the transcriptional activity of the receptor and mediate its anti-metastatic effects, we carried out an enrichment analysis against GO terms (biological functions) and KEGG pathways (molecular pathways) for genes identified as either differentially expressed or located near ERβ-bound chromatin.
As shown in Fig. 4A, testing differentially expressed genes in cells with endogenous ERβ revealed an over-representation of functions associated with tissue morphogenesis, organ development, differentiation and hormone-associated metabolic processes. These same processes were also enriched in the analysis of differentially expressed genes in cells with transfected ERβ, further supporting the functionality of the endogenous receptor in IBC cells (Fig. 4B). Interestingly, most of these enriched functions were associated with ERβ down-regulated genes, linking the repressive effect of the receptor on target gene expression to the control of differentiation (Supplementary Fig. 5A).
When the analysis included genes that were only located near ERβ-binding sites, more specific molecular functions related to cell adhesion, polarity, actin re-organization, and the RhoGTPase pathway were detected (Fig. 4C and Supplementary Fig. 5B). These cellular functions, known to control tissue morphogenesis, may relate to the processes of development and differentiation enriched during the analysis of differentially expressed genes (Fig. 4A, B). They also strengthen our previous findings correlating ERβ expression in IBC cells with the inhibition of actin-based cell migration and its associated RhoC GTPase signaling15.
To verify the relevance of the identified functions, we employed pathway analysis considering the direct targets of the receptor- common differentially expressed genes in cells with endogenous and transfected ERβ that also harbor ERβ-binding sites. We again noted the enrichment of biological terms, indicating the involvement of the receptor in development and differentiation (Fig. 4D and Supplementary Fig. 5C). Additionally, there was a noteworthy over-representation of terms associated with metabolic processes, including fatty acid metabolism, which are essential components and drivers of aggressive breast cancers. These metabolic themes had not previously been detected as ERβ-regulated during transcriptomic analysis of ERβ in breast cancer (Fig. 4D and Supplementary Fig. 5D). Similar functions were observed in an additional KEGG pathway analysis, reinforcing the essential role of metabolism in mediating the effects of ERβ in IBC (Supplementary Fig. 5E).
To examine whether the enriched signaling pathways and biological processes that are associated with the ERβ transcriptional activity in IBC cells overlap with those in non-IBC breast cancer, we considered the genes that are located near overlapping and non-overlapping ERβ binding sites in LY500307-treated ERβ-expressing IBC (KPL4) and previously published ERβ-expressing TNBC (MDA-MB-231) cells13. We also examined non-overlapping differentially expressed genes in the same cells. We observed biological processes that are enriched in both KPL4 and MDA-MB-231 cells (epithelial cell differentiation, organ development), but also GO terms that are uniquely enriched in KPL4 (stem cell development, protein transport) and MDA-MB-231 (TGFβ signaling, cell division, DNA damage response) cells (Fig. 4E, F and Supplementary Fig. 5F–H). This is consistent with recent reports that identified differences in enriched hallmark gene sets between IBC and TNBC, including those regulating protein secretion in IBC and regulation of cell cycle and DNA damage response (CCND1, CHEK2, TP53, FEN1) in TNBC23. These results also suggest that the tumor suppressor function of ERβ in different forms of breast cancer can be mediated by unique in each type transcriptional activities.

ERβ regulates the expression of miRNAs in inflammatory breast cancer
In addition to protein-coding mRNAs, we assessed small RNAs in IBC cells with different levels of ERβ expression by sequencing the small RNA fraction of ERβKO and ERβΚΟ + ΕRβ KPL4 cells using miRNA libraries. As shown in Supplementary Fig. 6A, multi-dimensional scaling analysis of small RNA-seq data clearly separated the replicates of the two different conditions. According to the volcano plot of Fig. 5A and the heat map in Supplementary Fig. 6B, a significant number of miRNAs were altered in ERβ knockout cells. Among the top five downregulated miRNAs in these cells are the miR-548-3p, miR-873-5p, miR-381-3p, miR-100-5p, and miR-216a-5p, which have been associated with tumor suppressor functions in breast cancer, consistent with the anti-tumor activity of ERβ in the same IBC cells24–28.
To strengthen the involvement of ERβ in the regulation of these miRNAs, we correlated the miRNA expression data with the ERβ binding sites containing the ERE motif identified in our ChIP-seq analysis of SUM149 cells. These binding sites were within and close to the promoters of miRNAs, from −4000 bp upstream to +100 bp downstream of the gene, using TSSs derived from complementing the available information from CAGE analysis of MDA-MB-453 cells that is used to identify TSSs with the TSSs of host genes for intragenic miRNAs (Fig. 1)29,30. We identified 52 significantly differentially expressed miRNAs located near ERβ binding sites, suggesting that these miRNAs are direct transcriptional targets of ERβ (Fig. 5B). We confirmed their differential expression and pattern of alteration by TaqMan PCR assay (Fig. 5C).
To define the network of transcription factors interacting with ERβ to regulate miRNAs, we scanned the cognate ERβ binding sites for enriched DNA binding motifs. Our screen revealed a similar list of DNA-binding elements (KLF, TCF, ZNFs) to those enriched in sites of nearby protein-coding genes, with the exception of ERE, which was not among the top-enriched motifs. This suggests that tethering with other transcription factors is the primary mechanism of interaction of ERβ with chromatin in the regulatory regions of miRNAs (Fig. 5D).
To delineate the functions of the receptor associated with the transcriptional regulation of miRNAs, we searched our RNA-seq data for differentially expressed mRNAs in ERβ-expressing KPL4 cells that are known targets of the ERβ-regulated miRNAs, using the TarBase v9 and microT-CDS 2023 databases. These databases provide both experimentally verified and predicted miRNA-gene interactions. As shown in Supplementary Fig. 7, we identified many miRNA-targeted mRNAs that follow an inverse expression pattern with that of miRNAs upon upregulation of ERβ using Spearman’s correlation (ρ < −0.8). We then divided the significantly deregulated genes in ERβ-expressing cells based on the presence of nearby ERβ binding sites and their interaction with miRNAs, distinguishing those regulated by direct ERβ binding from those that were regulated by miRNAs (Fig. 5E). Among these genes are previously identified ERβ targets (GPR141) as well as factors implicated in tumor immunity (CCL22) and breast cancer metastasis (Fig. 5E)15,31.
To better understand the miRNA-mediated effects of ERβ, we carried out KEGG pathway analysis of the target genes of the top five altered miRNAs. We noted upon knockout of ERβ, a significant enrichment of signaling pathways involved in breast cancer metastasis, including those associated with metabolism and actin-based cell migration, similar to the functions enriched during the analysis of target genes (Fig. 5F, G and Fig. 4). We also observed functions not in the list of top altered GO terms in the pathway analysis of direct target genes of ERβ (Fig. 4), including the circadian rhythm, endocytosis and cell cycle. These may represent unique functions that are mediated through miRNAs, underscoring the indirect gene regulation as an important aspect of the mechanism of ERβ action.

ERβ target genes are associated with clinical outcome in IBC
Given the significant association of ERβ with reduced metastasis in IBC patients15, we hypothesized that target genes downstream of the receptor mediating its anti-metastatic effects would similarly correlate with clinical outcomes. For this, we subjected differentially expressed genes in IBC cells with varying ERβ levels or after treatment with the ERβ agonist LY500307, located near ERβ binding sites, to rigorous analysis using a dataset with gene expression and clinical information from 25 IBC patients and 58 non-IBC patients32.
We performed Kaplan-Meier analysis by dividing patients into “low” and “high” expression subgroups based on optimal cutoff values. At a 10% False Discovery Rate (FDR), we identified fifteen target genes significantly associated with overall survival in IBC patients (Fig. 6A). Of these, nine were direct targets in KPL4 cells and six in SUM149 cells. These include the ERβ-repressed gene, paired-like homeodomain transcription factor 1 (PITX1), previously described as an ERα transcriptional target and metastasis-associated factor in breast cancer33,34 correlating with poor overall survival in IBC patients (FDR = 5.6 × 10−2) (Fig. 6B and Supplementary Data 3–5). Similarly, another ERβ-downregulated gene, HOMER3, which is a pro-metastatic factor and mediator of growth factor signaling, has been associated with metastasis and clinical outcome in TNBC35. Our findings link higher expression of HOMER3 to worse survival in IBC cohort (FDR = 8.7 × 10−2) (Fig. 6C and Supplementary Data 3–5).
Our findings also link the high expression of SERPINA1 (an ERα-direct target), which is associated with favorable prognosis in ERα+ and ERα + /ΗΕR2+ breast cancer patients to better clinical outcomes in IBC (FDR = 1.8 × 10−2) (Fig. 6D and Supplementary Data 3–5)36. The small heat shock protein B8 (HSPB8) is another estrogen-responsive gene associated with better survival in IBC patients (FDR = 2.7 × 10−2) (Fig. 6E and Supplementary Data 3–5)37. Additionally, protein tyrosine phosphatase μ (PTPRM) and signal peptide-CUB-epidermal growth factor-like domain-containing protein 2 (SCUBE2), recognized as favorable prognostic factors in invasive breast cancer and tumor suppressors, through their reversal of EMT and inhibition of invasion correlated with better survival in patients with IBC (FDR = 2.9 × 10−2) (Fig. 6F and Supplementary Fig. 8A)11,38,39. Among ERβ target genes that correlate with clinical outcome in patients with IBC, only the estrogen-regulated gene TFF1 was associated with better clinical outcomes in both IBC and non-IBC breast cancer patients (Fig. 6A and Supplementary Fig. 8B, C). TFF1 was previously found to inhibit the migration and invasion of breast cancer cells, and its expression was associated with better prognosis in TNBC40. In contrast, several genes that were identified as ERβ targets in IBC cells were associated with clinical outcomes in patients with non-IBC breast cancer (Supplementary Fig. 9A). Among these are the estrogen-regulated gene RET, associated with endocrine resistance and metastasis, and UBTD1, which relates to breast cancer immunity and response to immunotherapy41,42. To strengthen the clinical importance of the identified ERβ targets, we analyzed them using TCGA data from breast cancer patients of either subtype. We identified significant associations between the expression of several ERβ targets and survival of patients (Supplementary Fig. 9B). These genes include those associated with clinical outcomes in IBC (HSPB8) and others that are unique to TNBC, such as LMO7, an activator of RhoGTPase signaling and regulator of breast cancer metastasis, and the estrogen-dependent and pro-metastatic factor CREB5 (Supplementary Fig. 9B)43,44. In addition, we found genes associated with clinical outcomes in unselected breast cancers, including the tumor suppressor TP63 and the well-characterized ERα target GREB1 (Supplementary Fig. 8D, E and Supplementary Fig. 9B). Notably, the ERβ targets associated with clinical outcomes in both IBC and non-IBC encode for factors with a wide range of functions, indicating the diverse nature of pathways that are employed by ERβ to mediate its actions in cancer cells. Το further verify the functional association of the clinically important genes with ERβ, we cloned regions containing the identified by ChIP-seq ERβ binding site of SERPINA1, HSPB8 and TFF1 into a luciferase reporter construct and assessed the transcriptional activity in luciferase reporter assays, using as a positive control the promoter of TRIM31 that was verified to bind ERβ in a ChIP-qPCR experiment (Fig. 3H). In addition, we analyzed the expression of these genes in IBC cells with different ERβ levels. We observed increased expression of luciferase in the presence of the heterologous ERE-containing promoter in ERβ-expressing cells, indicative of ERβ binding and, at the same time, enhanced expression of the target genes, confirming the direct involvement of ERβ in transcriptional upregulation of favorable prognostic factors in IBC (Fig. 6G, H).

ERβ target genes regulate cell migration in IBC
We previously linked the anti-metastatic activity of ERβ with the inhibition of RhoGTPase signaling and actin-based cell migration15. We have now focused on identifying target genes of ERβ that mediate these effects of the receptor. We searched the list of the top ERβ-repressed genes for factors that drive cell migration by regulating the actin cytoskeleton. Among the top sixty ERβ downregulated genes, ten are reported as drivers of actin cytoskeleton remodeling (RELN, PREX1, DIAPH2, GPR141, AGAP11, ADCY2, ELMO1, KIRREL1, GRIK3, SDC2), including genes (PREX1, SDC2, ADCY2, ELMO1) that regulate RhoGTPases (Supplementary Data 3–5)45,46. We examined whether PREX1, a guanine nucleotide exchange factor (GEF) that directly activates RhoGTPases, ADCY2, which is stimulated by G proteins and SDC2, a transmembrane protein that is involved in cytoskeleton organization, induce migration and increase the activity of RhoC that drives actin-based cell migration and has been associated with metastasis in IBC in the absence of ERβ47,48. We confirmed the upregulation of PREX1, ADCY2 and SDC2 in ERβKO IBC cells and their downregulation after siRNA-mediated knockdown (Fig. 7A). Knockdown of PREX1 fully reversed the increased migration that was observed in absence of ERβ (adjusted P < 0.05), ADCY2 downregulation caused a significantly less profound regression of migration and SDC2 did not have any effect (Fig. 7B). Proportional to the effect on cell migration, downregulation of PREX1 fully inverted the increased RhoC activity in ERβΚΟ cells while ADCY2 knockdown displayed less potent inhibitory effect, strengthening the association of ERβ signaling with the inhibition of actin-based cell migration and associated metastasis as well as the importance of its target genes for the IBC biology and as potential therapy targets (Fig. 7C).

Discussion
Considering the previously reported association of non-genomic estrogen signaling with the development of more aggressive cellular phenotypes in IBC, we aimed to fully understand how estrogen impacts this aggressive form of breast cancer49. Due to the higher incidence of ERα negativity in IBC, we focused on ERβ to explain the effects of the hormone2. We initially linked the presence of ERβ in tumors with longer time to metastasis in patients and demonstrated the anti-metastatic activity of the receptor in preclinical models of the disease, reinforcing the previously developed notion of the tumor repressive role of ERβ in breast cancer13–15,20–22. We also associated the anti-metastatic effects of ERβ with the inhibition of actin-based cell migration and its associated RhoGTPase signaling15,16. However, we expected estrogen to mobilize additional pathways and engage a wide spectrum of biological processes to elicit effects in cancer tissues. Considering this hypothesis, we followed an unbiased genomics approach in the present study to corroborate the anti-metastatic function of ERβ and uncover unknown mediators of estrogen signaling that are biologically important and potentially viable therapy targets.
We applied a combined ChIP-seq and RNA-seq analysis in cell models that closely mimic the physiological activation of ERβ signaling, utilizing unique IBC cell lines with demonstrated endogenous expression of the receptor15, unlike previous studies assessing the ERβ DNA binding in breast and other types of cancer cells with artificially introduced receptor17,50–54. We compared gene expression in cells with endogenous and transfected ERβ and the DNA binding in transfected cells prior to and after activation of the receptor with a selective agonist. After cross-examining responsive genes against the chromatin binding profiles, we identified functional targets of the receptor that represent classical estrogen-responsive genes (IGFB4, GREB1), indicating the specificity of our methodology and the relevance of the classical hormone signaling for the anti-metastatic mechanism of ERβ action. Multiple common differentially expressed genes in cells with endogenous and transfected ERβ unequivocally demonstrate the functionality of the physiologically expressed receptor. This finding further addresses the controversy surrounding the expression of ERβ in breast cancer that has dominated the field since its discovery by demonstrating the presence of the receptor in cancer tissues55,56. Moreover, our experimental approach enabled us to define the basis of the ERβ transcriptional mechanism of action in IBC by identifying common DNA-binding sites and target genes across different cell models of the disease. We also complemented our approach for the first time by sequencing small RNAs and defining the ERβ-regulated miRNAs and their mRNA targets. By comparing these mRNAs with the altered genes close to ERβ binding sites, we were able to confidently separate the direct target genes of the receptor from those regulated by miRNAs.
Functional analysis of the identified targets dictates biological roles that account for the previously observed anti-metastatic activity of the receptor in our preclinical models of IBC15. Multiple enriched GO terms confirm the small GTPases and actin filament organization as mediators of active pro-metastatic signaling in IBC, which is inhibited by ERβ. These metastasis-associated processes were not among the top altered in non-IBC TNBC cells with active ERβ signaling13. Instead, ERβ appears to elicit anti-metastatic effects in TNBC by regulating TGFβ and the associated EMT, suggesting regulation of different signaling pathways in different subtypes of breast cancer11,13. The pathways of small GTPases and actin filament organization, which control cellular movement, together with the enriched functions of mesenchymal differentiation, cell projection, and extracellular matrix organization, link the antimetastatic effects of the receptor in IBC with a more physiological role in tissue differentiation during mammary gland development and after its massive expansion by the stimulation of ERα signaling (Fig. 7D)57. Although this role is evident from the analysis of our cancer cells, the lack of appropriate transgenic models to pinpoint its expression and function in physiological tissues has hindered progress in defining its biology in normal gland and cancer58.
In addition to the anticipated gene ontology terms, our extended pathway analysis revealed biological processes that were not previously attributed to ERβ. An increased number of enriched pathways are involved in hormone-associated metabolism, primarily in phospholipid metabolism. It is not clear whether these pathways mediate effects of the receptor specifically in highly metastatic IBC or represent a more general mechanism that is employed to adjust cancer metabolism to tumor differentiation. These questions call for additional investigations to determine whether alterations in fatty acid metabolism are required for the observed changes in cytoskeleton re-organization and metastasis or represent fundamental biochemical adaptations in cellular phenotypes with altered estrogen signaling. It is also equally important to explore how this altered ERβ signaling impacts the tumor microenvironment and how cellular interactions in this environment influence cancer progression and therapy response.
We identified several ERβ-repressed genes that function in oncogenic and metastatic pathways in IBC and verified the involvement of some of these genes in ERβ effects, including those that promote actin-based cell migration (ELMO1, GPR141, PREX1) and fatty acid metabolism (FABP5, CAV1, FASN)15,59. Among these genes are druggable molecules (FASN, FABP5) with their inhibitors already clinically tested in other conditions, including non-IBC tumors and other potentially viable therapeutic targets due to their expression on the cell surface (GPR141, CAV1)60,61. A comprehensive functional characterization of the identified ERβ target genes can assist their evaluation as molecular targets and, thus, determine whether their inhibition alone, or in combination with ERβ agonists, represents a novel strategy to effectively treat IBC.
Our findings, by linking cancer cell differentiation to endogenous ERβ, offer additional confidence about its presence and its proposed tumor repressive role in breast cancer. They also help uncover unknown aspects of the mechanism of action and the downstream signaling, which are essential for discovering novel viable molecular targets that, through more accurate calibration of estrogen signaling, will benefit the patients.

Methods

Cell culture
The KPL4 cells with different ERβ levels (Control, ERβKO and ERβΚΟ + ERβ) were developed and characterized as previously described62. Control and ERβ-expressing SUM149 cells were generated by infection with lentivirus containing the pLenti6/V5-FLAG-ERβ1 (full-length ERβ) construct as previously published11. Cells were maintained in HAM’s-F12 media (Corning) supplemented with 7% fetal bovine serum (FBS) (Sigma), 10 mM HEPES (Thermo Fisher), 5 μg/mL insulin (Sigma), 1 μg/mL hydrocortisone (Sigma) at 37 °C in 5% CO2. For ligand treatments, cells were cultured in phenol red-free media supplemented with 1% dextran-coated charcoal-treated FBS (Thermo Fisher) for 72 h prior to treatment with the ERβ agonist LY500307 (ApexBio) or DMSO (Sigma) for 2 h for ChIP-Seq οr for 24 h for RNA-Seq in the same media. All the cell lines were tested negative for mycoplasma, authenticated by short tandem repeat (STR) profiling and used within three passages after thawing. All reagents for cell culture and ligands are listed in Supplementary Data 1.

Immunoblotting
Immunoblotting was performed as previously described11. Briefly, cells were lysed using RIPA buffer (Thermo Fisher) with protease and phosphatase inhibitors (Roche, Sigma). Following separation in SDS-PAGE, proteins were transferred onto a nitrocellulose membrane (Thermo Fisher). Membranes were blocked with 5% dry milk and probed with the ERβ antibody (PPZ0506, Thermo Fisher) overnight at 4 °C and anti-β-actin antibody (Abcam). Proteins were visualized with an ECL detection kit (GE Healthcare).

Migration assay
For the wound-healing assay, 3 × 104 cells were seeded in a 96-well plate and allowed to form a monolayer before being treated with siRNAs (Supplementary Data 2). The wound was formed by the IncuCyte woundmaker, and cell density inside and outside the wound was imaged in an IncuCyte S3 Live-Cell Analysis System (Sartorius) constantly for up to 3 days and evaluated using the Scratch Wound Analysis Software.

RhoC activation assay
RhoC activity was assessed as previously described15. Briefly, cells were lysed in buffer containing 25 mM HEPES, pH 7.5, 150 mM NaCl, 1% NP-40, 10 mM MgCl2, 1 mM EDTA, 2% Glycerol and protease inhibitors for 20 min on ice. Lysates were centrifuged for 10 min at 14,000 × g at 4 °C and supernatants were incubated with Rhotekin RBD agarose beads overnight at 4 °C. After washing, proteins bound to beads were dissociated via resuspension in 2X SDS-PAGE sample buffer and analyzed by immunoblotting.

RNA extraction, real-time quantitative reverse transcription (RT) PCR
Total RNA was isolated using the Qiagen RNeasy Mini kit (Qiagen) and quantified using the NanoDrop 1000 spectrophotometer (Thermo Fisher). The quality of RNA was evaluated by capillary electrophoresis on an Agilent 2100 Bioanalyzer (Agilent Technologies). For qPCR, the mRNA was reversed-transcribed to cDNA using the High Capacity cDNA Reverse Transcription kit (Applied Biosystems). Real-time PCR was performed with the TaqMan gene expression master mix (Thermo Fisher) and gene-specific TaqMan assays (Supplementary Data 2) on a QuantStudio 6 Flex (Applied Biosystems) using the QuantStudio™ Real-Time PCR software. For miRNA analysis, miRNA was reverse transcribed from 20 ng total RNA using the TaqMan™ MicroRNA Reverse Transcription Kit (Thermo Fisher) and miRNA-specific RT-PCR primers (Supplementary Data 2). Quantitative real-time PCR was then performed in triplicates using miRNA-specific TaqMan assays (Supplementary Data 2) and TaqMan™ Universal PCR Master Mix (no AmpErase™ UNG). All quantitative data were normalized to GAPDH and ACTIN-β for mRNA and RNU6B for miRNA and the threshold cycle (Ct) was used for quantification.

Chromatin immunoprecipitation
Chromatin immunoprecipitation was performed as previously described15. Briefly, Protein-DNA complexes were captured by crosslinking with 11% formaldehyde for 10 min at room temperature. Nuclei were purified from cell pellets following washes in PBS, snap-freezing in liquid nitrogen and consecutive treatments with ice-cold hypotonic buffer (10 min at 4 °C) and NP-40/Triton X-100-based lysis buffer with protease inhibitors. Nuclei were then washed in ice-cold wash buffer, and chromatin was sheared by sonication and clarified by centrifugation. During immunoprecipitation, chromatin-bound FLAG-tagged ERβ was incubated with monoclonal anti-FLAG M2 antibody (Sigma) or ERβ antibody (PPZ0506, Thermo Fisher) and mouse IgG (Santa Cruz) and precipitated with protein G magnetic beads (Thermo Fisher). Following elution, protein-DNA complexes were decrosslinked at 65 °C overnight, and after treatment with RNase and proteinase K, the DNA was purified, and its enrichment was measured by real-time PCR using specific primers flanking the ERβ-binding sites (Supplementary Data 2).

Promoter analysis and luciferase assays
ERE-containing DNA fragments that are located upstream and within the promoter of SERPINA1, HSPB8, TRIM31 and TFF1 TSS were PCR amplified from human genomic DNA with forward and reverse primers that are listed in Supplementary Data 2. Amplicons were cloned between Kpn1 and Nhe1 sites into pGL4.10 vector (Promega). Cells were transfected using Lipofectamine 2000 with the reporter plasmids and a Renilla luciferase construct to normalize for transfection efficiency. Luciferase activity was measured with Dual Luciferase Reporter Assay (Promega) on a Luminometer (Glomax Navigator, Promega) as previously described63. Relative light units were calculated as the ratio of firefly luciferase to renilla luciferase activity.

Preparation of ChIP-sequencing libraries
Libraries for the ChIP DNA were prepared and sequenced by the Wistar Institute Genomics Core Facility. ChIP libraries were prepared using the NEBNext Ultra II DNA Library Prep Kit for Illumina (NEB #E7645). Quality control of the final libraries was performed on Agilent 2100 Bioanalyzer (Agilent Technologies), and the libraries were sequenced on NextSeq500 (Illumina). Sequencing of the ChIP-Seq libraries was performed, producing 151 bp single-end reads.

Preparation of total RNA and miRNA sequencing libraries
Libraries for total RNA and miRNA were prepared using the Illumina Stranded Total RNA Prep Ligation with Ribo-Zero Plus (Illumina) and the QIAseq miRNA Library Prep Kit for Illumina (QIAGEN), respectively. Quality control of the final libraries was performed on an Agilent 2100 Bioanalyzer, and the total and miRNA-seq were carried out on NextSeq500 (Illumina) in the Wistar Institute Genomics Core Facility. For the total RNA, at least 50 million paired-end reads of 74 bp in length were generated per sample using the High Output Kit. For the miRNAs, sequencing yielded ~15 million 75 bp single-end reads per sample using the Mid Output Kit.

RNA-seq data analysis
Quality control, quality trimming and adapter removal were performed with the combined use of FASTQC, Trim Galore and Minion tools. Pre-processed RNA-seq reads were aligned to the GRCh38 assembly with Ensembl v104 annotation, using GMAP aligner parameterized for splice-aware execution. The RSEM quantification tool was used to derive relative abundance estimates. Gene-level differential expression (DE) analysis was performed in R with the EdgeR package64. Briefly, weakly expressed genes were removed using the FilterByExpr function in R-imported raw read count matrices. Following the calculation of effective library sizes and dispersions, differential expression was assessed in pairwise contrasts with Quasi-Likelihood F-testing. Significantly differentially expressed genes were defined as those exhibiting an absolute log2(Fold-Change) (log2FC) ≥ 0.585 (i.e., ≥1.5-fold difference) and p (Benjamini-Hochberg adjusted, FDR) < 0.05. Package pheatmap65 was employed to generate heatmaps of z-score-normalized log2(CPM) moderated values64 for contrasts and their intersection. For miRNAs, small RNA-seq libraries were subjected to automatic quality control, pre-processing and quantification for miRNAs using DIANA-mAP workflow that utilizes FASTQC, Trim-Galore, Minion and DNApi for adapter-contaminant inference and removal and quality trimming, and performs quantification using miRDeep2 against precursor/mature miRNA sequences from miRBase v22.1. Differential expression analysis of mature miRNAs was performed as in long RNAs.

ChIP-seq data analysis
Samples were controlled for quality with FASTQC and processed with Trim-Galore and Minion to remove adapter sequences, retaining reads longer than 25 nt. Genomic alignment (GRCh38 assembly) was performed with bowtie2 at default settings66. Alignment files were processed with samtools and bedtools to retain uniquely mapped reads and filter out problematic “ENCODE blacklist” regions67. Peaks were called using MACS2 (FDR < 0.05) for individual pairs of samples, using deepTools’ plotEnrichemnt to obtain Fraction of Reads in Peaks (FriP). Replicate MACS2 results were combined with MSPC with default parameters (-r bio -w 1e-4 -s 1e-8)68,69. Briefly, Fisher’s meta-analysis was applied in overlapping peaks between replicates to obtain combined (FDR-corrected) significance levels for each peak, retaining (a) peaks exhibiting strong significance (FDR < 10−8) even in one replicate, and (b) weaker peaks (initial MACS2 FDR < 10−4) that occur in more than one replicate and achieve a combined FDR < 10−8. Quality control metrics are provided in Supplementary Table 1. For the significant peaks, counts of uniquely aligned reads were obtained with bedtools. Package pheatmap was used to create heatmaps of z-score-normalized log2(counts+1) for experimental and input libraries after subtracting the reads from the control or vehicle-treated KPL4 and SUM149 cells in each replicate using the EdgeR’s function calcNormOffsetsforChIP to normalize each treatment replicate with its corresponding control. Significant peaks were annotated using Ensembl v104, by defining promoters as the regions from 1 kb upstream to 100 bp downstream of the TSS, and gene bodies as the regions from 101 bp downstream of TSS until the gene end. Gene body peaks were further classified according to whether they reside in exonic or intronic regions, while peaks not overlapping promoters and gene bodies were classified as intergenic. Markov models of background sequences that do not overlap significant peaks were created with MEME suite (v5.4.1), separately for loci overlapping (i) promoters, (ii) gene bodies, (iii) regions from 1 mb to 1 kb upstream of TSS and (iv) regions up to 1 mb downstream of gene ends. FIMO70 was used to scan for significant individual motif occurrences (p < 10−4), using transcription factor binding motifs from JASPAR 2020 vertebrate collection.

RNA-seq and ChIP-seq integration
To identify direct targets of ERβ, significantly differentially expressed mRNAs and miRNAs were overlapped with ChIP-seq peaks and FIMO-derived motifs residing in promoters (−1 kb to +100 bp and −4 kb to +100 bp of the TSS for mRNAs and miRNAs, respectively), gene bodies, loci upstream of promoters and downstream of gene ends. The TSSs for miRNAs were derived from integrating the information from the available Cap analysis of gene expression (CAGE) analysis of MDA-MB-453 breast cancer cells with the data of intragenic miRNAs considering the TSSs of overlapping mRNAs29. The intersection of RNA-Seq with ChIP-seq data was performed for the contrasts Control vs. ERβΚΟ, ERβKO vs. ERβKO + ERβ and their intersection (concordantly up-/down-regulated genes) for mRNAs. For miRNAs, the intersection was performed for the contrast ERβKO vs. ERβKO + ERβ. Heatmaps of significantly deregulated genes overlapping significant ChIP-seq peaks were created using pheatmap package. ERβ direct target genes were verified by performing integration of ChIP-seq peaks with differentially expressed genes with the Binding and Expression Target Analysis (BETA) method using <0.01 as Rank Product Threshold and <0.05 as FDR threshold for differentially expressed genes.

Term enrichment analysis
Over-representation analyses were conducted against KEGG pathways and Gene Ontology terms (Molecular Functions, Biological Processes and Cellular Components) separately. Briefly, via ClusterProfiler R package71, different sets of genes were subjected to hypergeometric testing, using Ensembl Gene Identifiers and setting the significance threshold at FDR < 0.05.

Analysis of clinical datasets
Analysis of clinical datasets has been performed in accordance with our IRB protocols (1) Houston Methodist Research Institute “Evaluating the prognostic significance of ERβ in breast cancer” (RN:PRO00039865) and (2) MD Anderson Cancer Center “ERβ in IBC” (RN: 2019-1140). MAS5-level log2-transformed microarray (Affymetrix HG-U133A chip) expression data were obtained from a cohort of IBC (25) and non-IBC (58) patients that received neoadjuvant chemotherapy (anthracycline- and taxane-based chemotherapy) followed by surgery at MD Anderson Cancer Center. IBC patients had a similar histological classification (12% grade II and 88% grade III) with non-IBC patients (22.5% grade II and 77.5 grade III). Probesets were annotated to recent Ensembl Gene Identifiers and Gene Symbols using the annotation R package “hgu133a.db”. Survival analyses were conducted separately for IBC and non-IBC patient groups. For each gene, Kaplan-Meier analysis was conducted by separating “low” and “high” expression subgroups, using Maximally Selected Rank Statistics (MSRS), as implemented in the R package “survminer”, to obtain the optimal cutoff value. Briefly, function surv_cutpoint was employed to test all possible cutoffs, allowing a minimum 30% proportion of observations per group (miniprop = 0.3), and to perform log-rank testing for each cutoff. The value that produced the maximum test statistic (i.e., yielded the best separation between subgroups) was chosen as the optimal cutoff point. The Log-rank test was used to assess the null hypothesis of no difference in the survival distribution of the “low” and “high” expression subgroups, and respective Kaplan-Meier plots were generated (R packages “survival” and “ggplot2”). Benjamini-Hochberg adjustment was used to control the FDR at the 0.05/0.1 levels.

Supplementary information