IL6/STAT3 Signaling Hijacks Estrogen Receptor α Enhancers to Drive Breast Cancer Metastasis

Abstract

The cytokine interleukin-6 (IL6) and its downstream effector STAT3 constitute a key oncogenic pathway, which has been thought to be functionally connected to estrogen receptor α (ER) in breast cancer. We demonstrate that IL6/STAT3 signaling drives metastasis in ER⁺ breast cancer independent of ER. STAT3 hijacks a subset of ER enhancers to drive a distinct transcriptional program. Although these enhancers are shared by both STAT3 and ER, IL6/STAT3 activity is refractory to standard ER-targeted therapies. Instead, inhibition of STAT3 activity using the JAK inhibitor ruxolitinib decreases breast cancer invasion in vivo. Therefore, IL6/STAT3 and ER oncogenic pathways are functionally decoupled, highlighting the potential of IL6/STAT3-targeted therapies in ER⁺ breast cancer.

Keywords: FOXA1; IL6; STAT3; breast cancer; estrogen receptor; metastasis; mouse intraductal xenograft model; pioneer factor.

Figure 1. IL6/STAT3 signaling increases primary tumor growth and metastasis.

(A) Representative H&E staining and immunohistochemistry (IHC) for phosphorylated STAT3 (pSTAT3^Y705) and ER in luciferase-expressing T47D primary MIND tumors expressing human IL6 or empty vector (EV). Scale bar = 50 μm. (B) Growth of the primary T47D MIND tumors as determined by changes in radiance relative to day 1 after cell injection. 5 mice per condition, 17 glands for EV and 15 glands for IL6. Mean of the fold change in radiance per gland is shown with error bars representing SEM. ****padj<0.0001 (two way ANOVA, Sidak’s multiple comparisons test). (C) Relative metastatic burden in mice having primary tumors expressing IL6 or EV (n=5) as measured by radiance. Error bars represent SEM. ***padj<0.001 (two-way ANOVA, Sidak’s multiple comparisons test). (D) Representative images showing radiance from micrometastases in bones and lungs used in panel C. (E,F) Kaplan-Meier (KM) plots showing the association between BCSS (E) or OS (F) and IHC staining of pSTAT3^Y705 and IL6 in breast tumors from the Nottingham cohort (Ahmad et al., 2018; Aleskandarany et al., 2016). Log-rank p values are indicated for IL6⁺/pSTAT3⁺ vs all other cases. See also Figure S1.

Figure 2. STAT3 associates with the ER/FOXA1 complex in response to activation by phosphorylation in cell lines and patients.

(A) ER and STAT3 qPLEX-RIME in MCF7 cells in response to 30 min treatment with recombinant human IL6. Four replicates of the ER or STAT3 RIME and one pooled IgG control RIME were included in each 9plex tandem mass tags experiment. Red points, padj<0.05, qPLEXanalyzer (Papachristou et al., 2018). (B) ER qPLEX-RIME in nine patient samples, i.e. five primary ER⁺ breast tumors and four ER⁺ pleural effusions. Non-specific proteins previously identified in clinical breast tumors (Papachristou et al., 2018) were removed. In addition to ER (red) and STAT3 (blue), known ER-associated proteins are highlighted in light blue. The peptide coverage of STAT3 is shown in the insert. The schematic illustration was created with BioRender.com. See also Figure S2 and Table S1.

Figure 3. STAT3 chromatin binding overlaps extensively with binding of ER and FOXA1 in cell lines and patients.

(A) ChIP-seq in T47D cells in response to 30 min IL6 treatment (n=4). Groups of binding sites are defined based on the positional overlap between pSTAT3 and ER binding sites in the IL6 condition. Color bars indicate scale for normalized tag densities per 20 M subsampled reads. (B) Shared pSTAT3-ER binding sites from panel A are divided into constitutive or increased (padj<0.05, DiffBind (Stark and Brown, 2011)) ER binding sites. UCSC genome browser (http://genome.ucsc.edu, Dec. 2013 GRCh38/hg38) (Kent et al., 2002) screen shots for a representative gained and common ER binding site are shown to the right. (C) Positional overlap between ER and pSTAT3 chromatin binding based on ChIP-seq in ER⁺/pSTAT3⁺ primary breast tumors. The heatmaps show binding in the shared pSTAT3-ER sites for each patient. Color bars indicate scale for normalized tag densities per 20 M subsampled reads. The schematic illustration was created with BioRender.com. See also Figures S3 and S4.

Figure 4. Shared pSTAT3-ER-FOXA1 sites control the IL6-induced gene program.

(A) Left, differentially expressed genes (padj<0.05, DESeq2 (Love et al., 2014)) in response to one hour IL6 treatment of T47D cells as determined by RNA-seq (n=6). Right, significantly enriched Hallmark terms from MSigDB (Liberzon et al., 2015) in IL6-induced genes. (B) Enrichment of binding sites from T47D cells and patients in the vicinity of IL6-induced genes compared to an equal number of constitutively expressed genes in T47D cells. Lines illustrate the cumulative percentage of sites within a given distance from the TSS. (C) Percentage of breast cancer ATAC-seq peaks occupied by both pSTAT3 and ER or only pSTAT3 in T47D cells that are connected to IL6 target gene promoters through promoter-enhancer connections (Corces et al., 2018) (Red line). The distribution of connections between the same number of random ATAC-seq peaks and IL6 target gene promoters is shown for comparison in blue (10,000 iterations to estimate the distribution of the test statistic under the null hypothesis). *Approximate random permutation p value <0.0001. (D) Identified functional promoter-enhancer links between IL6-induced genes and pSTAT3-ER occupied breast cancer ATAC-seq peaks. See also Figure S5.

Figure 5. IL6/STAT3 signaling through pSTAT3-ER-FOXA1 enhancers is functionally independent of ER/FOXA1.

(A) Venn diagram showing the overlap between IL6-dependent (1,282 sites) and FOXA1-dependent (2,303 sites) ATAC-seq peaks (n=3). IL6-dependent sites were defined as those with padj<0.05 (DiffBind (Stark and Brown, 2011)) comparing siControl-IL6 to siControl-veh. FOXA1-dependent sites were defined as those with padj<0.05 comparing siFOXA1-IL6 to siControl-IL6. (B) ATAC-seq and ChIP-seq for pSTAT3, ER, FOXA1 and H3K27ac in IL6-and FOXA1-dependent sites from panel A (n=3-4). FOXA1 ChIP-seq is from Figure 3A. (C) Enrichment of IL6- and FOXA1-dependent sites from panel A in the vicinity of IL6-induced genes compared to constitutively expressed genes. Lines illustrate the cumulative percentage of sites within a given distance from the TSS. (D) Gene expression as reads per kb per million (RPKM) for IL6-induced genes with TSS within 50 kb of IL6-dependent ATAC-seq peaks (113 genes) in response to IL6 and STAT3, ER or FOXA1 knock down (n=4). *p<10^-6, **p<10^-9 (two-sided Wilcoxon test comparing veh and IL6 treatment for the different knock downs). Boxplot definitions: Center line, median; box limits, upper and lower quartiles; whiskers, 1.5x interquartile range; any outliers have been hidden. Notch shows the approximate 95% confidence interval for the median. (E) Regulation of IL6-induced genes (1 h treatment, padj<0.1, DESeq2 (Love et al., 2014)) in STAT3^-/- MCF7 cells and in response to fulvestrant treatment (100 nM for 3 h prior to the start of the IL6 treatment) in WT cells (n=4). (F) Expression of FOXA1-induced genes with TSS within 50 kb of FOXA1-dependent ATAC-seq peaks (369 genes) as in panel D. *p<10^-4, **p=2.2x10^-308 (two-sided Wilcoxon test comparing veh and IL6 in siControl to veh and IL6 in the different knock downs, respectively). See also Figure S6.

Figure 6. IL6/STAT3 signaling is resistant to endocrine therapy, but the JAK inhibitor ruxolitinib inhibits activation of STAT3 and invasion in resistant PDX.

(A) Regulation of Hallmark early estrogen responsive genes (n=200) (Liberzon et al., 2015) and IL6-induced genes (n=328 for T47D and n=67 for MCF7) in matched samples from 42 patients treated with letrozole for at least four months in the neoadjuvant setting (Selli et al., 2019). The mean regulation of 200 random genes is included for comparison with the error bar representing the standard deviation of 1,000 iterations. *p=1.8e-07 and **p< 2.2e-16 (Fisher’s exact test, two-sided). The schematic illustration was created with BioRender.com. (B) Representative H&E, pSTAT3^Y705 and Ki67 IHC staining for a WT ER⁺/PR^- PDX (AB555P) (Bruna et al., 2016) explant treated with ruxolitinib (500 nM) for two days ex vivo. Scale bars = 50 μm. (C) Quantification of pSTAT3^Y705 and Ki67 IHC staining in the PDX explant from panel B using the Halo TM (Indica labs) Image analysis platform. Error bars represent SEM between different tumor pieces (n=9-10) from different parts of the PDX tumor cultured on the same sponge. *p=0.099, **p=1.1e-06 (Student’s t-test, two-sided). (D) Relative mRNA expression of MYC, ELF3 and BCL3 in the PDX explant from panel B. Average dots per cell for these three genes and the housekeeping gene TATA-binding protein (TBP) were determined by RNAscope in situ hybridization. Error bars represent SEM between different tumor pieces (n=7-8) from different parts of the PDX tumor cultured on the same sponge. *p= 0.098, **p=0.0074, ***p=0.00013 (Student’s t-test, two-sided). (E) Primary tumor growth of endocrine-resistant intraductal PDX in response to fulvestrant (Fulv) or ruxolitinib (Ruxo) treatment (initiated one week after injection). Points represent average growth of tumors and error bars represent SEM. n=19, 22, 22, 16 tumors from 8 mice per treatment arm. *p<0.05, ****p<0.0001 (Two way ANOVA, followed by Tukey’s multiple comparison test). Insert shows the relative fold change in radiance at endpoint. (F) Quantification of expression of protein markers in the PDX from panel E as determined by IHC. Error bars represent SEM, n=4-5. (G) H&E staining and representative images of protein marker expression in adjacent sections in the PDX from panel E. (H) Quantification of Ki67 expression from panel G. Error bars represent SEM, n=4-5. (I) Quantification of tumor invasion in the mammary gland in the PDX from panel E. The invasive areas of the tumor, characterized by loss of the myoepithelial layer of the duct and tumor cell invasion into the surrounding tissue, were quantified by automated cell detection, manually curated and annotated. The percentage area of invasion relative to total tumor area (invasive and in situ) was determined. Error bars represent SEM, n=4-5. The presence of one outlier in the fulvestrant arm was identified and removed using Graphpad Prism outlier function with Q=2%. *p<0.05, **p<0.01, ***p<0.001, ****p<0.0001 (One way ANOVA, followed by Tukey’s multiple comparison test). See also Figure S7.

Figure 7. Schematic model of IL6/STAT3 and ER/FOXA1 function in breast cancer.

The two potent oncogenic pathways, ER/FOXA1 and IL6/STAT3, drive breast tumor growth and metastasis through functionally independent mechanisms using shared enhancers. Created partly with BioRender.com. FHRE, Forkhead response element; ERE, Estrogen response element; STAT-RE, STAT response element.