Diverse clonal fates emerge upon drug treatment of homogeneous cancer cells

Abstract

Even among genetically identical cancer cells, resistance to therapy frequently emerges from a small subset of those cells^1-7. Molecular differences in rare individual cells in the initial population enable certain cells to become resistant to therapy^7-9; however, comparatively little is known about the variability in the resistance outcomes. Here we develop and apply FateMap, a framework that combines DNA barcoding with single-cell RNA sequencing, to reveal the fates of hundreds of thousands of clones exposed to anti-cancer therapies. We show that resistant clones emerging from single-cell-derived cancer cells adopt molecularly, morphologically and functionally distinct resistant types. These resistant types are largely predetermined by molecular differences between cells before drug addition and not by extrinsic factors. Changes in the dose and type of drug can switch the resistant type of an initial cell, resulting in the generation and elimination of certain resistant types. Samples from patients show evidence for the existence of these resistant types in a clinical context. We observed diversity in resistant types across several single-cell-derived cancer cell lines and cell types treated with a variety of drugs. The diversity of resistant types as a result of the variability in intrinsic cell states may be a generic feature of responses to external cues.

Extended Data Fig. 1 ∣. FateMap reveals between-clone fate type diversity in treatment-naive cells, albeit to a lesser degree compared to resistant cells.

a. (left) UMAP of all barcoded treatment-naive cells. Total 16,432 cells (8,420 split A and 8,012 cells in split B) are colored by clusters determined using Seurat’s FindClusters command at a resolution of 0.6 (i.e. “Seurat clusters, resolution = 0.6”). (right) On the UMAP, we recolored each cell by its expression for a select subset of genes that were identified as differentially expressed in drug resistant cells via the Seurat pipeline (Cell counts available in Supplementary Table 11). b. Five representative examples demonstrate that a clone (cells sharing the same barcode) is constrained largely in a specific transcriptional cluster such that cells within a clone are more transcriptionally similar to each other than cells in other clones. c. Average pairwise correlation between cells within a clone was estimated based on the expression levels of the top 500 most variable genes. Each point represents the average value for Spearman’s correlation coefficient for all possible pairs of cells within a clone. For each clone, a paired control was created by randomly sampling an equivalent number of cells from the entire population. Higher average correlation coefficient in clones indicates higher transcriptional similarity among cells within a clone, as compared to cells that are not clones. Wilcoxon signed rank exact test (paired, two-sided) was used to compare the difference in average correlation coefficient. d. Fraction of variance explained by the experimental data and randomized data for the top 50 principal components (PCs). The number of PCs needed to explain the actual variance in data (indicated by the dotted line) is a measure of the degrees of freedom of variability of a given dataset. There was an increase (see Extended Data Fig. 1e for statistical testing) in the number of PCs needed to explain the variance in data from resistant cells (43 PCs) as compared to naive (30 PCs) and primed cells (23 PCs), suggesting that there is an increase in overall variability in samples when cells transition to becoming drug resistant. Primed cells were identified as cells where at least 40% of pre-resistant markers identified in (Emert et al. 2021) are higher than their average expression level. e. Average number of PCs needed to explain the variance in resistant, naive and pre-resistant cells. Error bars represent standard deviation over 100 simulations of randomized data. Mann-Whitney U-Test was used to estimate a p-value for pairwise difference in means. f. Comparison of Euclidean distances between clusters across resistant and naive populations of melanoma cells for varying numbers of clusters. We used the first 50 principal components to calculate the Euclidean distance between cells across clusters. We used Wilcoxon signed rank exact test (paired, two-sided) for statistical comparisons. g. Comparison between resistant and naive populations for total number of clusters, given fixed number of cells and shared nearest neighbor (snn) resolution. We used Wilcoxon signed rank exact test (paired, two-sided) for statistical comparisons of average number of clusters across resolutions. h. UMAPs of representative twin clones (sharing the same barcode) across the two splits A and B. The twins largely end up with the same transcriptional fate. This observation suggests that cells have similar transcriptional states prior to drug treatment.

Extended Data Fig. 2 ∣. Whole genome sequencing of treatment-naive and drug resistant fate type clones.

a. We performed a pairwise hypergeometric test for variants in all clones to determine statistical significance of variant overlap between clones. This was calculated with the following parameters (M = all CADD > 15 variants, n = # Variants in Sample 1, N = # Variants in Sample 2, X = # Variants in intersection of Samples 1 & 2). P-values are plotted on the heatmap where the p-value represents the probability of observing at least as large an overlap as observed if the two clones in fact had independently randomly selected variants from the full list of CADD > 15 variants. P-values below 0.05 represent two clones that are not genetically independent. b. Heatmap of genes with deleterious variants (CADD > 15) that were present with allele frequencies between 25% and 75% in both naive and resistant clones, colored by their CADD deleteriousness score. For genes that include multiple, unique variants, the variants were collapsed into one row, where the variant with the highest CADD score was plotted for that sample. The curated gene set represents the lack of variation in (Shaffer et al. 2017; Garman et al. 2017). Differentially expressed genes from the FateMap dataset show eight genes with variation. c. The expression patterns of the eight genes from the DEG list from FateMap with heterogeneously present genetic variants, visualized on UMAP (Cell counts available in Supplementary Table 11). d. To evaluate for acquired genetic resistance to therapy in the resistant clones, we next plotted variants on a heatmap (colored by their CADD deleteriousness score) if there was a significant difference in the allele frequencies of variants between naive and resistant clones by Fisher’s exact test (P < 0.05). Variants were only included if they were not present in any naive clones. The curated gene set represents the lack of acquired variants in genes from [2017 paper genes], [FateMap DEG], [Clone Genes], [Top 500 most Variable Genes], [Known Epigenetic Modifiers]. The [“all variants with CAAD > 15”] includes all variants with CADD c-scores over 15. e. We analyzed 143 genes classified as epigenetic modifiers for deleterious variants (CADD > 15) within naive clones. The chart shows the number of genes with variants in a subset of naive clones (2 genes) and in all naive clones (10 genes). f. Heatmap of deleterious variants in epigenetic modifier genes, colored by their CADD deleteriousness score.

Extended Data Fig. 3 ∣. Isolation, longitudinal profiling and functional mapping of drug resistant clones.

a. Schematic for longitudinal tracking and profiling of drug resistant colonies. Colonies were isolated, expanded and maintained over 4 to 6 weeks. Paired initial and late samples were sequenced at a bulk-level. b. Paired initial and late samples display minimal phenotypic drift in principal component (PC) space for top 500 most variable genes. Insets show brightfield images of representative samples. c. Euclidean distance (in PC1 and PC2) measured between paired initial and late samples and equivalent number of random initial-late pairs of samples. Lower Euclidean distance in true pairs as compared to random pairs implies that paired initial and late samples are transcriptionally more similar (closer in PC space) than any pair of initial and late samples. d. Scree plot depicting cumulative variance explained by PCs. Dotted line represents that most of the variance can be explained by the first 25 PCs alone. e. Euclidean distance (in first 25 PCs) measured between paired initial and late samples and equivalent number of random initial-late pairs of samples. Lower Euclidean distance in true pairs as compared to random pairs implies that paired initial and late samples are transcriptionally more similar (closer in PC space) than any pair of initial and late samples. f. Euclidean distance measured between paired early and late samples and equivalent number of random initial-late pairs of samples. Euclidean distance was measured in PC1 and PC2 space for top 200, 500 and 1000 variable genes. g. Euclidean distance measured between paired initial and late samples and equivalent number of random initial-late pairs of samples. Euclidean distance was measured in the PC space created by the first 25 PCs for top 200, 500 and 1000 variable genes. h. Mapping of invasiveness onto the single-cell RNA sequencing dataset from FateMap by comparing genes differentially expressed between the two slowest and the two fastest invading resistant colonies (UMAP colored for similarity score). The slowest invading colonies have a high similarity score for cluster 15 (and to some extent 4 and 6), while the fastest invading colonies have a high similarity score for cluster 8 (and to some extent 1).

Extended Data Fig. 4 ∣. FateMap on BRAF and NRAS mutant melanoma cell lines reveals between-clone fate type diversity.

a. (left) For another single-cell derived melanoma cell line WM983B E9-C6, we traced representative resistant cells in Adobe Illustrator and created cartoon schematics based on visual inspection of orientation and density. (right) Brightfield images of resistant colonies exhibiting different types of morphologies. b. We applied the Uniform Manifold Approximation and Projection (UMAP) algorithm within Seurat to the first 50 principal components to visualize differences in gene expression. Cells are colored by clusters determined using Seurat’s FindClusters command at a resolution of 0.5 (i.e. “Seurat clusters, resolution = 0.5”) (13,869 and 11,249 total cells respectively for split A and B). c. On the UMAP, we recolored each cell by its expression for a select subset of genes that were identified as differentially expressed via the Seurat pipeline and marked different clusters. MLANA, which marks melanocytes, is found largely in clusters 1,3, and 5; IFIT2, which marks type-1 interferon signaling, is found largely in cluster 4; NGFR, which marks neural crest cells, is found largely in cluster 2 and 4; AXL, which is a canonical resistance marker, is found largely in cluster 5 and 7. d. Six examples to demonstrate that a clone (cells sharing the same barcode) is constrained largely in a specific transcriptional cluster such that cells within a clone are more transcriptionally similar to each other than cells in other clones. Some clones are larger in size than others, and some exist as singletons, meaning they survive vemurafenib treatment but do not necessarily divide while exposed to the drug. e. We quantified the preference for a specific cluster across all barcode clones (clone size>4). Specifically, we calculated the fraction of dominant clusters for each clone and found it to be significantly higher (Wilcoxon test, two-sided, unpaired, p-value = 1.49e-15) than that for randomly selected cells. The analysis plotted here is for a cluster resolution of 0.5. f. Painting of singletons and colonies onto the UMAP demonstrated that singletons and colonies belonged to distinct regions and clusters. g. UMAPs of representative twin clones (sharing the same barcode) across the two splits A and B. The twins largely end up with the same transcriptional fate type, invariant of the clone size. This observation suggests that cells are predestined for distinct resistant fate types upon exposure to vemurafenib. h. NRAS mutant cell line WM3623 treated with three different doses of trametinib (10 nM, 20 nM, and 40 nM). Representative brightfield images after 2.5 and 5 weeks of drug treatment are shown for each dose. i. (left) UMAP of all barcoded 3623 cell line cells treated with Trametinib. 6,397 cells are colored by clusters determined using Seurat’s FindClusters command at a resolution of 0.6 (i.e. “Seurat clusters, resolution = 0.6”). (right) On the UMAP, we recolored each cell by its expression for a select subset of genes that were identified as differentially expressed in drug resistant cells via the Seurat pipeline. j. On the UMAP, we recolored each cell by its expression for a select subset of genes that were identified as differentially expressed in drug resistant cells via the Seurat pipeline. k. Five representative examples demonstrate that a clone (cells sharing the same barcode) is constrained largely in a specific transcriptional cluster such that cells within a clone are more transcriptionally similar to each other than cells in other clones. l. UMAPs of representative twin clones (sharing the same barcode) across the two splits A (6,397 cells) and B (7,538 cells). The twins largely end up with the same transcriptional fate type. This observation suggests that drug resistant cells are derived from the same clones having similar transcriptional states and are constrained in the gene expression space. One of the clones appears to be a dominant clone and gives rise to a large fraction of sequenced cells.

Extended Data Fig. 5 ∣. FateMap on an NRAS mutant melanoma cell line reveals between-clone fate type diversity.

a. NRAS mutant cell line WM3451 treated with three different doses of trametinib (20 nM, 40 nM, and 50 nM). Representative brightfield images after 2.5 and 5 weeks of drug treatment are shown for each dose. b. (left) UMAP of all barcoded 3451 cells treated with Trametinib. 5,789 cells are colored by clusters determined using Seurat’s FindClusters command at a resolution of 0.6 (i.e. “Seurat clusters, resolution = 0.6”). (right) On the UMAP, we recolored each cell by its expression for a select subset of genes that were identified as differentially expressed in drug resistant cells via the Seurat pipeline. c. Five representative examples demonstrate that a clone (cells sharing the same barcode) is constrained largely in a specific transcriptional cluster such that cells within a clone are more transcriptionally similar to each other than cells in other clones. d. Average pairwise correlation between cells within a clone was estimated based on the expression levels of the top 500 most variable genes. Each point represents the average value for Spearman’s correlation coefficient for all possible pairs of cells within a clone. For each clone, a paired control was created by randomly sampling an equivalent number of cells from the whole population. Higher average correlation coefficient in clones indicates higher transcriptional similarity among cells within a clone, as compared to cells that are not clones. Wilcoxon signed rank test (paired, two-sided) was used to compare the difference in average correlation coefficient. e. UMAP of all barcoded WM3451 P2G7 cells treated with Trametinib. Cells are colored by whether they are a singleton (i.e. clone size = 1). f. UMAPs of representative twin clones (sharing the same barcode) across the two splits A (5,789 cells) and B (7,473 cells). The twins largely end up with the same transcriptional fate type. This observation suggests that drug resistant cells are derived from the same clones having similar transcriptional states and are constrained in the gene expression space.

Extended Data Fig. 6 ∣. FateMap on treatment-naive primary human melanocytes reveals between-clone diversity.

a. (left) UMAP of all barcoded naive primary melanocyte cells. Cells are colored by clusters determined using Seurat’s FindClusters command (“Seurat clusters, resolution = 0.6”). (right) On the UMAP, we recolored each cell by its expression for a select subset of genes that were identified as differentially expressed in drug resistant cells via the Seurat pipeline. b. Six representative examples demonstrate that a clone is constrained largely in a specific transcriptional cluster such that cells within a clone are more transcriptionally similar to each other than cells in other clones. c. Average pairwise correlation between cells within a clone was estimated based on the expression levels of the top 500 most variable genes. Each point represents the average value for Spearman’s correlation coefficient for all possible pairs of cells within a clone. For each clone, a paired control was created by randomly sampling an equivalent number of cells from the whole population. Higher average correlation coefficient in clones indicates higher transcriptional similarity among cells within a clone, as compared to cells that are not clones. Wilcoxon signed rank test (paired, two-sided) was used to compare the difference in average correlation coefficient. d. UMAP of all barcoded naive primary melanocyte cells. 2,868 cells are colored by whether they are a singleton (i.e. clone size = 1). Cluster 10 is enriched for singletons and displays high expression of S100B, a marker identified to be associated with single cell colonies by FateMap. e. UMAPs of representative twin clones (sharing the same barcode) across the two splits A (2,868 cells) and B (3,333 cells). The twins largely end up with the same transcriptional fate type. This observation suggests that primary melanocyte cells derived from the same clone have similar transcriptional states and are constrained in the gene expression space.

Extended Data Fig. 7 ∣. FateMap on a triple negative breast cancer cell line reveals between-clone fate type diversity.

a. Nuclei scans (DAPI-stained) of resistant colonies emerging from treatment of the single-cell derived triple negative breast cancer cell line MDA-MB-231-D4 with 1nM paclitaxel. b. For the MDA-MB-231-D4 cell line, we traced representative resistant cells in Adobe Illustrator and created cartoon schematics based on visual inspection of orientation and density. c. Brightfield images of resistant colonies exhibiting different types of morphologies. d. We applied the Uniform Manifold Approximation and Projection (UMAP) algorithm within Seurat to the first 50 principal components to visualize differences in gene expression. 6,535 cells are colored by clusters determined using Seurat’s FindClusters command (“Seurat clusters, resolution = 0.5”). e. We observed silencing of the transcribed barcodes in a subset of colonies, as revealed by epifluorescence imaging of the GFP signal. The colony on the left is strongly expressing a GFP signal while the colony on the right has a very dim GFP signal. f. Cells with assigned barcodes were evenly distributed throughout the UMAP with no clear bias for any specific resistant fate types. g. Four examples from split A (6,535 cells) demonstrate that a clone is constrained largely in specific UMAP regions such that cells within a clone are more transcriptionally similar to each other than cells in other clones. h. Four examples from split B (8,745 cells) demonstrate that a clone is constrained largely in specific UMAP regions such that cells within a clone are more transcriptionally similar to each other than cells in other clones. i. We quantified the preference for a specific cluster across all barcode clones (clone size>4). Specifically, we calculated the fraction of dominant clusters for each clone and found it to be significantly higher (Wilcoxon, unpaired, two-sided) than that for randomly selected cells. The analysis plotted here is for a cluster resolution of 0.5. j. We found that our UMAP had superclusters defined by cell cycle (S, G1, G2M). Of the 3,720 clonal DEGs, 63 are cell cycle genes. We therefore regressed out cell cycle genes and the cell-cycle-genes-regressed data with UMAP. k. We quantified the preference for a specific cluster across all barcode clones after cell cycle regression (clone size>4). Specifically, we calculated the fraction of dominant clusters for each clone and found it to be significantly higher (Wilcoxon, unpaired, two-sided) than that for randomly selected cells. The analysis plotted here is for a cluster resolution of 0.5. l. UMAPs of representative twin clones across the two splits A and B. The twins largely end up with the same transcriptional fate type, invariant of the clone size. This observation suggests that cells are predestined for distinct resistant fate types upon exposure to chemotherapy drug paclitaxel.

Extended Data Fig. 8 ∣. FateMap reveals differences in clonal fate type outcomes between continuous and discontinuous therapy.

a. Schematic of the experimental design where we exposed single-cell-derived WM989 A6-G3 melanoma cells to continuous and discontinuous doses of targeted therapy drug vemurafenib. b. UMAP of all barcoded cells. 17,634 cells are colored by clusters determined using Seurat’s FindClusters command (“Seurat clusters, resolution = 0.6”). c. Pellet morphology for continuous (7,238 cells) and discontinuous (10,396 cells) treatment cells. Cells derived from discontinuous dosage have a larger and darker (more pigmented) pellet. This suggests that during discontinuous dosage, melanocytic cells (which are pigmented in nature) proliferate. d. On the UMAP, we recolored each cell by its expression for a select subset of genes that were identified as differentially expressed in drug resistant cells via the Seurat pipeline. e. UMAP of all barcoded cells. Cells are colored by type of dosage. f. UMAPs of representative twin clones (sharing the same barcode) that arise during discontinuous drug treatment. The twins largely end up with the same transcriptional fate type and have varying proliferative capacities. g. In discontinuous dosage, 68% of clones having high MLANA expression (log2 Expression > 2, in at least 50% of cells in a given clone) are proliferative (i.e. have clone size > 1). In continuous dosage, only 20% of clones having high MLANA expression are proliferative. h. (left) Total number of cells analyzed consisted of 60.5% discontinuous dosage samples and 39.5% continuous dosage samples. (right) The number of unique barcodes (i.e. resistant clones) displays a 3.6 fold increase in discontinuous dosage sample as compared to the continuous dosage sample. i. UMAPs of representative twin clones across the two splits of continuous and discontinuous dosing. Some twins end up in the similar transcriptional fate type while others tend to switch fate type.

Extended Data Fig. 9 ∣. Changing the therapy type to trametinib eliminates an additional resistant fate type present in the vemurafenib treatment.

a. UMAP where the resistant cells are colored by the associated therapy drug type, with dark blue representing vemurafenib (9,457 cells) and light blue representing trametinib (8,569 cells). Arrows represent UMAP regions that are present only in vemurafenib or trametinib. b. UMAP is split by each drug type, with colors representing clusters determined using Seurat’s FindClusters command(“Seurat clusters, resolution = 0.5”). Arrows represent UMAP regions that are present only in vemurafenib or trametinib. c. Painting of singletons and colonies onto the UMAP, colored by the condition, demonstrated that singletons largely belong to vemurafenib and are present predominantly in the MLANA-high cluster. Colonies are dispersed more across the UMAP with no particular region enriched for either condition except for the NGFR-high cluster. d. Imaging of nuclei (DAPI-stained) of resistant colonies emerging from treatment of WM989 A6-G3 cells to either vemurafenib or trametinib. The number of singletons in trametinib treated cells appear to be much less than those treated with vemurafenib, consistent with the sequencing data from FateMap. e. Quantification of the total number of colonies and singletons from each drug type of imaging data across biological replicates. This analysis demonstrated that while the total number of colonies are similar across the two drug types, there is a relative increase (~2.45-fold; n = 3 biological replicates) in the number of singletons in the case of vemurafenib. Error bars represent standard error of the mean. f. UMAPs are recolored for each cell by its expression for gene MLANA, which is a marker for cluster 3 relatively enriched in vemurafenib (as shown with arrows in A and B). g. A pie chart to demonstrate that of all the clones (barcodes) present in vemurafenib-treated split in cluster 3, only 4.8% were also present in the trametinib-treated split. h. A cumulative density contour plot capturing the types of fate switches that the MLANA-high cluster 3 clones from the vemurafenib-treated split adopt in the trametinib-treated split. i. Three representative examples of UMAP regions where twins from the MLANA-high cluster 3 in the vemurafenib-treated split adopt in the trametinib-treated split. j. UMAPs are recolored for each cell by its expression for gene NGFR, which is a marker for cluster 4 relatively enriched in trametinib (as shown with arrows in A and B). k. Composition of clones of different sizes within NGFR-high cluster 4 for both trametinib- and vemurafenib-treated splits. l. A pie chart to demonstrate that of all the clones (barcodes) present in the trametinib-treated split in cluster 4, 20.7% were also present in the vemurafenib-treated split. m. A cumulative density contour plot capturing the types of fate switches that the NGFR-high cluster 4 clones from the vemurafenib-treated split adopt in the trametinib-treated split. n. Two representative examples of UMAP regions where twins from the NGFR-high cluster 4 in trametinib-treated split adopt in the vemurafenib-treated split. o. UMAP for combined vemurafenib and trametinib treatment conditions recolored for each cell by its expression of the gene VCAM1, which is enriched in cluster 6. p. Painting of singletons and colonies onto the UMAP for the NGFR-high cluster 4, colored by the condition, showing a relative enrichment of cells from trametinib as compared to vemurafenib. This panel also demonstrates that both singletons and colonies occupy cluster 4 from each of the two conditions. q. We performed antibody stainings for NGFR on colonies emerging from treatment of the same number of starting melanoma cells with either vemurafenib or trametinib. Consistent with FateMap, we found an increased number of NGFR-positive resistant cells in trametinib treated cells as compared to the vemurafenib treatment. r. UMAP split by each drug condition (trametinib (8,569 cells) or vemurafenib and trametinib (7,023 cells)), with colors representing clusters determined using Seurat’s FindClusters command at a resolution of 0.5 (i.e. “Seurat clusters, resolution = 0.5”). s. UMAP recolored for combined resistant cells from trametinib (light blue) and vemurafenib and trametinib (dark blue). The cells from two conditions are interspersed into each other on the UMAP.

Extended Data Fig. 10 ∣. Inhibition of histone methyltransferase DOT1L results in the emergence of additional resistant proliferative clones and a reduction in singletons.

a. For each barcode identified by sequencing, we plotted its abundance in corresponding splits A (DMSO control) and B (DOT1L inhibition). Those present in both control and DOT1L splits are colored in dark blue, and those present only in either A (control) and B (DOT1L) are colored in cyan. Those present in both (dark blue; 171) exhibited a strong correlation, suggesting that their ability to survive and become resistant is invariant of drug dose. For those present only in either (cyan), we found them to be much more abundant in DOT1L (B, 43 barcodes) than DMSO control (A, 7 barcodes), suggesting that new barcodes, otherwise unable to survive in the control condition, become drug-resistant in the DOT1L inhibited condition. A total of one biological replicate. b. (left) Combined resistant cells in the control (9,343 cells) and DOT1L (7,044 cells) conditions obtained from UMAP applied to the first 50 principal components. Cells are colored by clusters determined using Seurat’s FindClusters command(“Seurat clusters, resolution = 0.8”). (right) UMAP is split by each condition. c. UMAP where the resistant cells are colored by the associated condition (control vs DOT1L). The arrow represents the UMAP region present predominantly in the control region and missing from the DOT1L-associated UMAP region. d. Quantification of singletons and colonies showed that while the number of resistant colonies is higher in DOT1L, it is accompanied by a reduced number of singletons cells compared to control. e. Painting of singletons and colonies onto the UMAP, colored by the condition, demonstrated that singletons largely belong to the control condition and are present predominantly in cluster 2 (MLANA-high). Colonies are dispersed more across the UMAP with no particular region enriched for either condition. f. Imaging of the nuclei (DAPI-stained) of resistant colonies emerging from vemurafenib treatment of WM989 A6-G3 cells, either for control or cells lacking DOT1L. g. Quantification of the total number of colonies and singletons from each fate type across n = 3 biological replicates demonstrated a relative increase (3.65-fold; n = 3 biological replicates) in total colonies and reduction in total singletons in the DOT1L and control conditions, respectively. Error bars represent standard error of the mean. h. UMAP is recolored for each cell by its expression for the gene MLANA, a marker for cluster 2, which is relatively enriched in control (as shown with an arrow). i. A pie chart to demonstrate that of all the clones (barcodes) present in the control condition split, only 3.1% were also present in the DOT1L inhibitor pretreatment split. j. Two representative examples of UMAP regions where twins from the MLANA-high cluster in the control condition go in the DOT1L condition. A cumulative density contour plot capturing the types of fate switches that MLANA-high cluster clones from control adopt in the DOT1L inhibitor-pretreated condition. k. A cumulative density contour plot capturing the types of fate switches that the MLANA-high cluster 2 clones from the control condition split adopt in the DOT1L inhibitor pretreatment split. l. Distribution of cells across clusters for control (top) and DOT1L inhibitor-pretreated (bottom) conditions for clone size>2.

Fig. 1:. FateMap reveals that between-clone fate type diversity arises from a single cell upon therapy treatment.

a, Schcmatic of single-cell-derived WM989 A6-G3 melanoma cells exposed to the targeted therapy drug vemurafenib, which formed resistant colonies in 3–4 weeks. Colonies were mixed together and single-cell sequenced. b, Uniform manifold approximation and projection (UMAP) applied to the first 50 principal components to visualize gene expression differences. The 8,212 cells are coloured by clusters determined using the FindClusters command: “Seurat clusters, resolution - 0.6*. $n=1$ of 2 biological replicates. c, Cells on the UMAP recoloured by the expression of a subset of differentially expressed genes. Genes with similar UMAP expression profiles are listed below each panel. ACTA2 is found largely in Seurat cluster 8; IFIT2 is found largely in cluster 12; VCAM1 is found largely in cluster 15; NGFR is found largely in cluster 7; and MLANA is found largely in clusters 0 and 3. d, Schematic of FateMap labelling of cells with unique DNA barcodes before vemurafenib exposure. WM989 A6-G3 cells were transduced with the FateMap barcode library at a multiplicity of infection (MOI) of approximately 0.15. WPRE, woodchuck hepatitis virus post-transcriptional regulatory element; EFS, elongation factor 1α short. e, Barcodcd cells from d were exposed to vemurafenib for 3–4 weeks and resultant colonies were analysed by scRNA-seq and barcode sequencing. f, Testing whether resistant cells sharing a barcode (a resistant clone) are more transcriptionally similar to each other than other clones. g, Clones, irrespective of their size, are largely constrained in specific clusters. h, Quantification of preference for specific clusters across all clones (clone size > 4; representative clones in yellow). Wilcoxon test (unpaired, two sided), P < 2.2 × 10⁻¹⁶. SNN, shared nearest neighbour. i, RNA FISH of genes marking resistant types. Consistent with FatcMap, we found resistant colonies that were selectively positive for each of the three markers tested, and others that were negative for all of these markers.

Fig. 2:. Differences in gene expression between clones correspond to differences in morphology, proliferation and invasiveness.

a, Classification of colonies as singletons, small colonies or large colonies. Clusters exhibited different proliferative capacities. b, Colony size distributions for each fate type. Unpaired, two-sided Mann–Whitney U test; intervals based on P value thresholds. ****P < 0.0001, ***P < 0.001, **P < 0.01, *P < 0.05. n = 1 of 2 biological replicates. c, Left, schematics based on visual inspection of morphology, orientation and density. Right, bright-field images of resistant colonies exhibiting different morphologies. d, Schematic of isolation and expansion of vemurafenib-resistant or trametinib-resistant colonies, a subset of which were then analysed by bulk RNA sequencing (RNA-seq), categorized for morphology and measured for invasiveness. e, Resistant cells were seeded at 3,000 cells per well, allowed to form spheroids over 96–120 h and then embedded in a collagen matrix. Red and cyan mark the invading and core boundary, respectively. RC, resistant colony. f, Invasiveness of resistant colonies emerging from treatment with trametinib was quantified by computing the ratio of the area enclosing the red and cyan boundaries as shown in e. Each dot represents one spheroid. g, Mapping of morphology onto the FateMap data by comparing genes differentially expressed from morphology to resistant colonies. Similarity score on UMAP represents the degree of overlap of differentially expressed genes between bulk-sequencing data and each cluster. Resistant colony fate type 2 maps predominantly to cluster 7, whereas fate type 3 maps to cluster 15. Singletons, as identified by imaging for the SOX10 gene, map to clusters 0 and 3.

Fig. 3:. Resistant fate types emerge following targeted therapy in patients as evidenced by spatial transcriptomic profiling.

a, Overview of all 29 punch biopsies from four patients that were sequenced using the GeoMx Digital Spatial Profiler (DSP) system for spatial transcriptomics. A total of 93 ROIs were selected for sequencing based on visual inspection of DNA (SYTO 13, blue), CD45 (red) and S100B (green) staining. Two patient samples were coupled with matched pre-treatment biopsies (marked untreated). b, Dot plots of counts per million (cpm)-transformed data for selected markers of resistant fate types as identified from in vitro FateMap. Each dot is a single ROI coloured by that region’s qualitative S100B staining level, faceted by the punch biopsy of each region. Specific punch biopsies highlighting nearby regions with different expression profiles are highlighted for two patients, both before and after treatment.

Fig. 4:. Cells are predestined for distinct resistant fates upon exposure to therapy.

a, Schematic of FateMap twin experimental designs. We transduced WM989 A6-G3 cells (MOI ≈ 0.15) with the FateMap barcode library. After 3–4 cell divisions, we sorted the barcoded population, split the cells (into A and B), treated each with vemurafenib and performed scRNA-seq and barcode sequencing on the colonies (Supplementary Table 11). b, Unique barcode abundance is identified by gDNA sequencing in splits A and B. Those present in both splits are dark blue (87), and those present in only one (14 each) are cyan. n = 1 of 2 biological replicates. c, Top, Venn diagram of the overlap between barcode clones present in both splits (dark blue) compared with those present only in either A or B (cyan). Bottom, comparison of the observed overlap between the shared barcodes (twins) surviving across splits with random survival chance (simulated 1,000 times). d, UMAPs of representative twin clones (sharing the same barcode) across the two splits A (8,212 cells) and B (7,262 cells). The resistant twins largely end up with the same transcriptional fate type, invariant of the clone size. e, Large clones superimposed on the UMAP, with each colour representing a unique resistant clone, f, Mixing coefficient is used to calculate the pairwise transcriptional relatedness of clones (see Methods). Higher mixing coefficient corresponds to higher transcriptional relatedness of clones (perfect mixing, 1; no mixing, 0). Representative example UMAPs are provided. g, Mixing coefficient for twin clones across splits A and B is presented with representative examples on the UMAP. h, Box plots showing cumulative mixing coefficients between clones within splits A (133) and B (102) (grey), non-twin clones across A and B (66) (grey), and twin clones (12) (blue). Unpaired, two-sided Wilcoxon test; P value for non-twin compared with twin clones is 4.513 × 10⁻⁸.

Fig. 5:. Changing the therapeutic dose results in stereotypic resistant fate type switching and altered transcriptional profiles.

a, Left, resistant colonies emerging from treatment of 25,000 WM989 A6-G3 cells with two vemurafenib doses (1 μM and 100 nM). Right, total colonies from each dose across n = 3 biological replicates; error bars represent s.e.m. b, Schematic of FateMap twin experimental designs for different vemurafenib doses. We transduced WM989 A6-G3 cells with the barcode library. After 3–4 cell divisions, we sorted the barcoded population, divided it into splits A and B, treated each with vemurafenib and performed FateMap. We list three possible outcomes (cell counts are available in Supplementary Table 11). c, Combined low-dose (13,400 cells) and high-dose (9,457 cells) resistant cells obtained from UMAP applied to the first 50 principal components. Cells are coloured by clusters determined using Seurat’s FindClusters command. d, UMAP with resistant cells coloured by dose. Arrows represent regions present at only one of the doses. e, The UMAP in d is split accoring to dose, with colours representing clusters determined in c. Arrows represent clusters that are present at only one dose. f, UMAPs recoloured according to expression of NGFR and MLANA, markers for clusters that are enriched in one of the two doses. g, Left, UMAP cluster coloured for cluster 9 (high for NGFR). Right, pie chart showing that 25.3% of the NGFR-high clones present in high dose were also detected in the low dose. h, Left and centre, representative examples of where twins from the NGFR-high cluster following high dose are located following low dose. Fate type switch 1 twins (38 out of 46) had similar fate types, as did fate type switch 2 twins (8 out of 46). Right, cumulative density contour plot of fate type switches from high dose to low dose.