Estrogen regulates divergent transcriptional and epigenetic cell states in breast cancer

Abstract

Breast cancers are known to be driven by the transcription factor estrogen receptor and its ligand estrogen. While the receptor's cis-binding elements are known to vary between tumors, heterogeneity of hormone signaling at a single-cell level is unknown. In this study, we systematically tracked estrogen response across time at a single-cell level in multiple cell line and organoid models. To accurately model these changes, we developed a computational tool (TITAN) that quantifies signaling gradients in single-cell datasets. Using this approach, we found that gene expression response to estrogen is non-uniform, with distinct cell groups expressing divergent transcriptional networks. Pathway analysis suggested the two most distinct signatures are driven separately by ER and FOXM1. We observed that FOXM1 was indeed activated by phosphorylation upon estrogen stimulation and silencing of FOXM1 attenuated the relevant gene signature. Analysis of scRNA-seq data from patient samples confirmed the existence of these divergent cell groups, with the FOXM1 signature predominantly found in ER negative cells. Further, multi-omic single-cell experiments indicated that the different cell groups have distinct chromatin accessibility states. Our results provide a comprehensive insight into ER biology at the single-cell level and potential therapeutic strategies to mitigate resistance to therapy.

Figure 1.

TITAN, a topic modeling approach to scRNA-seq analysis, identifies distinct gene sets (topics) specific to estrogen signaling gradients. (A) Schematic representation of experimental setup and TITAN: Topic Inference of Transcriptionally Associated Networks in scRNA-seq. Breast cancer cell and PdXO models are subject to 100 nmol/l β-estradiol and/or 100 nmol/l progesterone stimulation, triggering signaling network changes across cells. Using Latent Dirichlet Allocation (LDA), TITAN identifies latent transcriptional networks, or topics, in the scRNA-seq dataset. Topics are linked by a distribution of scores that relate genes to topics and topics to cells. These scores are then used to infer transcriptional networks. TITAN identifies distinct gene gradients in comparison to standard clustering approaches which are binary distributions. (B) PCA-based (top) and TITAN-based (bottom) UMAP visualization of MCF-7 cells treated with estrogen for different time points. (C) Comparison of Z-scores normalized to control in time-dependent estrogen response networks for TITAN, PCA-based clustering (Seurat), SCENIC and CountClust.

Figure 2.

Divergent estrogen driven topics identified in breast cancer cell lines and organoids. (A) Heatmap of TITAN normalized cell topic scores for MCF-7 cells treated with estrogen for different time points. (B) Feature plot representation of Topic 8, representative of topics increased by longer estrogen treatments and Topic 13, representative of topics decreasing with longer estrogen treatments. (C) Scatter plot representation of T-47D cells expressing topic 20 and topic 4 per duration of estrogen treatment. (D) Scatter plot representation of topic 20, with expression levels of TFF3 (red) and topic 4 with expression levels of HMMR (blue), genes associated specifically with each respective topic. (E) UMAP dimensionality reduction visualization of MCF-7 (14,788 cells), T-47D (10,424 cells), ZR-75–1 (9,433 cells), TAMR (3,750 cells), HCI003 (1,615 cells) and HCI011 (13,470 cells) breast cancer models colored by estrogen treatment duration (F) UMAP dimensionality reduction visualization of the same cells as (E) by expression of FOXM1 (purple) and ESR1 (orange) driven topics. (G) Ridge-plot representation of ESR1 (red) and FOXM1 (blue) topic expression distribution in breast cancer models treated with estrogen for 48 h.

Figure 3.

Estrogen and progesterone induce both unique and shared gene expression signatures. (A) PCA based UMAP dimensionality reduction of T-47D (10,668 cells) treated with estrogen (E2), progesterone (PG) alone or in combination (E2 + PG) for 3 and 48 h, colored by treatment. (B) Heatmap visualization of TITAN normalized cell topic scores for T-47D of the same dataset. (C) Scatter plot representation of topic 1 (y-axis) which represents networks regulated by progesterone treatment compared to topic 17 (x-axis) which represents a shared gene set activated by both hormones colored by type of treatment. (D) Violin plot representation of topic 17 (left) representative of gene sets activated by both hormones, 11 (center) representative of gene sets activated primarily by estrogen and 1 (right), representative of gene sets activated by progesterone. (E) Transcription factor (TF) enrichment analysis of genes involved in topics 17 (left), 11 (center) and 1 (right) using publicly available ChIP datasets from CistromeDB (29,30), significance of overlap P-values is generated by hypergeometric test.

Figure 4.

FOXM1 is activated through phosphorylation and regulates corresponding topics upon estrogen stimulation. (A) Distribution of pFOXM1 protein levels quantified as mean fluorescence intensity within 37,911 nuclei in MCF-7 cells plotted by cell cycle phase and estrogen (20 nmol/l) treatment length imaged by high-content quantitative imaging. Boxplot hinges correspond to the 25th–75th and whiskers correspond to 1.5 × IQR (inter-quartile range) of the hinge. Outlying points are plotted individually. (B) Relative gene expression of FOXM1 targets in MCF-7 cells with 20 nmol/l estrogen treatment for 6 h. SDH gene used as internal control. Error bars represent s.e.m. in three biological replicates.

Figure 5.

Divergent ESR1 and FOXM1 topics identified in breast cancer patient scRNA-seq datasets and show subtype preference. (A) Ridge-plot representation of ESR1 (left) and FOXM1 (right) topic expression distribution in 36,798 epithelial cells from 20 breast cancer patients colored by tumor subtype. (B) Scatter plot representation of the ESR1 topic 36 (y axis) compared to the FOXM1 topic 39 (x axis) colored by tumor subtype in 20 breast cancer samples. Data used in this figure are from Wu et al. (13).

Figure 6.

Multiomic analysis reveals divergent gene expression and chromatin states upon estrogen stimulation. (A) Schematic of experimental design and analysis of multiome single cell profiles in MCF-7 and T-47D cell lines with and without 20 nmol/l β-estradiol (E2) treatment. RNA profiles undergo TITAN topic modeling and paired ATAC profiles undergo cisTopic topic modeling. Results are then correlated. (B) UMAP projection of MCF-7 (left, 6,438 cells) cell profiles colored by treatment. Projections were calculated using TITAN, cisTopic, and TITAN + cisTopic profiles (see Methods). UMAP projection of TITAN + cisTopic profiles colored by ESR1 and FOXM1 defined TITAN topics (far right). Color legend shared with panel (D). (C) Topic expression heatmap correlation between TITAN and cisTopic topics. TITAN ESR1 (Topic 37, blue rectangle) and FOXM1 (Topic 4, red rectangle) topics are highlighted. UMAP projection of TITAN + cisTopic profiles colored by cisTopic topics 20 (left) and 13 (right) showing open chromatin regions specific to TITAN topics (bottom). (D) Scatterplot of FOXM1 and ESR1 TITAN topic values in MCF-7. Main figures are colored by treatment with R2 correlation and linear fit shown as a black line. Inset panels show binned FOXM1 and ESR1 topic enriched cells. These bins were defined as cells having ≥20% quantile value for the respective bin and ≤75% quantile value of the opposite bin. This resulted in 1127 and 1,152 MCF-7 cells for FOXM1 and ESR1 topics respectively. Color legend shared with panel (C). (E) Volcano plot of differential expression and cistrome defined enrichment scores of ChIP-seq peaks between FOXM1 (right) and ESR1 (left) topics for binned MCF-7 cells defined in panel (D). The top 10 genes or CistromeDB peak sets are shown for FOXM1 and ESR1 bins. Names displayed show genes or target protein and cistrome ID. Positive Log2 fold-change reflects ESR1 bin enrichment. Points plotted in red are significant (adjusted P-value ≤ 0.05). (F) Coverage plot of aggregate ATAC profiles for binned MCF-7 cells defined in panel (D). Genome track covers the gene body of PGR (left) and CENPF (right) with 5kbp up and downstream. Black arrow represents the transcription direction. Side plot displays the differential gene expression of PGR. Peaks subpanel shows defined open regions on aggregate data, peaks are then colored by overlap with FOXM1 and ESR1 binding motifs. Links subpanel displays cis-coaccessibility linked between peaks, displaying peaks which are correlated. Gray highlight boxes overlap peaks which show at least nominal significance (logistic regression, P-value ≤ 0.05). *** denotes adjusted logistic regression P-value ≤ 0.01.