Single-cell RNA-sequencing analysis of estrogen- and endocrine-disrupting chemical-induced reorganization of mouse mammary gland

Abstract

Menopause is a critical window of susceptibility for its sensitivity to endocrine disrupting chemicals due to the decline of endogenous estrogen. Using a surgical menopausal (ovariectomized) mouse model, we assessed how mammary tissue was affected by both 17β-estradiol (E2) and polybrominated diphenyl ethers (PBDEs). As flame retardants in household products, PBDEs are widely detected in human serum. During physiologically-relevant exposure to E2, PBDEs enhanced E2-mediated regrowth of mammary glands with terminal end bud-like structures. Analysis of mammary gland RNA revealed that PBDEs both augmented E2-facilitated gene expression and modulated immune regulation. Through single-cell RNA sequencing (scRNAseq) analysis, E2 was found to induce Pgr expression in both Esr1⁺ and Esr1^- luminal epithelial cells and Ccl2 expression in Esr1⁺ fibroblasts. PBDEs promote the E2-AREG-EGFR-M2 macrophage pathway. Our findings support that E2 + PBDE increases the risk of developing breast cancer through the expansion of estrogen-responsive luminal epithelial cells and immune modulation.

Keywords: Breast cancer; Computational biology and bioinformatics; Endocrinology.

Fig. 1

Preliminary whole-mount and histological analysis showed differences between vehicle, E2, and E2 + PBDE-treated mice. a Overall scheme of experimental design. BALB/cj female mice (9 weeks old) were ovariectomized. Mice were treated (n = 10 per treatment) at 19 weeks of age and were euthanized one week after treatment. b Toluidine blue stained whole-mounted mammary glands showed that vehicle mice had regressed branching trees as characterized by thin ductal structure by comparing to mammary gland structure of intact mice. Treatment of OVX mice with E2 or E2 + PBDE showed a regrowth of gland structures, characterized by thicker and larger ductal structures and greater prevalence of TEB-like structures. Scale bars = 0.4 cm. c TEB-like structures had multiple layers of cells with Ki67-positive staining. d The TEB count was significantly increased in mice fed with E2 + PBDE compared with the mice with E2 only. The data are expressed as mean ± SEM of the mean using duplicate assays. Tukey’s multiple comparisons test was performed. *P < 0.05. e Duct quantification by KRT18 immunohistochemistry staining. The number of large, small, and total ducts were counted among vehicle treated (n = 4), PBDE treated (n = 6), E2 treated (n = 5), and E2 + PBDE treated (n = 5) mice. A threshold of 1906 µm² was used to classify ducts as large or small. Ducts that fell below this threshold were classified as small, and ducts that had areas larger than this value were classified as large. f Overview of scRNAseq approach using mouse mammary gland samples. Cells from two separate experiments were sequenced on the 10X Genomics platform and then pooled together for analysis

Fig. 2

ScRNAseq yielded 11 different clusters representing different populations of luminal, fibroblast, immune, and other cells. a T-SNE plot of cells from the combination of two independent experiment sets from vehicle, E2 and E2 + PBDE group. Cells were color coded according to clusters. b Heatmap showing specific genes (top 5 differentially expressed genes in each cluster, top 20 gene lists shown in Supplementary Data 1). c Dendogram showing cluster relationship to other clusters based on gene expression data. d Violin plots of gene expression of selected genes unique for each cluster. Krt18 (luminal cell marker), Krt14 (basal cell marker), Ptprc (immune cell marker), Col1a1 (fibroblast marker), Pecam1 (endothelial cell marker), were used to identify distinct cell populations. e T-SNE plots of cells separated by treatment group. f Cell distributions from two independent experiments by percentage for all three treatment groups among all 11 total clusters. Symbols indicate individual cell distributions for each of the two experimental replicates (triangle, circle)

Fig. 3

Luminal cells were re-clustered to reveal three subpopulations. a Feature plot highlighting clusters C4 and C5, which were putatively identified as luminal epithelial cells based on established epithelial markers, such as Krt18 and Epcam. b T-SNE plot of re-clustered luminal cells shows four distinct clusters. c Heatmap showing specific genes (top 5 differentially expressed genes in each luminal cluster). d Dendrogram showing luminal cluster relationship to other luminal clusters based on gene expression data. e Selected feature plots of genes used to classify the subpopulations of luminal cells into luminal progenitor cells, intermediate luminal cells, and mature luminal cells. f Heatmap on the MSigDB analysis based on the average normalized expression values of genes in the four luminal cell clusters. L1 had the most upregulated pathways relative to the other luminal clusters. g T-SNE plots of luminal cell clustering separated by treatment group. h Cell distributions from two independent experiments by percentage for all three treatment groups among all four luminal cell clusters. Symbols indicate individual cell distributions for each of the two experimental replicates (triangle, circle)

Fig. 4

PBDEs enhance E2 effects by acting on key regulators of mammary gland development. a Violin plot and feature plots of Esr1 expression in the luminal cell clusters. Vehicle-treated cells showed higher expression of Esr1 relative to E2 and E2 + PBDE-treated cells. b Violin plot and feature plots of Pgr expression in the luminal cell clusters. Pgr⁺ cells and Pgr expression increased in E2-present cells relative to vehicle-treated cells. c Feature plots showing co-expression of Esr1 and Pgr. The color code represents the binned and scaled expression data. Pgr expression is not specific to only Esr1⁺ cells. d Immunofluorescence analysis of ERα (red), PR (green), and DNA stain using DAPI (blue). Arrow indicates each cell types; ER+ only (red), PR+ only (green), and DP (yellow). Scale bars = 100 μm. e Quantification of Ki67+ ducts on whole mammary gland sections. IHC slides of mammary glands from vehicle (n = 4), E2 (n = 4), and E2 + PBDE (n = 6) were stained for Ki67. Transverse represents ducts that have a clear lumina, and end buds were defined as ducts that have no clear lumina and are densely packed. High Ki67 levels were defined as ducts that have most of its cells positively stained for Ki67, while low Ki67 levels were defined as ducts that are primarily Ki67−, but do have a few stained cells. f Representative images of ducts classified as high Ki67, low Ki67, and none for quantification. Scale bars = 50 µm. g Feature plot of Areg expression in luminal cell clusters. Areg expression was mostly specific to the intermediate and mature luminal clusters (L1, L2). h Bar plots of results from two independent experiments for luminal cells expressing Areg. PBDE treatment in addition to E2 further increased the number of Areg⁺ cells. Symbols indicate individual cell distributions for each of the two experimental replicates (triangle, circle)

Fig. 5

ECM associated and fibroblast cells were re-clustered to show a specific subpopulation of Ccl2⁺ cells. a Feature plot highlighting Epcam, Cd45, Cd31, Ng2, Cd52, and Krt14 negative cells to negatively select for ECM and fibroblast cells. b T-SNE plots of re-clustered ECM/Fibroblast clusters. c Bar plots showing the distribution of cells from two independent experiments in each ECM/Fibroblast cluster. Symbols indicate individual cell distributions for each of the two experimental replicates (triangle, circle). d Dendrogram showing the relationship of each ECM/fibroblast cluster with other clusters. e Heatmap showing the top 5 differentially expressed genes for each of the 7 ECM/Fibroblast clusters. f Heatmap on the MSigDB analysis based on the average normalized expression values of genes in the 7 ECM/Fibroblast cell clusters. Cluster F2 had the most upregulated pathways relative to all other fibroblast clusters. g Violin plots of selected genes emphasizing the unique gene signature of Esr1⁺ ECM/Fibroblast cells. Clusters F0, F1, and F2 had more cells with higher expression of Esr1 and Ccl2, whereas Egfr was observed in nearly all fibroblast clusters

Fig. 6

Re-clustered immune cell clusters C6, C8, and C9 resulted in 11 subclusters. a Feature plot highlighting C6, C8, and C9, which were putatively identified as immune cells based on high Ptprc expression in these clusters. b T-SNE plots of re-clustered immune cells. c Heatmap showing the top 5 differentially expressed genes for each of the 11 immune cell clusters. d Dendrogram highlighting the relationship among each of the immune cell sub-clusters. e Feature plots of selected genes to further classify immune cells as T cells (Cd3e, Gata3), NK cells (Klrb1c), dendritic cells (Cd24a), B cells (Cd19), macrophages (Fcgr1), and fibrocytes (Col1a1). f Feature plot of C1qa, which emerged as a differentially expressed gene between macrophage clusters I5 and I6. g Bar plots showing the distribution of cells from two independent experiments in each of the immune cell clusters. Symbols indicate individual cell distributions for each of the two experimental replicates (triangle, circle)

Fig. 7

Proposed model summarizing two different mechanisms of M2 macrophage recruitment and polarization by E2 or E2 + PBDE treatment. E2 enriches ECM/fibroblast cells expressing Ccl2, which in turn recruits M2 macrophages. PBDE, in the presence of E2, increases Areg expression in Esr1⁺Pgr⁺ luminal cells, which promotes M2 macrophage polarization through the E2–AREG–EGFR–M2 macrophage pathway