Long-term Multimodal Recording Reveals Epigenetic Adaptation Routes in Dormant Breast Cancer Cells

Abstract

Patients with estrogen receptor-positive breast cancer receive adjuvant endocrine therapies (ET) that delay relapse by targeting clinically undetectable micrometastatic deposits. Yet, up to 50% of patients relapse even decades after surgery through unknown mechanisms likely involving dormancy. To investigate genetic and transcriptional changes underlying tumor awakening, we analyzed late relapse patients and longitudinally profiled a rare cohort treated with long-term neoadjuvant ETs until progression. Next, we developed an in vitro evolutionary study to record the adaptive strategies of individual lineages in unperturbed parallel experiments. Our data demonstrate that ETs induce nongenetic cell state transitions into dormancy in a stochastic subset of cells via epigenetic reprogramming. Single lineages with divergent phenotypes awaken unpredictably in the absence of recurrent genetic alterations. Targeting the dormant epigenome shows promising activity against adapting cancer cells. Overall, this study uncovers the contribution of epigenetic adaptation to the evolution of resistance to ETs.

Significance: This study advances the understanding of therapy-induced dormancy with potential clinical implications for breast cancer. Estrogen receptor-positive breast cancer cells adapt to endocrine treatment by entering a dormant state characterized by strong heterochromatinization with no recurrent genetic changes. Targeting the epigenetic rewiring impairs the adaptation of cancer cells to ETs. See related commentary by Llinas-Bertran et al., p. 704. This article is featured in Selected Articles from This Issue, p. 695.

Figure 1.

Genetic profile of tumor awakening in the clinical setting. A, High-depth profiling (median 105.47×) of ER⁺ breast cancer (BC; estrogen receptor–positive breast cancer) late relapses using a custom targeted panel. The simplified treatment scheme of patients is shown on the left. The heat map shows the mutations in ET resistance drivers in ER⁺ breast cancer passing the filters for allele depth ≥20, Alternate F1R2 + F2R1 ≥4, allele frequency ≥0.1, and consequence level of moderate or high. Time to relapse (years), recurrence in the data set, allele frequency, and relapse site are indicated. Significant genes are indicated based on dN/dS analysis from the q-value of neutrality test at the gene level (*qglobal_cv ≤0.1). B, Clinical histories of patients 2–5. The table shows age and response time to ET for each patient (letrozole). C, Scatter plots of VAF from whole-genome sequencing (WGS) data. Pairwise comparisons were done for pretreatment (diagnostic biopsies) versus progression (surgical biopsies). All patients were managed with primary endocrine therapy until progression. Labeled genes passed two filters: bona fide breast cancer drivers and ET resistance drivers in ER⁺ breast cancer and FATHMM significant score >0.6 (predicted damaging). Detected variants are labeled and color-coded according to detection at diagnosis (teal), progression (magenta), or both (gray). The highlighted gene (TP53) is annotated as a variant detected in ET resistance drivers in ER⁺ breast cancer according to the comprehensive ET-resistance driver gene list compiled based on Bertucci et al. (14). Marginal histograms of VAFs are shown on the sides of each plot. D, Spatial transcriptomics analysis of patients 1–3. On the left, representative images of regions of interest (ROI) from patient 3, pre- and post-treatment, are shown with the relevant staining. Green, pan cytokeratin (CK⁺); yellow, immune cells (CD45⁺); purple, stroma. On the right side, GeoMx UMAPs of previously identified pre-adapted SWNE up and down signatures from (2), and G₂–M checkpoint signatures are shown for patients 1–3 (CK⁺ segment). D1 L biopsy was not suitable for spatial transcriptomics analysis due to poor specimen quality and was excluded from further examinations. S1L: surgical biopsy in the left breast; S1R: surgical biopsy in the right breast; R1R: loco-regional relapse after surgery in the right breast.

Figure 2.

TRADITIOM genetic analysis and lineage composition. A, Cell counts for MCF7 and T47D HYPERflasks for −E2 (circle) and TAM (diamond) arms at their respective time of collection (teal: latency–time between the onset of the treatment and cell cycle arrest in the whole cell population; yellow: dormancy, magenta: awakening, early progression. α−ζ: MCF7 awakening carbon copies, A–F: T47D awakening carbon copies). B, TRADITIOM Live set-up: 12 replicates were seeded in a 48-well plate and imaged two times a week for 150 days using Incucyte Zoom (9 scanning windows per well) to monitor awakening dynamics. Awakening was defined as wells reaching a confluency of 50%. Minor differences in initial plating (violin plot) do not predict awakening times. C, Heatmap of MCF7 high-frequency barcodes (frequency ≥10%) among UT (untreated) samples’ endpoints and TRADITIOM carbon copies (replicates) for both TAM and –E2 arm (AI) at the time of dormancy and awakening (T30–60, AI30–60: dormancy time points, cells treated for 30 or 60 days with tamoxifen or estrogen deprivation, respectively). D, Phylogenetic tree for MCF7 TRADITIOM samples based on SNVs and indels (VAF >10%). Data depict TRADITIOM samples from POT (pretreatment), awakenings (early progression), and TEPs (late progression) of both treatment arms [TAM and AI (–E2)] and their untreated counterpart samples (UT) cultured in parallel for 170 days. E, VAF heatmap for ET resistance driver genes in ER⁺ breast cancer, compiled based on Bertucci et al. (14) derived from TRADITIOM WGS, illustrating only variants with significant changes in VAF across samples (Fisher exact test, P < 0.01). F, Heatmap of T47D high frequency barcodes (frequency ≥10%) among UT (untreated) samples’ endpoints and TRADITIOM carbon copies (replicates) for –E2 (AI) arm at the time of dormancy and awakening. G, Phylogenetic tree for T47D TRADITIOM samples based on SNVs and indels (derived from the targeted sequencing using a custom targeted panel, VAF >10%). H, VAF heatmap for ET resistance driver genes in ER⁺ breast cancer, compiled based on Bertucci et al. (14) derived from the targeted sequencing using a custom targeted panel, illustrating only variants with significant changes in VAF across samples (Fisher exact test, P < 0.01). Data depict T47D TRADITIOM samples from POT (pretreatment), awakenings (early progression), and TEPs (late progression) of both treatment arms (TAM and –E2) and their counterpart untreated (UT) samples cultured in parallel for 180 days.

Figure 3.

Adaptation is driven by divergent transcriptional reprogramming. A, Tamoxifen resistance analysis of the TAM TEPs (late progression) to increasing doses of 4-OHT is depicted in the left. Growth rates of –E2 (AI) TEPs in response to treatment with different drugs: Tamoxifen (Tam, 4-OHT), fulvestrant (Fulv), CDK7 inhibitor (CDK7i), CDK4/6i (palbociclib) are depicted on the right. Representative graphs are shown as normalized confluency fold change upon 7 days of compound treatment (n = 3). B, Principal component analysis (PCA) of bulk RNA-seq expression data for all MCF7 TRADITIOM samples (POT (pretreatment), latency (time between onset of treatment and dormancy), dormancy, awakening (early progression), and TEP (late progression) of both treatment arms [TAM and –E2(AI)] and their counterpart untreated (UT) samples cultured in parallel for 170 days. C, Heat map from bulk RNA-seq data depicting a subset of TRADITIOM MCF7 dormancy signature with top 50 significantly up- and downregulated genes in TAM and –E2 (AI) treated samples during dormancy in comparison with POT (pretreatment). D, GSEA for TRADITIOM dormancy signature.

Figure 4.

The adaptive journey of individual lineages at single-cell level. A, Schematic cartoon of TRADITIOM LSC experimental design. A low-complexity (100 barcodes for MCF7 and 200 barcodes for T47D) barcoded population was seeded in T75 flask format and exposed to ET (TAM and –E2 for MCF7 and –E2 for T47D) in a nonperturbed system (no serial passaging). Cells were collected at the indicated time points (duplicates of T0, 1 month, 2 months, and 4 awakening samples for MCF7 and duplicates of T0 and 1 month for T47D) and analyzed by scRNA-seq. Cells were imaged once a week till awakening with 108 scanning windows covering each flask. B, Proliferation dynamics of TRADITIOM LSC for a representative sample [AI (−E2) aw1] determined by weekly imaging along with the local confluency heat map at the time of collection (awakening). Magenta lines follow the growth dynamics of main awakening areas from time zero (T0) to overt expansion that represents the collection point for scRNA-seq analysis. Violin plot shows the confluency of each scanning area at T0 (onset of estrogen deprivation). The T0 confluency of the awakening area is highlighted in purple. C, Heat map of winner barcodes’ frequency for TRADITIOM LSC carbon copies for –E2 (AI) arm at the time of awakening, in the POT population (pretreatment), at the start of the experiment (T_0, time zero) and at early [1 month (1 mo)] and late dormancy [2 months (2 mo)] stages derived from either genomic barcode sequencing (g) or scRNA-seq (sc). D, UMAP projections of MCF7 TRADITIOM LSC –E2 (AI) arm of T0, early [1 month (1 mo)] and late dormancy [2 months (2 mo)] and awakening samples. E, Adaptive journey of winner and non-winner lineages (others) of each TRADITIOM LSC –E2 (AI) awakening sample from T0 to awakening (early progression). Pie charts depict the occupancy of miscellaneous UMAP clusters for each lineage. The approximate awakening time of each carbon copy (replicate) is shown with arrows. F, UMAP projections of T47D TRADITIOM LSC samples at T0 and dormancy [1 month (1 mo)]. G, Dot plot indicates similarity of transcriptional space occupied by MCF7 cells under −E2 treatment (AI) with those under TAM (tamoxifen) and with T47D counterpart under −E2 treatment. Marker genes for clusters in T0, dormancy, transition, and awakening from MCF7 TAM and T47D −E2 UMAPs were checked for enrichment in cluster marker genes from MCF7 −E2 (AI) UMAP. Size of the dot represents the significance of enrichment in context to –log₁₀(P value; correction method: FDR, P ≤ 0.05). The color of the dots corresponds to the respective color of clusters at T0, dormancy, and awakening annotated by colored bars next to cluster names on the x and y axes (T0: gray, dormancy: yellow, awakening: magenta).

Figure 5.

Failed awakenings. A, TRADITIOM Live awakening topography analysis depicting the awakening dynamics. Twelve carbon copies (cc; replicates) were seeded in a 48-well plate and imaged two times a week for 150 days using IncuCyte Zoom (9 scanning windows per well) to monitor awakening dynamics. Awakening was defined as wells reaching a confluency of 50% (FA: failed awakening; D: dormant; AW: awakening; GA: global awakening; LA: localized awakening; 1LA: 1 localized awakening; >1LA: more than 1 localized awakening). B, IncuCyte time-lapse images from a failed and a bona fide awakening are shown as an example. C, Proliferation dynamics of a representative MCF7 [AI (−E2) awakening 2] and T47D TRADITIOM LSC sample determined by weekly IncuCyte time-lapse imaging until awakening. Pink lines indicate failed awakenings. Magenta lines follow the growth dynamics of main awakening areas. D, MCF7 Geminin-mCherry NLS-GFP cells were treated with estrogen deprivation (−E2) for 3 months to establish a detailed understanding of long-term dormancy–awakening dynamics. Image sets were analyzed using Essenbio Sartorius software from daily imaging. 35% of replicates (n = 60) had dormant persister cells/small colonies until day 88. The proportion of Geminin-mCherry-positive (S/G₂–M; pink lines) is indicated normalized to total count which was quantified by NLS-GFP (gray lines). E, Distribution of Ki-67 expression levels in winner and non-winner (others) lineages associated with either G₁ or S/G₂–M states across T0, dormancy, and awakening samples of MCF7 TRADITIOM LSC cell-cycle regressed data set. F, Cell-cycle regressed UMAPs of the subset of cells in S/G₂–M state of TRADITIOM MCF7 LSC data set. G, Pathway enrichment analysis comparing failed awakenings (S/G₂–M cells in dormancy: cluster 0, 3) to bona fide awakenings (S/G₂–M cells of awakening winner lineages: cluster 4).

Figure 6.

Targeting the dormant epigenome. A, Clustered heat maps of histone posttranslational modifications of super-SILAC mass spectrometry for TRADITIOM MCF7 and T47D samples [time zero (T0), latency (time between treatment onset and dormancy entry), dormancy, awakening (early progression), and TEPs (late progression)]. Significantly enriched (dormancy 30 days vs. TEPs, two-tailed t test: *, P < 0.01; **, P < 0.001; ***, P < 0.0001) modifications are depicted in bold, and the ones found to be associated with dormancy are highlighted in yellow. B, Schematic representation of small-molecule inhibitor experiments. Inhibitors against G9a (H3K9me2), EZH2 (H3K27me3), and KMT5B/C (H4K20me3) were used either alone or in combination. Start time of the inhibition was either at the beginning of estrogen deprivation to target persister pool generation or at 30 days of estrogen deprivation (dormancy) to target established dormant cells. C, Proliferation dynamics of MCF7 cells in E2-supplemented conditions (+E2) after treatment with inhibitors against EHMT2, EZH2, KMT5B/C, dual combinations of each and vehicle. D, Proliferation dynamics of MCF7 cells in estrogen-deprived conditions (−E2) after treatment with inhibitors against G9a, EZH2, KMT5B/C, dual combinations of each and vehicle. Proliferation dynamics of T47D cells in E2-supplemented (+E2; E) and deprived (−E2) conditions (F) after treatment with inhibitors against G9a, EZH2, KMT5B/C, dual combinations of each and vehicle (one-way ANOVA with Dunnett correction: *, P < 0.05; **, P < 0.01; ***, P < 0.001; ****, P < 0.0001). Error bars represent standard deviation (n = 3). G, Relapse-free survival (RFS) curves for ER⁺ breast cancer patients stratified based on the expression of the epigenetic dormancy signature (high vs. low EHMT2/EZH2/KMT5C expression). Left: no adjuvant treatment; middle: adjuvant endocrine therapy (TAM/AI); right: AI adjuvant treatment. Multivariate analysis for clinically relevant prognostic biomarkers is shown in the onset table. H, Proliferation dynamics of MCF7 dormant cells (pretreated for 30 days with –E2) after treatment with inhibitors against G9a, EZH2, KMT5B/C, dual combinations of each and vehicle (one-way ANOVA with Dunnett correction: *, P < 0.05; **, P < 0.01; ***, P < 0.001; ****, P < 0.0001). Error bars represent standard deviation (n = 3). I, Model: endocrine therapy-induced dormancy is characterized by a consistent epigenetic reprogramming involving a global increase in histone repressive marks (H3K9me2, H3K27me3, and H4K20me3). The dormant epigenome is unstable and through a progressive loss of the histone repressive marks (erosion), cells resume proliferation in a process that mimics patient relapse (awakening). Epidrugs (G9a/EZH2/KMT5B/C inhibitors) can interfere with epigenetic reprogramming and block the formation of persister dormant clones. During adaptation, dormant cells engage in sporadic cycling (failed awakening) while under therapeutic stress possibly forcing cells into a subsequent round of epigenetic reprogramming that could also be antagonized with epidrugs.