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 (Figure 1a). Treatment with 1 µM TAM inhibited cell growth (Figure1b) through decreasing the cells in S phase whereas increasing the cells in G1 phase (Figure 1c and 1d), suggesting that TAM treament induced G1 arrest. The growth of cells was almost completely inhibited until week 5 but recoveered thereafter (Figure 1b). The cell cycle population of TAM-treated cells was also dysregulated until week 5 but became comparable to control cells after week 6 (Figure 1c and 1d). These results showed the process by which cells survived and restored growth potential in the absence of ER signaling.
We next analyzed the bulk RNA-seq data of TAM-treated MCF-7 cells to identify the changes in gene regulatory network and compare with the changes 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 acquisition 15. 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 (Supplementary Figure 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 Figure 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 by principal component analysis (PCA). The PC1 axis classified the expression profiles of all DEGs into two groups (Figure 1f). The separation in log2FC values was the greatest between weeks 4 and 5, which preceded the recovery of cell growth for 1 week to 2 weeks (Supplementary Figure 3). In addition, the separation in log2FC values between weeks 1 and 5 was greater than that between weeks 6 and 12, indicating that gene expression became more stable at the later stages.
We then investigated the relationships between 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 (Figure 1f). The enrichment analysis of each group was carried out using the Reactome pathway (Figure 1g) and KEGG (Figure 1h) databases. Cluster A genes showed a rapidly decreasing expression pattern compared with cluster F genes. Cluster A was enriched in genes related to both receptor tyrosine kinase signaling and fatty acid metabolism. On the contrary, clusters C and D consisted of genes showing consistent up-regulation. Genes in cluster C responded early to TAM than genes in cluster D, and some of these cluster C genes were related to TGFβ signaling or extracellular matrix–receptor interaction, such as TGFB2, SMAD3, and MMP9, whereas others 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 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. 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, and autophagy. The interferon and FoxO signaling pathways exhibit anti-survival functions in cancer cells exposed to anti-cancer agents16, whereas autophagy contributes to cell survival under normal conditions17, 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 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 resistance 18, 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 (before starting TAM treatment), week 3 (beginning of complete cell growth inhibition), week 6 (end of complete cell growth inhibition), and week 9 (at the acquisition of TAM resistance) (Figure 2a). RNA-seq analysis of 1,108 single cells yielded 577 high-quality single-cell gene expressions (Figure 2a; see the Materials and Methods section). The Pearson correlation coefficient of 11,413 gene expression values between individual cells decreased at week 3 and then gradually increased at weeks 6 and 9 (Figure 2b). This changing pattern of correlation coefficient suggests 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 software19 to correct for the bias caused by the difference in the cell cycle stage. PI staining showed the same trend in bulk cell populations (Figure 1c), indicating that the percentage of cells in the S phase and G2/M phase decrease only when cells cannot grow (Figure 1b and 2c). All cells were mapped on a three-dimensional UMAP plot and projected in two dimensions (Fig. 2d and Supplementary Figure 4). Single-cell data were roughly divided into two groups: week 0 and the others (weeks 1–9). Cells were widely distributed in space at weeks 3 and 6 but were localized in two separate regions at week 9. These cells could be clustered into six subpopulations in the UMAP plot (Figure 2e). Cells in subgroups 1 and 6 belonged to the week 0 group, and these cells were almost diminished at week 3. By contrast, cells in subgroups 2 and 3 newly emerged at week 3. Finally, these cells were converted into two groups: one containing subgroup 4, and the other containing subgroup 5.
Cluster-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 Figure 5a). Subgroups 1 and 6 were the major groups at week 0, and marker genes in these subgroups included typical ER pathway target genes such as AREG20 and GREB121. This result clearly showed that transcription activity of the ER was down-regulated in the other subgroups. On the other hand, marker genes in other subgroups, especially subgroup 2, were rather nonspecific. 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 (Figure 2f and Supplementary Figure 5b). Subgroups 1 and 6 contained highly expressed gene modules 1, 2, and 12, which function on estrogen-dependent gene expression, unfolded protein response, and amino acid and nucleotide metabolism. On the other hand, gene expression module 5 was relatively low. Subgroups 2 and 3 were the major subpopulations in weeks 3 and 6. Both these subgroups showed high expression levels of genes in module 5, some of which are involved in interferon signaling, TGFβ signaling, and tight junctions. These enriched terms were quite similar to the early responsive cluster C in the bulk RNA-seq experiment (Figure 1g and h). Subgroup 4, whose population was increased at week 9, showed high expression levels of genes in modules 4, 6, and 7, 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 6. 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 week 9, contained highly expressed genes in modules 8–11. 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 11 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. On the basis of the aforementioned results, we uncovered 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 the migratory phenotype with high expression levels of adhesion molecules and RTK signaling-related genes (subgroup 4).
Trajectory analysis of TAM resistance
To confirm the cell transition trajectory into two different types of resistant subgroups, we conducted pseudotime analysis (Figures 3a–c). The pseudotime of each cell calculated from the gene expression data was correlated with the sampling time after starting the continuous TAM treatment (Figure 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 (Figure 3e). To estimate important molecules involved in the emergence of subgroups 4 and 5, we analyzed DEGs along with an estimated cell trajectory. A total of 273 and 79 genes were detected as highly expressed genes in subgroups 4 and 5, respectively (Figure 3f and g). Then, we analyzed the upstream regulators of these genes using the ChIP-Atlas database22 (Figure 3h and 3i). Upstream analysis of the trajectory to subgroup 4 (Figure 3f and 3h) showed that only 4 of the top 10 factors represented ChIP-seq data obtained from MCF-7 cells, and most of the others were obtained from the triple-negative breast cancer (TNBC) cell line (Figure 3h, shown in green). These results also suggest that most of upregulated genes in subgroup 4 are controlled by bromodomain-containing proteins, BRD4 and BRD2, which recognize acetylated histones and act as super enhancers 23,24. 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 (Figure 3j, top and middle), as shown in bulk RNA-seq data (Figure 1). 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 (Figure 2f).
In subgroup 5, 7 of the top 10 factors were obtained from ER-positive breast cancer or normal cells (Figure 3i), suggesting that subgroup 5 retained the transcriptional network of parental MCF-7 cells. ChIP-seq analysis using anti-WDR5 antibody showed the best q-value and fold enrichment score. WDR5 has been associated with histone methylation25. Although WDR5 was not up-regulated in bulk RNA-seq data (Figure 3j, bottom), other histone methylases might regulate genes related to subgroup 5. Other up-regulated candidates estimated from ChIP-Atlas database included MYC, TAF1, PML, and BRD4. Among these genes, MYC, TAF1, and PML were up-regulated before week 4, as shown by bulk RNA-seq analysis (Figure 3h, bottom). These data indicate that MYC, TAF1, and PML1 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
Finally, we constructed a mathematical model that reproduces the process of acquiring TAM resistance, based on cell trajectories obtained by pseudotime analysis, to estimate the core processes that contribute to the growth and differentiation of resistant cell populations (Figure 4a). This model consists of four major cell subpopulations: cells initially sensitive to TAM (XS, clusters 1 and 6 in Figure 3e), pre-resistant cells (XP, clusters 2 and 3), resistant cells showing high expression of carbon metabolism-related genes (XR1, cluster 5), and resistant cells with highly adhesive and invasive phenotypes (XR2, cluster 4). We hypothesized that TAM induces the death of XS and XP cells but not that of resistant cell populations and sigmoidally promotes cell state transition from XS to XP and from XP to XR1 or XR2 in response to a period of TAM treatment. This hypothesis assumes that cell state transition is caused by the accumulation of genetic or epigenetic changes induced by continuous TAM treatment. We obtained 20 model parameter sets, which could reproduce two experimental time-course datasets, total cell growth rate in the presence of TAM (Figure 1b) and rate of cell subpopulation (Figure 3e), simultaneously (Fig. 4b and c and Supplementary Figure 6a) using the BioMASS computational framework26.
Comparing each parameter set, we found two remarkable features of the well-fitting parameter sets. First, the growth rate of subpopulation XR2 (rate constant of parameter v12) was significantly greater than that of XR1 (v9) (Supplementary Figure 6b). This finding is consistent with the result that the subpopulation of XR2 expressed some cell division-related genes (Figure 2f). Second, a parameter determining the steepness of sigmoidal change in v10 during the time course was greater than the steepness of sigmoidal change in v7 (Supplementary Figure 6c). This implies that the cell transition speed from XP to XR2 showed switch-like alteration triggered by the accumulation of epigenetic alterations. This 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 (Figure 2f), and pseudotime analysis, in which XR2 was the most differentiated subtype compared with other subtypes (Figure 3e). These results indicate the potential of a constructed model and obtained parameter sets in reproducing not only the growth rate and ratio of cell populations during the time course but also the biological features of subpopulations.
On the basis of the model and parameter sets, we then performed sensitivity analysis of each single reaction to examine whether a change in each single reaction affects the size of the resistant cell population. We focused on the integral growth rate after week 3, which determines the total increase in cell number after the acquisition of TAM resistance, and calculated the sensitivities of each reaction on that (Figure 4d). The results indicated that v12, growth velocity of XR2, was the most critical phenomenon for the increase in resistant cell numbers. On the contrary, neither the reverse transition from resistant cell types to XP (or from XP to XS) nor cell death caused by TAM could be a negative regulator of the increase in resistant cell numbers alone. Thus, we showed that the increase in TAM-resistant cells is strongly controlled by the growth rate of XR2. Finally, we examined the synergistic inhibitory effect of two biological events in the model, cell growth and forward state transition to two different subtypes, on TAM-resistant cell growth (Figure 4e). Combined inhibition of cell growth rate of both resistant subtypes caused synergistic regression (growth rate <1) at broad inhibitory ranges than the other conditions (Figure 4e, bottom right). However, combination inhibition including at least one transition to the resistant subpopulation induced potent regression when the cell transition was completely reduced. 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 itself.