E2 treatment inhibited the proliferation of GS3 in vivo and in vitro
A consistent regression of GS3 was observed by the treatment of tumor-implanted mice with E2 (1 mg) for 28 days (Fig. 1a). This phenotype was preserved across serially transplanted GS3 tumors (up to 14 passages so far). GS3 tumor pieces that were simultaneously transplanted with E2, showed no growth for half a year; meanwhile GS3 tumor pieces implanted alongside a placebo pellet (no E2) grew beyond 1000 mm3 (Fig. 1b). A repeating pattern of regression and growth of the GS3 tumor volume was observed when subjected to intermittent E2 treatment. After three rounds of intermittent E2 treatment on GS3, an E2 independent growth developed (Fig. 1c). To examine 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. 1d).
E2 downregulated the expression of ERα and cell cycle proliferation genes in GS3
ERα and ERβ genes in GS3 were wild-type and not amplified. The only other reported estrogen-suppressive PDX, WHIM16, has an ESR1 amplification [17, 24]. 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. 2a), indicating the participation of ERα in this suppression process. H&E staining showed that cell density (the number of cells per square) decreased after E2 treatment (Fig. 2b). IHC indicated that the proportion of ERα+ cells decreased, and ERα staining intensity was reduced after E2 treatment depending on treatment terms (Fig. 2b). Western blotting and RPPA analysis showed that the E2 treatment resulted in lower ERα protein level in vivo (Fig. 2c, d). E2 treatment decreased the ESR1 expression of GS3-organoids at the mRNA level in vitro as well (Fig. 2e). IHC also showed that progesterone receptor (PR) + cells appeared in E2-treated tumors (Fig. 2b) and that E2 treatment increased the levels of PR protein and of mRNA (Fig. 2d, e, respectively). However, the number of Ki-67+ cells in IHC (Fig. 2b) and the CDK1 and cyclin-B1 expressions at the protein level in RPPA (Fig. 2d) decreased after E2 treatment.
To investigate the overall effects of E2 treatment on GS3, bulk RNAseq was performed using 5-, 10-, and 53-day E2-treated tumors as well as control tumors. GSEA analysis showed that the hallmark estrogen response early/late gene sets were upregulated after E2 treatment (Table S2). Although the ESR1 expression level decreased (possibly in part due to reduction of ERα+ cells, as indicated by IHC), expression levels of ER target genes, such as PGR, GREB1, and TFF1, increased depending on the E2 treatment time (Fig. 2f). In addition, cell cycle progression genes, such as CDK1, TOP2A, E2F2, and MKI67, were down-regulated in E2-treated tumors (Fig. 2f). The organoids isolated from longer E2-treated tumors expressed a higher level of ER target genes (Fig. 2g). In addition, the expression of two tumor suppressor genes, IL24 and GADD45A, was upregulated in E2-treated PDXs, as well as in organoids that were isolated from E2-treated tumors, depending on treatment time (Fig. 2f, g). 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 (as an indication of the presence of functional active ERα), and an induced cell cycle arrest and growth suppression.
scRNAseq to study effects of E2 on GS3
To better understand the mechanism of estrogen-induced growth suppression and the potential contributions of ERα+ cells and ERα− cells, we performed scRNAsEq. Following a treatment which lasted 7 days using E2 (1 mg) or a placebo pellet, GS3 tumors were harvested and digested into a single cell suspension. Although results of longer E2/placebo treatments showed larger differences of gene expressions according to the bulk RNAseq analysis, the viability of single cells isolated from GS3 tumors decreased after the tumors shrunk. Since the presence of a large number of dead cells would affect the quality of single cell analysis, we decided to treat GS3-PDX for 7 days in which the single cells’ viability was more than 80%. Single cell samples were prepared for scRNAseq using a 10x Genomics platform. Single cell preparations from two tumors from each treatment were combined and processed for scRNAseq (Fig. 3a). After filtering out mouse-derived cells (8% of total cells) based on their species-specific sequences, 20,467 genes were detected from 10,631 cells: 3,927 cells from placebo-treated tumors and 6,704 cells from E2-treated tumors. These were all epithelial cells that expressed human genes. The principle component analysis was performed on 2,000 HVGs. A UMAP was generated to summarize and visualize the data in a two-dimensional subspace, which led to the identification of 4 major cell clusters (C0–C3) in GS3 (Fig. 3b), where each dot represented one cell. The number of cells in each cluster is specified in Table S3.
Characteristics of 4 single cell clusters
The relationships between 4 clusters were visualized by a heat map using the top 20 differentially expressed genes (DEGs) and hierarchical clustering based on the average expression of HVGs (Fig. 3c, d). Top 20 DEGs in C0 and C1 included commonly known ER target genes. Those top 20 DEGs in C2 included genes of the hallmark TNFA signaling via NFKB gene set, and those in C3 included cell cycle progression genes. Although cells with a mitochondrial read rate > 20% were already excluded, C0 and C2 included cells with higher expression of mitochondrial genes (Fig. 3c). The gene expression pattern of C0 was relatively similar to that of C2, while C1 was similar to C3 (Fig. 3d). GSEA analysis of each cluster indicated that there were three important hallmark gene sets in GS3: the hallmark estrogen response early/late and G2M checkpoint (Table S4). The hallmark estrogen response early/late gene sets were upregulated in C0 and downregulated in C3. In contrast, the hallmark G2M checkpoint gene set was significantly downregulated in C0 and upregulated in C3 (Table S4).
Impact of E2 on gene expression at the single-cell level
Cells from E2- or placebo-treated tumors were clearly placed in different clusters (Fig. 3e). More than 90% of the cells in C0 and C1 were E2-treated cells and more than 95% of the cells in C2 and C3 were placebo-treated cells (Fig. S3a). The cell cycle scoring was performed using the Seurat package according to the developer’s vignette. C0/C2 included more G1 phase cells while C1/C3 included more G2M phase cells (Fig. 3f, Fig. S3b). GSEA analysis of the direct comparison between C0 and C2 showed that the hallmark estrogen response early/late gene sets were upregulated in C0, where most of the cells were from E2-treated tumors. Comparing C1 to C3 directly, the hallmark G2M checkpoint was downregulated in C1, where the majority of the cells were also from E2-treated tumors. The hallmark estrogen response early/late gene sets were upregulated in C1 compared to C3 as well (Table S4).
Comparison of E2-treated cells vs. placebo-treated cells
Most ESR1+ cells were distributed between C0 and C2 (Fig. 4a). The percentage of ESR1+ cells decreased by 10% in E2-treated cells compared to placebo-treated cells. However, the percentage of cells expressing ER target genes, such as PGR, PDZK1, SERPINA1, and GREB1, dramatically increased in E2-treated cells (Fig. 4a). Furthermore, the percentage of cells expressing cell cycle progression genes, such as MKI67, CCNE2, CCNA2, and E2F2, decreased in E2-treated cells (Fig. 4b). Distribution of cells remaining at the G1, S, or G2M phase in GS3 changed after 7 days of E2 treatment. The percentage of cells which progressed to the G2M phase decreased (-4.7%) and the percentage of cells that arrested at the G1 phase increased (+ 12.5%) in E2-treated cells (Fig. 4c).
Crosstalk between ESR1+ cells and ESR1− cells
The percentage of cells which were 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. The decreasing of G2M phase cells 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. Interestingly, the hallmark estrogen response early/late gene sets were upregulated not only in ESR1+ cells, but also in ESR1− cells after E2 treatment. Moreover, the hallmark G2M checkpoint gene set was significantly downregulated in ESR1− cells after E2 treatment (Table 1). Therefore, our single cell analysis revealed that E2 treatment upregulated the expression of estrogen response genes and reduced the proliferation of both ESR1+ and ESR1− cells in GS3.
Table 1
Hallmark gene sets analysis of ESR1+ cells and ESR1− cells
|
Upregulated in E2 (vs. placebo)
|
P-value
|
Downregulated in E2 (vs. placebo)
|
P-value
|
ESR1− cells
|
ESTROGEN_RESPONSE_LATE
|
1.22E-16
|
G2M_CHECKPOINT
|
1.45E-30
|
|
ESTROGEN_RESPONSE_EARLY
|
2.99E-15
|
E2F_TARGETS
|
2.83E-22
|
|
HTNFA_SIGNALING_VIA_NFKB
|
1.06E-07
|
TNFA_SIGNALING_VIA_NFKB
|
1.74E-12
|
ESR1+ cells
|
ESTROGEN_RESPONSE_EARLY
|
2.69E-19
|
TNFA_SIGNALING_VIA_NFKB
|
1.59E-20
|
|
ESTROGEN_RESPONSE_LATE
|
7.36E-18
|
P53_PATHWAY
|
7.78E-13
|
|
ANDROGEN_RESPONSE
|
5.45E-05
|
APOPTOSIS
|
1.07E-11
|
E2-induced IL24+ cells through ERα
The percentage of IL24+ cells in the E2-treated was 2.75-fold of that in placebo-treated cells (Fig. 5a). IL24+ cells included more ESR1+ cells compared to IL24− cells. In IL24+ cells, the percentage of cells that progressed to the G2M phase decreased (-12.4%) and those that arrested at the G1 phase increased (+ 9.7%) when compared to those in the IL24− cells (Fig. 5b). IL24+ cells expressed higher levels of mitochondrial genes than 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).
The GS3 tumors obtained E2 independence after the three rounds of intermittent E2 treatment. The decrease of ERα+ cells and Ki-67+ cells after E2 treatment was reversed after intermittent E2 treatment (Fig. 2b). To identify the key pathway driving 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 ER target genes), in which 40% of the genes behaved similarly to the placebo-treated sample (mainly cell cycle progression genes) (Fig. 5e). The IL24 expression level was upregulated in E2-treated tumors when compared to placebo-treated tumors and downregulated in intermittent E2-treated tumors when compared to (continuous) E2-treated tumors (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/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 2).
Table 2
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.06E-37
|
G2M_CHECKPOINT
|
1.05E-08
|
ESTROGEN_RESPONSE_LATE
|
7.60E-12
|
E2F_TARGETS
|
1.05E-08
|
ESTROGEN_RESPONSE_EARLY
|
4.46E-10
|
MYC_TARGETS_V1
|
6.05E-05
|
APOPTOSIS
|
6.03E-09
|
P53_PATHWAY
|
2.71E-03
|
Predication of the transitions between placebo-treated cells and E2-treated cells by trajectory analysis
To investigate the transition between placebo-treated cells and E2-treated cells, single cell trajectory analysis was performed using Monocle, an algorithm that explores relationships between cells through pseudo-time. There were 12 different ‘cell states’ identified during the trajectory analysis (Fig. 6a). The number of cells and distribution of cells treated with E2/placebo in each state are in Table S5 and Fig. S4a. The trajectory colored by four clusters in the UMAP (Fig. 3b) was shown in Fig. S4b. Trajectory analysis of single cells revealed three major branches associated with [I] placebo (State 1–2), [II] E2 (State 3–11), and [III] common E2 and placebo (State 12) (Fig. 6b). In the E2 branched group, ESR1+ cells and G1 phase cells progressively decreased from state 9 (the endpoint of the E2 branch) to state 3 (junction), and G2M phase cells increased inversely with ESR1+ cells and G1 phase cells (Fig. 6c, d). The hallmark estrogen response early/late gene sets were upregulated in both ESR1+ cells and ESR1‒ cells (Fig. 6e). Significantly, the hallmark G2M checkpoint gene set was downregulated, even in ESR1‒ cells (Fig. 6e). IL24+ cells and PGR+ cells were primarily distributed in state 9, which was the endpoint of the E2 branched group (Fig. 6f, Fig. S4c). There were very few IL24+ cells within state 9 at the G2M phase (Fig. 6g). As expected, in the placebo branched group, the hallmark estrogen response early/late gene sets were not upregulated, and the hallmark G2M checkpoint gene set was not downregulated in ESR1+ cells, when compared to the E2 branched group (Fig. 6e). In the common branch (containing cells from both placebo- and E2-treated tumors) that included both ESR1+ cells and ESR1‒ cells, there were no G1 phase cells and the hallmark G2M checkpoint was dramatically upregulated(Fig. 6e). MKI67+ cells (at G2M phase), including some ESR1+ cells, were mainly present in the placebo branch (State 1) and the common E2 and placebo branch (State 12) (Fig. S4c).
E2-induced apoptosis after the long-term E2 treatment
To investigate the relationship between tumor shrinkage and apoptosis, we detected apoptotic cells from formalin-fixed tumor tissues. A few apoptotic cells were observed in 28-day E2-treated tissue (Fig. 7a). There was a significant difference between placebo and 28-day E2-treated tissues (P < 0.0001, Fig. 7b). Western blotting analysis further demonstrated that long term (28-day) E2 treatment induced expression of apoptosis-associated proteins, as indicated by elevated levels of cleaved PARP (a hall mark of apoptosis) and the pro-apoptotic protein Bax, as well as by reduced levels of the pro-survival protein Bcl-xl. Cleaved-PARP expression was reduced after GS3 tumors gained E2 independence. (Fig. 7c).