Time-series transcriptome profiles of MCF-7 cells during continuous TAM treatment
We first investigated the effect of continuous TAM treatment on human breast cancer MCF-7 cells, whose growth depends on ER signaling (Fig. 1a). Treatment with 1 µM TAM initially inhibited cell growth (Fig. 1b) through decreasing the number of cells in S phase whereas increasing that in G1 phase (Fig. 1c and 1d), suggesting that TAM treatment induced G1 arrest. The growth of cells was almost completely inhibited until week 5 (W5) but recovered thereafter (Fig. 1b). The cell cycle of TAM-treated cells was also dysregulated until W5 but became comparable to control cells after W6 (Fig. 1c and 1d). These results showed the process by which the cells survived and restored their growth potential in the absence of ER signaling.
We next analyzed the bulk RNA-seq data of TAM-treated MCF-7 cells to identify changes in the gene regulatory network and compare them with those in cellular phenotype. We previously performed time-course bulk RNA-seq analysis of TAM-treated MCF-7 cells and identified gene sets that play critical role at the “tipping” point of resistance acquisition15. In this study, we re-analyzed this dataset in order to focus on time-dependent changes in the expression of each gene, correcting for the effect of the difference in library preparation processes between samples (see Materials and methods and Supplementary Fig. 1). A total of approximately 6,000 differentially expressed genes (DEGs) were identified between TAM-treated and control cells, of which approximately 3,000 were up-regulated (Supplementary Fig. 2).
Gene expression in TAM-treated cells was normalized to control cells at each time point, and log2 fold-change (log2FC) values were calculated. The log2FC values of all genes at week 0 were set as a theoretical value of zero. We then obtained time-course patterns of log2FC values of 6,982 DEGs at 13 time points and analyzed them by principal component analysis (PCA). The PC1 axis classified the expression profiles of all DEGs into two groups; from W0 to W4, and from W5 to W12 (Fig. 1e). The separation in log2FC values was the greatest between weeks 4 and 5, which preceded the recovery of cell growth (Supplementary Fig. 3). In addition, the separation in log2FC values between W1 and W5 was greater than that between W6 and W12, indicating that gene expression became more stable at the later stages.
We then investigated the relationships between dynamic gene expression patterns and gene functions. Cluster analysis of the z-score of log2FC values revealed six groups of genes with distinct expression patterns: cluster A, rapidly decreasing expression; cluster B, initially up-regulated and recovered at week 5; cluster C, rapidly increasing before cell growth rate recovery; cluster D, a gradual increase in expression concomitant with growth rate recovery; cluster E, initially up-regulated and then down-regulated; and cluster F, gradually decreasing (Fig. 1f and Supplementary Table 1). The enrichment analysis of each group was carried out using the Reactome pathway (Fig. 1g) and KEGG (Fig. 1h) databases. Cluster A (rapidly decreasing expression pattern) was enriched in genes related to both receptor tyrosine kinase signaling and fatty acid metabolism. Genes in cluster C (rapid up-regulation) were related to TGFβ signaling or extracellular matrix–receptor interaction, such as TGFB2, SMAD3, and MMP9, or encoded collagen or laminin proteins. This implies that the reorganization of the gene network regulating epithelial–mesenchymal transition or extracellular matrix secretion, both of which contribute to cancer malignancy, occurred before cell growth recovery in the TAM-treated condition. Genes related to ribosomes were also enriched in cluster C, indicating that ribosomal biogenesis may be up-regulated during continuous TAM treatment. On the other hand, cluster D, in which gene expression level increased gradually, was enriched in genes functioning in the thyroid hormone signaling pathway, HIF1 pathway, and glycolytic process, and downstream signaling of RAS, among others. These data show that a Warburg-like effect co-occurs with TAM resistance acquisition. Genes in cluster B showed a non-monotinic dynamic expression pattern characterized by a transient decline from week 1 to week 4, followed by recovery to the basal level. This trend was similar to the growth rate pattern observed in the TAM-treated condition. Cluster B contained numerous genes involved in cell cycle regulation such as CCND1 and E2F1, and DNA replication such as RAD51, DNA damage (DNA repair, transcriptional regulation of TP53, and base excision repair), RNA metabolism, kinesins, beta-catenin degradation, generic transcription pathways, and the Fanconi anemia pathway. Genes in cluster E were up-regulated only when cell growth was effectively inhibited by TAM; this pattern was opposite to that observed in cluster B. Cluster E was enriched in genes involved in interferon signaling, FoxO signaling, autophagy, as well as transcription pathways, exytocin signaling pathways, and the MAPK signaling pathway. The interferon and FoxO signaling pathways exhibit anti-survival functions in cancer cells exposed to anti-cancer agents17, whereas autophagy contributes to cell survival under normal conditions18, suggesting that genes in cluster E reflect both antitumor as well as adaptation mechanisms triggered by the TAM treatment. Cluster F was the largest dynamic cluster containing 2,088 genes and was characterized by a consistent decrease in gene expression. Enrichment analysis showed that cluster F contained genes related to multiple cellular functions including energy metabolism, growth hormone synthesis, Rho GTPase (Rap1) signaling, and crucially, the estrogen signaling pathway, suggesting that TAM treatment inhibits the estrogen-dependent gene expression mechanism, and TAM resistance observed in our experiment may be supported by an estrogen/ER-independent mechanism.
Single-cell RNA-seq analysis of MCF-7 cells under continuous TAM treatment
Because cell-to-cell heterogeneity of phenotypic features is a key mechanism of drug resistance19, we investigated TAM-induced changes in gene expression profiles at a single-cell level. On the basis of the results of the cell growth assay and bulk RNA-seq data, we focused on four time points: week 0 (W0, just before starting TAM treatment), week 3 (W3, beginning of the period of complete cell growth inhibition), week 6 (W6, end of the period of complete cell growth inhibition), and week 9 (W9, at the acquisition of TAM resistance) (Fig. 2a). RNA-seq analysis of 1,108 single cells yielded 577 high-quality single-cell gene expressions (Fig. 2a; see the Materials and Methods section). Averaged single-cell expression profiles were well correlated with bulk expression (Supplementary Fig. 4). The Pearson correlation coefficient of 11,413 gene expression values between individual cells decreased at W3 and then gradually increased at W6 and W9 (Fig. 2b). This changing pattern of correlation coefficient might suggest the selection of cells that can survive the TAM treatment and subsequently transition into multiple stable states.
To visualize cell-to-cell diversity in detail, we conducted uniform manifold approximation and projection (UMAP), one of the standard methods of dimensional reduction. Before drawing the UMAP plot, we calculated the probability score of cell cycle progression in each cell using Seurat 3 software20 to correct for the bias caused by the difference in the cell cycle stage (Fig. 2c). All cells were mapped on a three-dimensional UMAP plot and projected in two dimensions (Figs. 2d, 2e, 2f, and Supplementary Fig. 5). Single-cell data were roughly divided into two groups: W0 and the others (W3, W6, and W9). Cells were widely distributed in space at W3 and W6 but were localized in two separate regions at W9. These cells could be clustered into six subpopulations in the UMAP plot (Fig. 2f). Cells in subgroups 1 and 6 belonged to the W0 group, and these cells were strongly diminished by W3. By contrast, cells in subgroups 2 and 3 newly emerged at W3 and prevailed in those clusters until W6. Finally, at W9 these cells splitted into two groups: one containing subgroup 4, and the other containing subgroup 5.
Subgroup-specific gene modules and their functions
We first investigated marker genes in each subgroup. The top five genes showing the highest specificity scores were selected in each subgroup (Supplementary Fig. 6a). Subgroups 1 and 6 were the major groups at W0, and marker genes in these subgroups included typical ER pathway target genes such as AREG21 and GREB122. This result showed that transcriptional activity of ER for typical target genes was down-regulated in the other subgroups. Comparing to the marker genes in subgroups 1 and 6, the expression level of marker genes in other subgroups, especially that in subgroup 2, does not clearly distinguish the cells into the subgroups. Interestingly, almost all marker genes of subgroups 4 and 5 also showed high expression in subgroup 3, suggesting that the pre-resistant subgroup 3 could potentially mature into distinct resistant subgroups by rewiring the genetic network.
Next, we analyzed the genetic modules specifically expressed in each subgroup or each week (Fig. 2g, Supplementary Fig. 6b, and Supplementary Table 2). Subgroups 1 and 6 contained highly expressed gene modules j, k, and l, which are enriched in ESR-mediated signaling, unfolded protein response, and amino acid and nucleotide metabolism. On the other hand, gene expression of module a was particularly low for these subgroups. Subgroups 2 and 3 were the major subpopulations in W3 and W6. Both these subgroups showed high expression levels of genes in module a, some of which are involved in interferon signaling, TGFβ signaling, and tight junctions. These enriched terms showed a strong resemblance to the early responsive cluster C in the bulk RNA-seq experiment (Figs. 1g and 1h). Subgroup 4, whose population was increased at W9, showed high expression levels of genes in modules g, h, and i, as shown in the heatmap. These genes encoded cell adhesion-related molecules such as integrin β4 (ITGB4), laminin β2 (LAMB2), and zyxin (ZYX), and some genes were involved in ROCK activation mechanisms. These modules also include several terms related to signal transduction, such as the VEGF pathway and thyroid hormone signaling. In addition, some chromatin remodeling enzymes and lysine-specific histone demethylases were included in module h. These results indicate that TAM-resistant cells in subgroup 4 showed higher activities of cell adhesion and migration, with an altered signaling pathway and epigenetic status. Subgroup 5, which represented another major population during W9, contained highly expressed genes in modules b, c, d, and e. This result indicates that genes related to innate immune responses, oxidative phosphorylation, and translation are highly expressed in the cells in subgroup 5. In addition, module c contained genes related to carbon metabolism, especially the glycolysis/glycogenesis pathway, suggesting that cells in subgroup 5 exhibit unique metabolic adaptation to TAM-induced stress. Based on the aforementioned results, we found that TAM-resistant ER-positive breast cancer cells obtained from the same parental cell line could be divided into two types, one of which acquired the re-wired metabolic network (subgroup 5) and another acquired high expression levels of adhesion molecules with changing of epigenetic status (subgroup 4).
Trajectory analysis of TAM resistance
To confirm the cell transition trajectory into two different types of resistant subgroups, we conducted pseudotime analysis (Figs. 3a, 3b, and 3c). The pseudotime of each cell calculated from the gene expression data was correlated with the sampling time after starting the continuous TAM treatment (Fig. 3d). The pseudotime of cells in subgroup 4 was higher than that of cells in subgroup 5, suggesting that cells in subgroup 4, showing high expression of epigenetic modulators, are more divergent from parental cells than cells in subgroup 5 (Fig. 3e). To indetify important molecules involved in the emergence of subgroups 4 and 5, we analyzed DEGs along the estimated cell trajectory. A total of 273 and 79 genes were detected as highly expressed genes in subgroups 4 and 5, respectively (Fig. 3f, Supplementary Tables 3 and 4). Then, we investigated the transcriptional regulators of the highly expressed genes in subgroup 4 and 5 using the ChIP-Atlas database23 (Fig. 3g and 3h), which covers almost all public ChIP-seq data. Prediction of proteins that bind near the transcriptional start site of DEGs in the trajectory to subgroup 4 showed that only 4 of the top 10 factors represented ChIP-seq data from MCF-7 cells, and most of the others were obtained from the triple-negative breast cancer (TNBC) cell line (Fig. 3g, left, shown in green). These results also suggest that most up-regulated genes in subgroup 4 are controlled by bromodomain-containing proteins, BRD4 and BRD2, which recognize acetylated histones and act as super enhancers24,25. In addition, our results also suggest the possible involvement of oncogenic transcription factors SMAD3 and ERG in the trajectory to subgroup 4. Genes encoding these molecules were up-regulated before the cells acquired TAM resistance in bulk RNA-seq data (Fig. 3h, top and middle), and related terms to these molecules (“Signaling by TGF-beta receptor complex” and “MAPK singaling pathway”) were detected in cluster C and E in bulk RNA-seq data, respectively (Fig. 1g and 1h). This indicates that subgroup 4 genes have different epigenetic status, which is clearly distinct from that of parent MCF-7 cells; this result was consistent with the enrichment analysis of specific gene modules (Fig. 2g).
In subgroup 5, 7 of the top 10 factors were obtained from ER-positive breast cancer or normal cells (Fig. 3g, right), suggesting that subgroup 5 retained the transcriptional network of parental MCF-7 cells. ChIP-seq data from anti-WDR5 antibody showed the best q-value and fold enrichment score. Although WDR5 was not up-regulated in bulk RNA-seq data (Fig. 3h, bottom), various binding molecules of WDR5, methylated histone H3 lysine 4 or MYC for example26, might regulate genes related to subgroup 5. Especially, MYC was detected as a candidate estimated from ChIP-Atlas database for upstream regulator of increased genes in subgroup 5. Other candidates included TAF1, PML, and BRD4. Among these genes, MYC, TAF1, and PML were up-regulated before week 4, as shown by bulk RNA-seq analysis. These data indicate that MYC, TAF1, and PML may contribute to one of the emerging TAM-resistant subpopulations. Taken together, our analysis revealed key molecular candidates that drive two different TAM-resistant subgroups.
Mathematical modeling of the TAM resistance acquisition process
We next constructed aphenomenological mathematical model that reproduces the population-level dynamics of TAM resistance, based on cell trajectories obtained by pseudotime analysis, to estimate the relative contribution of TAM-mediated cell growth- and differentiation to the acquisition of resistance (Fig. 4a). This model comprises cell transformation among four major cell subpopulations as estimated by single-cell gene expression profiles: cells initially sensitive to TAM (XS, subgroups 1 and 6 in Fig. 3a), pre-resistant cells (XP, subgroups 2 and 3), resistant cells showing altered expression of metabolism-related genes (XR1, subgroup 5), and resistant cells with highly adhesive phenotype with epigenetic changing (XR2, subgroup 4). The state transitions between the four subpopulations were assumed to follow the graph structure estimated by pseudotime analysis (Fig. 3c): XS cells change to XP in the presence of tamoxifen, and cells in XP change to the two resistant states XR1 and XR2. The model also enables cell transitions in the reverse direction. In addition, we assumed that continuous TAM treatment accelerates the rate of forward cell transition rate in response to a cumulative history of TAM treatment and describe this acceleration as a sigmoid function of the integral of TAM. This assumption is based on the previous findings that cell state transitions require the accumulation of genetic or epigenetic changes, which from the tipping point of resistance15. We explicitly considered extrinsic noise-mediated cell-to-cell variability by fitting 20 independent model parameter sets to two experimental time-course datasets describing cell proliferation and differentiation dynamics using the BioMASS computational framework27. Specifically, our model reproduced the experimentally observed dynamic distributions of both total cell growth rate in the presence of TAM (Fig. 1b) and of the proportions of the four different cell subpopulation (Fig. 2g, heatmap in green),, which could reproduce, simultaneously (Figs. 4b, 4c, and Supplementary Fig. 7a).
We found two remarkable features of the well-fitting parameter distributions. First, the growth rate of subpopulation XR2 (rate constant of reaction v12) was significantly greater than that of XR1 (v9) (Supplementary Fig. 7b). This finding is consistent with the result that the subpopulation of XR2 expressed some cell division-related genes (Fig. 2g). Second, the parameter determining the steepness of the rection rate mediating acquisition of TAM resistance was greater for v10 than that for at in v7 (Supplementary Fig. 7c). This implies that the cell transition from XP to XR2 is more sensitive to cumulative TAM exposure time, which may be caused by the accumulation of epigenetic alterations. The finding is substantiated by the results of single-cell RNA-seq analysis, which showed that the genetic feature of XR2 displayed high expression levels of chromatin-modifying enzymes (Fig. 2g), and pseudotime analysis, in which XR2 was the most differentiated subtype compared with other subtypes (Fig. 3e). These results indicate that parameter inference using our phenomenological model that is based on cell population dynamics can help us to successfully pinpoint underlying mechanisms (differential subpopulation-specific rates) which furthermore are consistent with empirical observations.
Using the fitted parameter sets (Supplementary Table 5) as nominal values, we then performed local sensitivity analysis to examine how much a given change in each single reaction affects the mean-over-time growth rate after the 3rd week (Fig. 4d). The results indicated that v12, the growth rate of XR2, was the most critical factor affecting the mean-over-time growth rate. On the contrary, neither the reverse transition from resistant cell types to XP (or from XP to XS) nor cell death caused by TAM was found to significantly affect the growth rate. Finally, we examined the synergistic inhibitory effect of two key biological processes, cell growth and forward state transition to two different subtypes, on TAM-resistant cell growth (Fig. 4e). Combined inhibition of cell growth rate of both resistant subtypes caused additive growth inhibition (growth rate < 1) at broader inhibitory ranges than the other intervention pairs (Fig. 4e, blue line). However, under the conditions in which the growth of one resistant subtype is repressed, complete inhibition of the transition to another resistant subtype showed stronger regression than complete cell growth inhibition of the same subtype (Fig. 4e, comparison of the gray box). This result indicates that inhibition of cell state transition by, for example, epigenetic inhibitors, has the potential to be more effective than targeting the growth of the resistant subpopulation alone.
Inhibition of molecules mediating the generation of two resistant subpopulations induces regression in the pre-resistant stage
Finally, we experimentally confirmed the hypothesis derived from pseudotime analysis and mathematical modeling that simultaneous intervention of subgroup 4, which undergoes epigenetic alterations via chromatin modification, and subgroup 5, wherein PML acts as an upstream regulator, more effectively inhibits the growth of the resistant cell population. Histone demethylase KDM5B is a candidate epigenetic modulator for subgroup 4, which reportedly modulates resistance to endocrine therapies by increasing transcriptional heterogeneity7. The combination of KDM5 inhibition and PML knockdown suppressed the cell growth of TAM-treated cells (W3, W6, and W9) but not that of cells not treated with TAM (W0) (Fig. 5). Particularly, the growth inhibition in W3, the timing when the cell growth was completely inhibited by TAM (Fig. 1b), indicates that combined inhibition induces a decrease in cell numbers. Together, these results demonstrate that the inhibition of molecules important for the transition to resistant cell subgroups could induce regression before complete acquisition of TAM resistance.