SC31 and GS3 behaved oppositely regarding tumor growth with estrogen
Measurements of tumor volume showed that E2 promoted growth of SC31 and suppressed growth of GS3 (Fig. 1A). A consistent promotion of SC31 and regression of GS3 were observed by the E2 treatments in which an E2 pellet (1 mg) was implanted in mice with SC31/GS3 tumor initially (Fig. 1A, P = 0.016, P = 0.008, respectively), or after tumors grew up to 150–200mm3 (Fig. 1A, P = 0.04, P = 0.0006, respectively). These growth phenotypes were maintained after serially transplanted SC31/GS3 tumors (up to 13/14 passages so far). H&E staining showed that cell density (the number of cells per microscopic field) increased in SC31 and decreased in GS3 after E2 treatment (Fig. 1B). IHC indicated that the proportion of ERα+ cells, ERα staining intensity, and the number of Ki-67+ cells were increased in SC31 and decreased in GS3 after E2 treatment (Fig. 1B). However, IHC showed that PR+ cells appeared in both SC31 and GS3 after E2 treatment (Fig. 1B).
E2 downregulated the expression of ERα and cell cycle proliferation genes in GS3
Because the growth response to E2 in these two ER+ breast cancers were unexpectedly different, we performed molecular analyses to determine the mechanism of the different response in GS3, the tumor suppressed by E2. These analyses included exome sequencing, bulk RNA sequencing and in vitro validation. Human whole exome sequencing of GS3 did not detect ESR1 and ESR2 variants (Table S2). ERα and ERβ genes in GS3 were wild-type and not amplified. The only other reported estrogen-suppressing PDX model, WHIM16, has an ESR1 amplification [14, 15]. To verify the effect of E2 on the proliferation of GS3 in vitro, organoids were generated from untreated GS3 tumors. E2 suppressed the proliferation of GS3 organoids in an E2 concentration-dependent manner (Fig. 2A). These findings demonstrated that estrogen/ERα signaling could modulate ER+ tumor growth by unusual mechanisms. To determine which ER subtype is involved in E2-induced regression of GS3, we performed co-treatment of E2 and ERα-specific antagonist (MPP) or ERβ-specific antagonist (PHTPP) in vitro. E2-mediated inhibition of GS3-organoids could be reversed by the co-treatment of MPP, but not by PHTPP (Fig. 2B), indicating the participation of ERα in this suppression process. Western blotting and RPPA analysis showed that the E2 treatment resulted in lower ERα protein level in vivo (Fig. 2C, 2D). E2 treatment decreased the ESR1 expression of GS3-organoids at the mRNA level in vitro as well (Fig. 2E). E2 treatment increased the levels of PR protein and mRNA (Fig. 2D, 2E, respectively).
To investigate the overall effects of E2 treatment on GS3, bulk RNAseq was performed on tumors from mice treated with E2 for 5, 10 and 53 days. GSEA analysis showed that the hallmark estrogen response early and late gene sets were upregulated after E2 treatment (Table S3). Although the ESR1 expression level decreased (possibly in part due to reduction of ERα+ cells, as indicated by IHC), expression levels of estrogen-regulated genes, such as PGR, GREB1, and TFF1, increased after E2 treatment (Fig. 2F). In addition, cell cycle progression genes, such as CDK1, TOP2A, E2F2, and MKI67, were down-regulated in E2-treated GS3 tumors (Fig. 2F). The organoids isolated from longer E2-treated tumors expressed estrogen-regulated genes at a high level (Fig. 2G). In addition, the expression of two tumor suppressor genes, IL24 and GADD45A, were upregulated in E2-treated PDXs, as well as in organoids that were isolated from E2-treated tumors, depending on treatment time (Fig. 2F, 2G). Collectively, our examination of GS3 tumors after E2 treatment indicates that inhibition of GS3 growth is associated with a reduction of ERα+ cells, an increased expression of ERα target genes (indicating the presence of functional active ERα), and an induced cell cycle arrest and growth suppression.
Impact of E2 on gene expression at the single-cell level in SC31 and GS3
To better understand the mechanism of estrogen-regulatory effects and the potential contributions of ERα+ cells and ERα– cells within the tumors, we performed scRNAseq on individual cells from SC31 and GS3. After removing mouse cells from the single cell suspension based on species-specific DNA markers (8% of total cells), 19,568 genes were detected from 5,581 cells in SC31 (3,110 cells from placebo-treated tumors and 2,471 cells from E2-treated tumors), and 20,467 genes were detected from 10,631 cells in GS3 (3,927 cells from placebo-treated tumors and 6,704 cells from E2-treated tumors). A UMAP plot led to the identification of 8 major cell clusters in SC31 (Fig. 3A) and GS3 (Fig. 3B), where each dot represented one cell. The number of cells in each cluster is specified in Table S4. In SC31, C3 and C6 included more placebo-treated and G1 phase cells while C2 and C5 included more E2-treated and G2/M phase cells (Fig. 3A). In GS3, C1 and C7 included more E2-treated and G1 phase cells while C3, C5 and C6 included more placebo-treated and G2/M phase cells (Fig. 3B).
Characteristics of single-cell clusters in SC31 and GS3
The relationships between 8 clusters of each single-cell analysis were visualized by a heat map using the top 10 differentially expressed genes (DEGs) and hierarchical clustering based on the average expression of HVGs. In SC31, top 10 DEGs in C2, C4, and C5 included genes belonging to the hallmark G2M checkpoint gene set (Fig. 3C, Table S5). A dendrogram showed that C2 and C5, where approximately 62 % of the cells were from E2-treated tumors, were most similar (Fig. 3D), and cells in those clusters were mainly in G2/M or S phase (Fig. S4A). GSEA analysis showed that the hallmark mTORC1 signaling gene set was upregulated (P = 3.30E-41) in C6 (Table S5). In GS3, top 10 DEGs in C0 and C2 included estrogen response genes (Fig. 3E), and cells in those clusters were mainly G1 phase cells (Fig. S4B). The hallmark G2M checkpoint gene set was downregulated in these clusters (P = 2.02E-22 and P = 3.76E-25, respectively) according to GSEA analysis (Table S6). In C1, which included approximately 90% of G1 phase cells (Fig. S4B), the hallmark estrogen response early and late gene sets were upregulated (P = 5.13E-06 and P = 3.90E-07, respectively), and the hallmark G2M checkpoint gene set was significantly downregulated (P = 5.62E-25) at the same time (Table S6). C1 and C0, where the hallmark G2M checkpoint gene set was downregulated, were relatively similar according to a dendrogram (Fig. 3F).
Comparison of E2-treated cells vs. placebo-treated cells in SC31 and GS3
In SC31, GSEA analysis showed that the hallmark G2M checkpoint gene set and the hallmark estrogen response early and late gene sets were significantly upregulated (P = 8.99E-72, 5.43E-52, 5.43E-52, respectively) after E2 treatment (Table 1). Cells expressing ESR1 or estrogen-regulated gene PGR were distributed mainly in C1, C3, and C6 (Fig. S4C). The percentage of ESR1+ cells and PGR+ cells increased by 31% and 37% respectively in E2-treated cells compared to placebo-treated cells (Fig. 4A). A majority of MKI67+ cells were found in C3 and C6 (Fig. S4C), and the percentage of MKI67+ cells increased by 15% in E2-treated cells (Fig. 4A). In GS3, GSEA analysis showed that the hallmark estrogen response early and late gene sets were upregulated (P = 1.63E-17 and 4.55E-16, respectively), while the hallmark G2M checkpoint gene set was downregulated (P = 3.95E-13) after E2 treatment (Table 1). Most ESR1+ cells were present in C1, C2 and C6, and the percentage of ESR1+ cells decreased by 10% in E2-treated cells compared to placebo-treated cells (Fig. 4B, S4D). Cells expressing PGR+ were found in C1, C2 and C6, in the same distribution as ESR1+ cells, and the percentage of cells expressing PGR increased by 11% in E2-treated cells (Fig. 4B, S4D). Furthermore, most of cells expressing MKI67 were primarily in C6. The percentage of MKI67+ cells decreased by 4% in E2-treated cells (Fig. 4B, S4D). The percentage of cells expressing other estrogen-regulated genes, such as AREG and PDZK1, also increased in E2-treated cells in both SC31 and GS3 (Fig. S5A, S5B). The percentage of cells expressing other cell cycle progression genes, such as CCNE2 and CDK1, increased in E2-treated cells in SC31, while decreased in GS3 (Fig. S5A, S5B).
Distribution of cells in the G1, S, or G2/M phases changed after E2 treatment. In SC31, the percentage of cells which progressed to the G2/M phase increased (+6.6%) and the percentage of cells arrested at the G1 phase decreased (-14.5%) in E2-treated cells (Fig. S6A). In GS3, the percentage of cells which progressed to the G2/M phase decreased (-4.7%) and the percentage of cells arrested at the G1 phase increased (+12.5%) in E2-treated cells (Fig. S6B). Collectively, the estrogen-regulated gene, PGR, was expressed after E2 treatment in both SC31 and GS3, while the percentage of cells expressing ESR1 and MKI67 or entering G/2M phase increased in SC31 but decreased in GS3.
Table 1. Hallmark gene sets analysis of SC31 and GS3 (E2 treated vs. Placebo)
|
Upregulated in E2
|
P-value
|
Downregulated in E2
|
P-value
|
SC31
|
E2F_TARGETS
|
5.94E-84
|
INTERFERON_GAMMA_RESPONSE
|
5.43E-64
|
|
MYC_TARGETS_V1
|
2.88E-76
|
INTERFERON_ALPHA_RESPONSE
|
6.33E-64
|
|
G2M_CHECKPOINT
|
8.99E-72
|
TNFA_SIGNALING_VIA_NFKB
|
1.70E-33
|
|
ESTROGEN_RESPONSE_EARLY
|
5.43E-52
|
APOPTOSIS
|
6.64E-25
|
|
ESTROGEN_RESPONSE_LATE
|
5.43E-52
|
P53_PATHWAY
|
2.15E-18
|
GS3
|
ESTROGEN_RESPONSE_LATE
|
1.63E-17
|
TNFA_SIGNALING_VIA_NFKB
|
6.06E-21
|
|
ESTROGEN_RESPONSE_EARLY
|
4.55E-16
|
G2M_CHECKPOINT
|
3.95E-13
|
|
|
|
P53_PATHWAY
|
3.95E-13
|
|
|
|
APOPTOSIS
|
5.90E-12
|
|
|
|
E2F_TARGETS
|
6.23E-12
|
Footnote: This table included top five hallmark gene sets or those which P value was smaller than 1.0E-10.
Effects of E2 treatment on ESR1+ cells and ESR1– cells in SC31 and GS3
Since the percentage of ESR1+ cells changed oppositely in SC31 and GS3, we decided to compare ESR1+ cells and ESR1– cells by analyzing them separately. The percentage of ESR1+ and ESR1– cells in two tumors treated with E2 or placebo was shown in Table S7. In SC31, the percentage of cells progressing to the G2/M phase increased by 7.2% in ESR1+ cells and by 5.8% in ESR1– cells after E2 treatment. Furthermore, the percentage of cells in the G1 phase decreased by 15.9% and 10.9% in ESR1+ and ESR1– cells respectively (Fig. 4C). In GS3, the percentage of cells arrested at the G1 phase increased by 10.9% in ESR1+ cells and by 16.0% in ESR1– cells after E2 treatment. In addition, the percentage of cells progressing to the G2M phase decreased by 3.4% and 7.8% in ESR1+ and ESR1– cells, respectively. ESR1– cells included more G2M phase cells compared to ESR1+ cells in GS3 (Fig. 4D). However, the decreasing percentage of G2/M phase cells in GS3 after E2 treatment was greater in ESR1– cells compared to ESR1+ cells (Fig. 4D).
To clarify the difference of signaling pathways between ESR1+ cells and ESR1– cells, we performed GSEA analysis for ESR1+ cells and ESR1– cells separately. Firstly, we focused on placebo-treated cells to learn original characteristics of two tumors. In SC31, the hallmark G2M checkpoint gene set was upregulated in ESR1+ cells (P = 1.99E-08), and the hallmark mTORC1 signaling gene set was upregulated in ESR1– cells compared to ESR1+ cells (P = 1.56E-03, Table S8). In GS3, the hallmark TNFA/NFκB signaling gene set was upregulated in ESR1+ cells (P = 4.50E-34) and the hallmark G2M checkpoint gene set was upregulated in ESR1– cells (P = 7.91E-15, Table S8). Secondly, we focused on the comparison of placebo-treated vs. E2-treated cells. In SC31, the hallmark G2M checkpoint gene set and estrogen response early and late gene sets were upregulated not only in ESR1+ cells (P = 1.98E-60, 1.35E-06, 8.26E-08, respectively), but also in ESR1– cells (P = 1.64E-25, 1.13E-07, 4.72E-09, respectively) after E2 treatment (Table 2). In GS3, the hallmark estrogen response early and late gene sets were upregulated not only in ESR1+ cells (P = 2.69E-19, P = 7.36E-18, respectively), but also in ESR1– cells (P = 1.22E-16, P = 2.99E-15, respectively) after E2 treatment. Moreover, the hallmark G2M checkpoint gene set was significantly downregulated in ESR1– cells (P = 1.45E-30) after E2 treatment (Table 2). Therefore, our single-cell analysis revealed that E2 treatment upregulated the expression of estrogen response genes in all epithelial cells of both SC31 and GS3. While estrogen upregulated the expression of cell cycle proliferation genes of SC31 in both ESR1+ and ESR1– cells, it downregulated the expression of cell cycle proliferation genes of GS3.
Table 2. Hallmark gene sets analysis of ESR1+ cells and ESR1– cells (E2 treated vs. Placebo)
|
|
Upregulated in E2
|
P-value
|
|
Downregulated in E2
|
P-value
|
SC31
|
ESR1+ cells
|
E2F_TARGETS
|
2.91E-83
|
|
INTERFERON_GAMMA_RESPONSE
|
6.11E-64
|
|
|
G2M_CHECKPOINT
|
1.23E-72
|
|
INTERFERON_ALPHA_RESPONSE
|
6.81E-63
|
|
|
MYC_TARGETS_V1
|
1.10E-69
|
|
TNFA_SIGNALING_VIA_NFKB
|
1.06E-34
|
|
|
ESTROGEN_RESPONSE_EARLY
|
6.52E-57
|
|
APOPTOSIS
|
2.82E-23
|
|
|
ESTROGEN_RESPONSE_LATE
|
1.67E-51
|
|
P53_PATHWAY
|
3.16E-21
|
|
ESR1– cells
|
E2F_TARGETS
|
1.94E-62
|
|
INTERFERON_ALPHA_RESPONSE
|
4.80E-66
|
|
|
MYC_TARGETS_V1
|
1.94E-62
|
|
INTERFERON_GAMMA_RESPONSE
|
7.63E-66
|
|
|
G2M_CHECKPOINT
|
1.44E-48
|
|
TNFA_SIGNALING_VIA_NFKB
|
7.61E-34
|
|
|
ESTROGEN_RESPONSE_EARLY
|
2.49E-43
|
|
APOPTOSIS
|
3.63E-25
|
|
|
ESTROGEN_RESPONSE_LATE
|
4.75E-42
|
|
HYPOXIA
|
3.83E-16
|
GS3
|
ESR1+ cells
|
ESTROGEN_RESPONSE_EARLY
|
1.11E-19
|
|
TNFA_SIGNALING_VIA_NFKB
|
2.20E-22
|
|
|
ESTROGEN_RESPONSE_LATE
|
3.20E-18
|
|
P53_PATHWAY
|
3.69E-13
|
|
|
|
|
|
APOPTOSIS
|
5.57E-12
|
|
ESR1– cells
|
ESTROGEN_RESPONSE_LATE
|
5.61E-17
|
|
G2M_CHECKPOINT
|
3.05E-31
|
|
|
ESTROGEN_RESPONSE_EARLY
|
1.44E-15
|
|
E2F_TARGETS
|
8.25E-23
|
|
|
|
|
|
TNFA_SIGNALING_VIA_NFKB
|
5.48E-14
|
Footnote: This table included top five hallmark gene sets or those which P value was smaller than 1.0E-10.
E2-induced IL24+ cells through ERα only in GS3
Interleukin-24 (IL24) has been shown to be a tumor suppressor [16–20]. In SC31cells, no IL24 expression was observed in either E2-treated or placebo conditions. (Fig. 5A). In GS3 cells, the percentage of IL24+ cells in the E2-treated category was 2.75-fold of that in placebo-treated cells (Fig. 5A). More IL24+ cells expressed ESR1+ compared to IL24– cells (Fig. 5B). In IL24+ cells, the percentage of cells that progressed to the G2/M phase decreased (-12.4%) and those that arrested at the G1 phase increased (+9.7%) when compared to those in the IL24– cells (Fig. 5C). When comparing IL24 expression levels in total mRNA extracted from GS3 tumors, IL24 expression level was upregulated after E2 treatment in a time dependent manner, and its expression was inhibited by ICI (Fig. 5D). To clarify the difference of signaling pathways between IL24– cells and IL24+ cells, we performed GSEA analysis for IL24– cells and IL24+ cells separately. The hallmark TNFA signaling via NFKB, hallmark estrogen response early and late, and hallmark apoptosis gene sets were upregulated in IL24+ cells after E2 treatment. On the other hand, the hallmark G2M checkpoint gene set was downregulated in IL24+ cells after E2 treatment (Table 3).
Table 3. GSEA hallmark gene sets analysis of genes expressed on IL24+ cells vs. IL24– cells
Upregulated in IL24+ cells
|
P-value
|
|
Downregulated in IL24+ cells
|
P-value
|
TNFA_SIGNALING_VIA_NFKB
|
1.33E-17
|
|
G2M_CHECKPOINT
|
8.25E-09
|
ESTROGEN_RESPONSE_ EARLY
|
1.73E-09
|
|
E2F_TARGETS
|
8.25E-09
|
ESTROGEN_RESPONSE_ LATE
|
1.73E-09
|
|
MYC_TARGETS_V1
|
5.24E-05
|
APOPTOSIS
|
2.45E-06
|
|
P53_PATHWAY
|
2.46E-03
|
Footnote: This table included all hallmark gene sets from the results of GSEA.
Intermittent E2 treatment in GS3
To investigate possible changes in response to E2 in GS3 tumors, we exposed PDX mice to cycles of intermittent E2 treatment (on and off every 28 days). A repeating pattern of tumor growth inhibition (E2 exposed) and growth (E2 absent) was initially observed. However, after three rounds of intermittent E2 treatment on GS3, an E2 independent growth developed (Fig. 6A). The decrease of ERα+ cells and Ki-67+ cells evaluated by after E2 treatment was reversed after intermittent E2 treatment (Fig. 6B), ERα protein level by Western blotting, as well (Fig. 2C). The increase of PR+ cells by IHC after E2 treatment was not reversed after intermittent E2 treatment (Fig. 6B). To identify the key pathway driving this emerging E2 independence, we performed bulk RNAseq on GS3 tumors treated with placebo, E2 (for 28 days), or intermittent E2. The hierarchical clustering using 1,312 DEGs showed that approximately 60% of genes in the intermittent E2-treated sample had the same trend as E2-treated samples (mainly estrogen-regulated genes), in which 40% of the genes behaved similarly to the placebo-treated sample (mainly cell cycle progression genes) (Fig. 6C). After E2 treatment for 28 days, the hallmark estrogen response late gene set was upregulated and the hallmark G2M checkpoint gene set was downregulated: the same results as scRNAseq analysis (Fig. S7A, S7B). Interestingly, the hallmark Reactive Oxygen Species (ROS) pathway gene set was upregulated after the long-term E2 treatment, despite the effective pathway size being relatively small (Fig. S7A). Upregulation of IL24 expression level and the hallmark ROS pathway in 28-day E2 was reversed in intermittent E2 samples (Fig. 5D, S7B). Downregulation of the hallmark G2M checkpoint gene set in 28-day E2 was also reversed in intermittent E2 samples (Fig. S7B).