Cellular adaptation to cancer therapy along a resistance continuum

a, Drug-adapted lines generated for Ovsaho, a BRCA-mutant ovarian cancer cell line, using dose-escalation from 1 to 40 μM of olaparib. b, Cell viability assays (9 days of treatment) for the Ovsaho adapted lines showing progressive increase in resistance. Bar plots show the average of 6 technical replicates with standard deviation (SD) error bars. c, Same as in a for COV362, also a BRCA-mutant ovarian cancer cell line. d, Same as in b for COV362. e, Drug-adapted lines generated for A375, a melanoma cell line, using dose-escalation from 30 nM to 4 μM of dabrafenib. f, Cell viability assays (4 days of treatment) for the A375 adapted lines. g, Drug-adapted lines generated for PC9, a lung adenocarcinoma cell line, using dose-escalation from 8 nM to 1.2 μM of osimertinib. h, Cell viability assays (5 days of treatment) for the PC9 adapted lines. Cell viability assays were performed using CellTiter Glo 2.0 and values are normalized by DMSO control. Data are shown as means and SDs of 6 technical replicates.

a, Generation of drug-tolerant persisters for Kuramochi cell line by treatment with 320 μM of olaparib for 9 days. The surviving cells were allowed to recover for 9 days without treatment. Two rounds of treatment-recovery experiments were performed in triplicates and the viability of the recovering cells was tested. Cell viability curves are shown as the average of three replicates (n = 6 technical replicates) and error bars as SD. The bar plot shows no significant IC₅₀ increase for persisters relative to the parental population (C vs P320.1: P = 0.03; ns, not significant; two-tailed t-tests). b–d, Extended treatments at constant doses. Representative Kuramochi adapted lines (T5, T20 and T80) were maintained on olaparib treatment at the final concentrations to which they had adapted (5, 20 and 80 μM, respectively) for the number of days that they had taken to adapt to their subsequent populations at increased concentrations (T10, T40, T160, respectively) (see Fig. 1a), generating the extended T5-E, T10-E and T80-E populations (three independent replicates). Prolonged exposure at constant doses did not generate similarly resistant populations for the T5 and T20 extended treatments when compared to their dose-escalated counterparts (T10 and T40, respectively). The extended treatments did however lead to higher resistance relative to their parental populations. The cells from the T80 extended treatments were similarly resistant to T160, highlighting the role of treatment duration. Together, these observations indicate that both dose and time influence the fitness outcomes along the adaptive path. Cell viability was measured by CellTiter Glo 2.0 (each biological replicate with 6 technical replicates) and by cell counting (bar plots show the average of three biological replicates and error bars as SD) after 9 days of treatment at their respective escalated doses (10, 40 and 160 μM). IC₅₀ (left bar plots) and relative cell counts (right bar plots) were compared with their parental and dose-escalated populations (significant P-values are shown; ns, not significant; two-tailed t-tests). Dose response curves for the original T5, T10, T20, T40, T80 and T160 are from Fig. 1b as performed in 3 independent experiments. e–h, To test whether dose escalation facilitates drug adaptation in our model, we compared the resistance phenotypes between cells that were dose-escalated from 1 to 10 μM with drug-naive cells directly exposed to the same doses for the same duration. Dose-escalated cells were able to proliferate and reach confluency, contrasting with drug-naive cells directly exposed to higher doses (5 and 10 μM). The continuous long-term 10 μM treatment led to very few viable cells, which prevented us from testing their fitness with viability assays (h). The populations that were dose-escalated up to 2.5 μM were significantly more resistant than their continuously treated counterparts. Plots indicate dose-response curves for three biological replicates (R1, R2, R3; 6 technical replicates each) and bar plots of average IC₅₀ values (3 biological replicates with SD error bars). Bar plots (right) showing the average cell counts in f-h for the 3 biological replicates with SD error bars. Cell viability assays and cell counts (treated with 2.5, 5 and 10��μM of olaparib, respectively) after the treatment periods are shown and compared with the dose-escalated populations (P-values are indicated; ns, not significant; two-tailed t-tests). Microscopic photos show that T10 (dose-escalated) cells reached confluency while T10.C (continuous) did not (scale bar 50 μm).

a, Spearman’s correlations across the averaged transcriptomes of the Ovsaho olaparib-adapted lines. Untreated control (C) and T5 to T40 were collected for scRNA-Seq. b-c, UMAP of Ovsaho scRNA-seq data on the adapted lines colored by treatment (b) and by cluster as determined by Louvain clustering (c), highlighting the emergence of distinct cellular states. d, Bar plot representing the percentage of cells from each condition in each cluster, indicating cell state frequency changes across the adapted lines. e, Differentially expressed genes across Ovsaho clusters, highlighting the expression of shared genes with Kuramochi (Fig. 2a), COV362 (Extended Data Fig. 3j) and PDX resistance models (Fig. 4f). Color bar represents the scaled average expression for a particular cluster and circle sizes are the fraction of cells expressing the gene. Overall, genes highlight dedifferentiation and metabolic reprogramming along the resistance continuum. f, Same as in a for COV362. g-h, Same as in b and c for COV362, respectively. i-j, Same as in d and e for COV362, respectively. k, Drug-tolerant persisters generation by treatment of drug-naive cells with 10 and 320 μM of olaparib for 9 days. The adapted T10 and T320 lines were maintained at the same doses, respectively (left, scale bar 50 μm, representative of two biological replicates). UMAP of the scRNA-Seq data for the persister and adapted cell populations (right). Phenotypic features of persisters are shown in Supplementary Fig. 2. l, Differentially expressed genes across conditions (relative to C) highlighting shared and distinct marker gene expression among the persister and adapted lines.

a, Spearman’s correlations across the averaged transcriptomes of the A375 (melanoma) adapted lines that evolved resistance through dose-escalation from 30 nM to 4 μM of dabrafenib (Extended Data Fig. 1e,f). b, UMAP of scRNA-seq data on the individual lines. Colors and numbers indicate subpopulations determined by Louvain clustering. c-d, UMAP for all cells colored by treatment (c) and by cluster as determined by Louvain clustering (d), highlighting distinct transcriptional states across the adapted lines. e, Bar plot representing the percentage of cells from each condition in each cluster, indicating the emergence of cell states and frequency changes that suggest state selection along the resistance continuum. f, Differentially expressed genes across A375 clusters, highlighting the expression of lineage plasticity and metabolic reprogramming markers across the resistance continuum. Differentiation markers indicate multiple intermediate states, showing a progression from SOX10^high (clusters 0-3) to SOX10^low expression (clusters 4-7) and loss of MHC class II genes (CD74, HLDA-DRA/B) (Supplementary Table 4). Diverse lineage states arise along the resistance path, such as cells expressing both melanocytic (DCT, MITF) and neural-crest (NGFR) lineage markers (cluster 3); a population expressing highly MITF, RAP1B and TMSB4X (cluster 4); cells expressing neuronal markers (ROBO1, BMP5, SYT4) (cluster 5); and finally a mesenchymal-like invasive state (SOX9, FN1) (clusters 6-7, Supplementary Table 4). A gradual increase in the expression of ATF4 and NRF2 regulated genes toward the most resistant states (clusters 4 to 7), indicating extensive metabolic rewiring. Color bar represents the scaled average expression for a particular cluster and circle sizes are the fraction of cells expressing the gene. g-h, Same as in a and b for PC9 (lung adenocarcinoma), respectively. Adapted lines evolved resistance through dose-escalation from 8 nM to 1.2 μM of osimertinib (Extended Data Fig. 1g,h). A persister population (P1.2) was generated by treatment with 1.2 μM of osimertinib for 9 days. i-l, Same as in c-f for PC9, respectively, including a persister population (P1.2, i). PC9 cells transitioned to a KRT6A^low low expression and coexist in a combination of CD44^highCD24^low and CD44^lowCD24^high states (l) (Supplementary Table 5).

a, Scatter plot showing the average expression of CD24 and CD44 markers in each of the identified subpopulations colored by treatment condition (left) and state (right). b, Representative plots of the flow cytometry analysis of CD24 and CD44 staining across the Kuramochi adapted populations. Bottom left quadrants represent double negative staining for both markers. A transition from CD24^highCD44^low (top right) to CD24^lowCD44^high (bottom right quadrants) is observed. Numbers indicate the percentage of cells in each quadrant. Colored squares indicate the subpopulations sorted for the cell viability assays shown in (c), ATAC-Seq (Fig. 3) and CRISPR-based screen (Fig. 5). c, Cell viability assays of sorted subpopulations indicating differential resistance within (T10) and between (T10, T40, T320) treatment conditions. Sorted subpopulations were labeled according to their CD24/CD44 profiles. Bar plots show the IC₅₀ averages of two independent experiments (6 technical replicates each).

a, Western blot of the induced SOX17 in the T320 line compared to the parental control (C) cells (three technical replicates). Bar plot shows the quantification signal. Full western blot shown in Supplementary Fig. 6. b, Schematic of the constitutive SOX17 overexpression construct (top). Quantification by live cell imaging of the percent area covered by cells over a 20 day time period of T320 parental and T320 with SOX17 overexpression on DMSO (left) and on olaparib (5 μM) treatment (right). Average and SD of three biological replicates (three technical replicates each) shown. c, Western blot of SOX17 shRNA knockdown in parental control (C) cells compared to non-targeting scrambled shRNA (three technical replicates). Bar plot shows the quantification signal. Full western blot shown in Supplementary Fig. 7. d, Quantification of SOX17 shRNA knockdown in parental control (C) cells and adapted T5, T40 and T320 by live cell imaging (three biological replicates, three technical replicates each) of the percent area covered by cells normalized by DMSO controls on day five of treatment of 2.5 μM olaparib. Standard two-sided t-test was performed by comparing with the negative control cells (Neg. Ctrl) (SOX17 #1: P = 0.01; SOX17 #2: P = 0.02; SOX17 #3: P = 0.05; T40: P = 0.004; T320: P = 0.008; ns: not significant). e-i, EMT signature scores were calculated across treatment conditions for each resistance model. The set of genes annotated as “Hallmark epithelial-mesenchymal transition” from the MSigDB was used to calculate the scores with the UCell algorithm. Violin plots represent the score distributions with the average plus and minus SD. Numbers of cells are: (e) C (N = 453), T1 (N = 1223), T2.5 (N = 760), T5 (N = 398), T5 (N = 398), T10 (N = 1030), T20 (N = 640), T40 (N = 918), T80 (N = 910), T160 (N = 1348), T320 (N = 568); (f) C (N = 840), T5 (N = 713), T10 (N = 704), T20 (N = 731), T40 (N = 650); (g) C (N = 1148), T5 (N = 180), T10 (N = 278), T20 (N = 447), T40 (N = 272); (h) C (N = 442), T003 (N = 478), T006 (N = 906), T0125 (N = 306), T0250 (N = 733), T05 (N = 493), T1 (N = 605), T2 (N = 641), T4 (N = 455); (i) C (N = 169), P1.2 (N = 673), T0008 (N = 329), T0016 (N = 528), T0032 (N = 414), T0075 (N = 505), T0150 (N = 329), T0300 (N = 442), T0600 (N = 424), T1.2 (N = 395). j, Representative images (biological duplicates) of stained membranes from a transwell cell invasion assay. Colors represent the untreated cells (C), adapted lines T5 and T320 and the Slug-lines on different doxycycline (dox) concentrations (C-Slug₀ as negative and C-Slug₁ on 1 μg/ml of dox). Quantification of transwell cell invasion assay data of three images per cell line (right). One way ANOVA comparing Slug-induced (SNAI2) lines to the drug-naive (C) cell line. k, Quantification of induced SNAI2 gene expression from qRT-PCR (two independent runs in triplicates). One way ANOVA comparing Slug-induced lines to the drug-naive (C) cell line. l, Representative western blot (three technical replicates) of the induced Slug lines with different doxycycline concentrations ranging from 0.01 μg/ml to 1 μg/ml. We observed that C-Slug₁ produced a relevant phenotype, which was used for further experiments. Bar plot shows the quantification signal. Full western blot shown in Supplementary Fig. 8. m, Quantification by live cell imaging of DMSO controls For Slug-induced lines. Data represents percent area covered by cells over a 54 day time period. Average and SD of three biological replicates (three technical replicates each) are shown.

a, Heat maps depicting scRNA-Seq inferred copy number alterations (CNAs) for Kuramochi adapted and persister lines relative to untreated control (C). Gene expression signal is scaled along the chromosomes as calculated by inferCNV. Red represents genomic regions with inferred copy number gain and blue represents regions with inferred copy number loss. Colors and labels on the left represent the subpopulations identified by transcriptome clustering (Fig. 1d) showing correspondence between inferred CNAs and transcriptional states. b, Enrichment of genes within CNA regions identified by whole-exome sequencing overlapping with the genes present in the transcriptional modules detected for the adapted lines (Fig. 2a) in respect to random distributions (P-values determined by 10,000 permutations, two-sided permutation tests). Blue and red lines indicate the observed overlap of copy number loss and gain, respectively, and the gray lines indicate the significance threshold in an expected random distribution. c, Inferred CNAs defined for the adapted Ovsaho and COV362 lines relative to their untreated controls. Color annotations represent sample and transcriptional clusters identified in Extended Data Fig. 3c,h.

a, Correlation between the number of differentially accessible peaks and total number of detected peaks per sample. The lack of a positive correlation indicates that the number of differentially accessible peaks are not explained by peak detection bias (Spearman’s rho = −0.62, P = 0.07, two-sided Spearman Rank Correlation test). b, Heat maps showing the hierarchical clustering (rows) of the enriched TF motifs in the differentially more (green) and less (purple) accessible peaks (relative to C). Color bar represents the scaled p-values for TF motif enrichment (only transcription factors with P-value < 1e-5 were included, one-sided hypergeometric test, Benjamini-Hochberg corrected). The analysis controlled by peaks not overlapping CNAs (“no-CNA”) reveals the same patterns when compared to all peaks. c, Enrichment of differentially accessible ATAC-Seq peaks overlapping CNA regions (relative to the untreated control C) identified by whole-exome sequencing for representative adaptive lines (T10, T40 and T320) (P-values determined by 1,000 permutations, two-sided permutation tests). Blue and red lines indicate the observed overlap of copy number loss and gain, respectively, and the gray lines indicate the significance threshold from the expected random distribution. The proportion of peaks not overlapping CNAs ranged from 70% for less, and 80–90% for more accessible peaks, highlighting widespread global chromatin changes. d, Transcription factor (TF) motif enrichment analysis controlling for peaks not overlapping with CNAs, showing similar frequency of peaks containing TFs motifs involved with regulation of stress response (AP1, NRF2 and ATF4) and ovarian cancer TFs (SOX17, PAX8 and WT1). Bar plots represent the percentage of differentially accessible peaks for each condition (relative to C). Color bars represent the motif enrichment P-values (one-sided hypergeometric test, Benjamini-Hochberg corrected). e, Proliferation of C, T320 and T320W (which is the T320 population after drug removal for 120 days) on DMSO, 10, 40 and 320 μM of olaparib for 8 days. The T320W exhibits partial reversibility in resistance, especially at 320 μM. Proliferation was measured every 2 days using Cell TiterGlo 2.0 assay and values are normalized from day 0 (24 h after cell seeding). Data represents the average (dots) of 6 technical replicates and error bars indicate the SD. f, Box plots of the signal for less (purple) (N = 6771) and more (green) (N = 7332) accessible peaks shown in Fig. 3f. This shows that the signal across the less accessible regions (relative to C, Fig. 3d) both with and without treatment are overall stable when comparing T320 and T320W. For the more accessible regions, the signal decreases in T320W without treatment and recovers its accessibility during drug treatment. P-values were determined from two-sided Wilcoxon tests and the effect sizes are Cohen’s d. Boxes represent the interquartile range (IQR) and centre lines represent the median. Whiskers extend to the smallest and largest values within 1.5 times the IQR from the lower and upper quartiles, respectively. Data points beyond the whiskers are outliers.

a, Tumor growth curves for PDXs treated with four doses of talazoparib (shades of green) and vehicle controls (gray) in the persistence experiment. Each condition was collected in duplicates for scRNA-Seq. Two replicates of each condition were generated (indicated by “.1” and “.2”). b-d, Tumor growth curves for PDXs in first (b), second (c) and third (d) rounds of talazoparib treatment in the resistance experiment. A drug holiday in the first round is indicated in red. Two replicates for each condition were collected for scRNA-Seq: vehicles (gray), round 1 (solid pink) and round 3 (blue, with an additional vehicle treated R3.V). A replicate tumor from round 1 (dashed pink) was seeded in a new cohort to generate the round 2 treatment. e, Hematoxylin and eosin (H + E) staining of the entire section of V1, D3.1 (persistence) and R3.2 (resistance) (column 1) and an enlargement (column 2). Histology was annotated by a gynecologic pathologist. The V1 tumor has the classic appearance of serous carcinoma with nests of neoplastic cells with “slit-like spaces” separating strands of tumor cells. The D3.1 tumor still is recognizable as serous carcinoma while more poorly differentiated than V1. The R3.2 tumor has a biphasic pattern with some areas that retain a recognizable papillary pattern, while it is dominated by poorly differentiated/undifferentiated tumor structures with discohesive tumor cells. DAPI staining (blue), CD44 (green) and SOX17 (red) as single channels (columns 3-5), merged (column 6) and an enlargement (column 7). Consistent with scRNA-Seq data, the well-differentiated regions in V1 highly express SOX17 and lack expression of CD44, which is also observed in D3.1. Regions of similar histology in R3.2 have weaker expression of SOX17 and remain to lack expression of CD44 (columns 4-5). In contrast, the poorly differentiated regions in R3.2 strongly express CD44 and weakly express SOX17 (columns 6-7), thus indicating heterogeneous populations in different resistant states. Scale bars represent 100 μm except in column 1 (1000 μm). f, Heat map of the scaled gene average expression across samples of the dose-response experiment (persistence), highlighting differentially expressed genes that show a dose-dependent trend. g, Inferred CNAs defined for the persistence and resistance PDX experiments relative to their respective vehicle treated tumors. Color annotations represent sample and transcriptional clusters identified in Fig. 4b,d.

a, Rank plots showing the score differences for C_D, T10, T320 and T320_D. Higher absolute score differences indicate increased gene depletion relative to C. Highlighted genes are significantly depleted (FDR < 0.05, two-sided likelihood ratio test (LTR)) at the endpoint. Only genes with negative score differences are shown, indicating gene depletion. b, Overlap between the significantly depleted genes with high absolute score differences (>1 SD) relative to C. P-values were determined by one-sided hypergeometric tests. c, Heat maps of the differentially abundant metabolites (P < 0.05, two-sided t-tests) comparing olaparib sensitive (C) and resistant (T320) populations. Conditions represent cells on olaparib treatment (320 μM for 24 h) or cells without treatment. The experiment was performed in triplicates for each condition. Color scale represents scaled normalized metabolite abundance. Metabolites involved in nucleotide biosynthesis are highlighted in red. d, Dose-response curves of olaparib treatment or in combination with CB-839 after 6 days for representative populations within the resistance continuum (C, T5, T10, T40 and T320). Data points are normalized relative to vehicle-treated controls (for each respective line) and represent the average of 3 independent experiments (6 technical replicates per experiment) and SD. e, Bar plots showing the IC₅₀ determined from dose-response curves (d) after 6 days of treatment with olaparib, or the combination of olaparib with CB-839. Dots represent three independent experiments (average of 6 technical replicates). Cells were sorted based on CD24^low/CD44^high profiles (as in Extended Data Fig. 5b,c). f, Bar plot showing the average cell viability of cells treated with CB-839 (50 nM) normalized to vehicle controls (3 independent experiments, 6 technical replicates each). g, Heat map of the bliss synergy scores⁴¹ of the olaparib and CB-839 treatment combination. h, Representative western blot (biological duplicates) showing the expression of NRF2 and its target HO1 (HMOX1) with and without olaparib treatment. Full western blot is shown in Supplementary Fig. 10.