Clonal differences underlie variable responses to sequential and prolonged treatment

Abstract

Cancer cells exhibit dramatic differences in gene expression at the single-cell level, which can predict whether they become resistant to treatment. Treatment perpetuates this heterogeneity, resulting in a diversity of cell states among resistant clones. However, it remains unclear whether these differences lead to distinct responses when another treatment is applied or the same treatment is continued. In this study, we combined single-cell RNA sequencing with barcoding to track resistant clones through prolonged and sequential treatments. We found that cells within the same clone have similar gene expression states after multiple rounds of treatment. Moreover, we demonstrated that individual clones have distinct and differing fates, including growth, survival, or death, when subjected to a second treatment or when the first treatment is continued. By identifying gene expression states that predict clone survival, this work provides a foundation for selecting optimal therapies that target the most aggressive resistant clones within a tumor. A record of this paper's transparent peer review process is included in the supplemental information.

Keywords: barcoding; cancer systems biology; clonal tracing; drug resistance; scRNA-seq.

Figure 1:. Melanoma cell clones exhibit variable responses to a first and second treatment.

A) Shown are two possible scenarios for how clones could respond to sequential treatment with different agents. All resistant clones derived from one treatment could have the same response to a second treatment (top), or they could have different responses to a second treatment (bottom). B) WM989 BRAF V600E mutant melanoma cells with a nuclear GFP tag imaged after four weeks in CoCl₂ followed by four weeks in dabrafenib and trametinib (Dab/Tram, top left). The top left images show whole well scans of the clones. The colored boxes show zoomed images selected from the whole well. Orange shows a clone that grew in the second treatment, blue shows a clone that remained a similar size in the second treatment, and red shows a clone that died in the second treatment. CoCl₂ resistant cells tend to grow on top of each other causing some cropped images to look out of focus. Scale bars in whole well scans are 5 mm while scale bars in cropped scans are 500 μm.

Figure 2:. DNA barcoding allows for tracking the transcriptome of clones through multiple rounds of treatment.

A) Schematic showing experimental protocol for clone tracking and scRNA-seq of WM989 BRAF V600E mutant melanoma cells through treatment with with dabrafenib and trametinib (Dab/Tram), CoCl₂, or cisplatin (see Methods). B) Schematic showing how sequencing reads from the genomic DNA of cells that survived a treatment provide a quantitative readout of treatment response. C) Heatmaps showing 100 randomly subsampled and rank-ordered clones detected after initial treatment with Dab/Tram (left), CoCl₂ (middle), and cisplatin (right) followed by their abundance in the second round of treatment. Individual clones are colored by the ln(Reads Per Million (RPM) + 1) using data from sequencing clonal barcodes from gDNA of surviving cells. D) The number of total detected clones per condition are displayed in their assigned color. It should be noted that there are clones whose abundance were below the threshold of detection after treatment 1 that were later detected in treatment 2. The total number of unique clones that were detected after the first or second treatment are displayed below each heatmap.

Figure 3:. The treatment history of cells are reflected in their gene expression.

A) Schematic showing the theoretical possibilities of the first treatment driving endstate gene expression (top) or the second treatment driving endstate gene expression (bottom). Displayed are how these possibilities would look in UMAP space or plot as pairwise Pearson correlation between cell transcriptomes. B) (Left) UMAP (GSE253739) of the nine conditions that had received both a first and second treatment. (Center) UMAP with cells that received the indicated treatment first highlighted in red and those that received the treatment second highlighted in blue. (Right) Quantification of cell similarity using pairwise Pearson correlations of gene expression. For combination dabrafenib and trametinib treatment (Dab/Tram), pairwise Pearson correlations between cells that received Dab/Tram first (treatment 1, T1) were ~3.5-fold higher than the matched random sampling of cells (random control 1, R1), while correlations between cells that received Dab/Tram second (treatment 2, T2) were ~37-fold higher than the matched control (random control 2, R2) and ~10.5-fold higher than those that received Dab/Tram first. For cisplatin treatment, pairwise Pearson correlations between cells that received cisplatin first were ~6-fold higher than the matched random sampling of cells, while correlations between cells that received cisplatin second were ~28.5-fold higher than the matched control and ~5-fold higher than those that received cisplatin first. For CoCl₂ treatment, pairwise Pearson correlations between cells that received CoCl₂ first were ~13.5-fold higher than the matched random sampling of cells and ~2.5 fold higher than those that received CoCl₂ second, while correlations between cells that received CoCl₂ second was ~5.5-fold greater than the matched random control. All comparisons were statistically significant by a two-sided t-test with p < 2.2e-16. Violin plots display 10,000,000 subsampled data points per condition, but statistical comparisons and averages were calculated on non-subsampled data. Mean values are displayed below each graph.

Figure 4:. Cells within a clone maintain gene expression similarities through multiple rounds of treatment.

A) (Left) Location of the cells from a real clone (purple) and its matched random control (tan) from one simulation overlaid on the dabrafenib and trametinib (Dab/Tram) resistant UMAP (GSE253739). (Middle) Pairwise Pearson correlations of gene expression (top 2000 most variable genes) of cells from the real clone. The average of these pairwise Pearson correlations was 0.43. (Right) Pairwise Pearson correlations of gene expression of cells from the simulated clone. The average of these pairwise Pearson correlations was −0.03. For the pairwise Pearson correlation figures, the blocks above the diagonal are shaded according to the linear color scale to indicate the correlation while the numerical value is displayed in the associated block in the lower diagonal. The averages of the pairwise Pearson correlations of the cells within each clone that received a treatment were compared to those of random samplings of the same number and sized clones over 100 simulations. Clones analyzed ranged from 5 to 11,257 cells. Displayed are the first simulations for each condition. As follows are the number of simulations where the observed similarities were higher than random for each condition by a one-sided Wilcoxon Rank Sum Test with p < 0.05: B) Dab/Tram 100/100, C) Dab/Tram to Dab/Tram 100/100, D) Dab/Tram to CoCl₂ 100/100, E) Dab/Tram to cisplatin 100/100, F) CoCl₂ 99/100, G) CoCl₂ to Dab/Tram 100/100, H) CoCl₂ to CoCl₂ 100/100, I) CoCl₂ to cisplatin 100/100, J) cisplatin 100/100, K) cisplatin to Dab/Tram 100/100, L) cisplatin to CoCl₂ 100/100, M) cisplatin to cisplatin 100/100.

Figure 5:. Treatment with dabrafenib and trametinib induces resistance to CoCl₂ in clones marked by high IL6ST expression.

A) Schematic showing how clone 2 only survived CoCl₂ treatment after being treated with dabrafenib and trametinib (Dab/Tram). B) UMAP (GSE253739) of dabrafenib and trametinib resistant cells. Cells from clones that will not survive (sensitive to) CoCl₂ are colored blue, and cells from clones that were induced to survive CoCl₂ are colored in red (these did not survive CoCl₂ when applied as the first treatment). C) UMAP of dabrafenib and trametinib resistant cells showing relative expression of IL6ST. The color scale describes log-normalized data. D) Schematic for validation of induced resistance. Cells were treated with dabrafenib and trametinib for four weeks. Dabrafenib and trametinib resistant cells were then stained for expression of IL6ST protein and sorted into high (top 30%) and low (bottom 30%) populations. These IL6ST-high and -low populations were treated with CoCl₂ for an additional four weeks and then survival and growth was assessed by imaging. E) Representative DAPI images of IL6ST-low and -high dabrafenib and trametinib resistant cells after CoCl₂ treatment (one out of six analyzed whole-well images are displayed for each of the IL6ST-Low and IL6ST-High conditions). We binarized the DAPI stained nuclei images such that nuclei were white and the remaining pixels black. Scale bars are 2 mm. F) Violin plots showing the number of white pixels per well in IL6ST-low vs -high cells. Each data point is a separate well from the experiment. These numbers were normalized to a plate that was scanned the day after playing to eliminate differences from seeding. p < 5×10^−4 determined by a one-sided t-test where the IL6ST-high cells were expected to have greater survival than the IL6ST-low cells. A second biological replicate is shown in Supp. Fig. 6D.

Figure 6:. Prolonged treatment with the same agent causes population-level and clonal changes in pathway expression.

A) UMAPs (GSE253739) comparing cells sequenced after four and 6.5–8 weeks in treatment with (left) dabrafenib and trametinib (Dab/Tram), (middle) cisplatin, or (right) CoCl₂. B) We identified differentially expressed genes between 6.5–8 and four weeks of treatment with each agent. We then identified GO terms enriched after either the first round of treatment (blue) or after the second round of treatment (red) based on the log₁₀(q-value) and -log₁₀(q-value) respectively and displayed these values in the heatmap. This is a measure of the confidence of enrichment, not an activity score. GO terms that were not significantly enriched, or failed other thresholds, in either direction for a condition are colored in grey. Enriched pathways were grouped into categories of related function, and the pathways that make up each group are detailed in Supplemental Figure 7. Subpanels C-E perform perform the same analysis as performed in B with all of the cells from each condition (with that comparison displayed again) on a clone-by-clone basis for all clones with at least five cells in both rounds of treatment for Dab/Tram, cisplatin and CoCl₂, respectively. On the right of each subpanel is displayed the percentage of total sequenced single-cells per condition from each clone after each round of treatment. Pathways enriched in both directions are colored yellow.

Figure 7:. Clones with high EGFR expression exhibit enhanced survival when subjected to extended treatment with dabrafenib and trametinib.

A) Stacked colored bars represent the proportion of cells from each clone with at least five cells after initial treatment (four weeks) with dabrafenib and trametinib (Dab/Tram) in the (magenta) EGFR-high cluster or (green) NGFR-high cluster (top). Dotted lines separate clones that are all EGFR-high, mixed EGFR and NGFR expression, and NGFR-high. Aligned below are the same clones after the second round of treatment (6.5 weeks) with Dab/Tram. Clones that still have greater than five cells are again colored by their proportion of cells in either the EGFR- or NGFR-high cluster while clones that are no longer detected (less than five cells) are colored white. B) Experimental protocol for comparing the sizes of EGFR-high and NGFR-high clones. C) (Left) Example image (subset of replicate 2, from Supp. Fig. 11B) showing dabrafenib and trametinib resistant cells labeling nuclei with DAPI (grey) and labeling EGFR (magenta spots) and NGFR (green spots) using Hybridization Chain Reaction (HCR). Predominantly EGFR-high clones are circled in magenta, predominantly NGFR-high clones are circled in green, and mixed clones are circled in yellow. (Right) Cropped images highlighting that HCR labels individual transcripts with spots for EGFR and NGFR. Scale bars are 200 μm. D) Quantification of replicate 1 clone sizes based on total DAPI fluorescence within circled clones from Supp. Fig. 11A. Data is rank-ordered based on clone size and clones are colored based on their classification: EGFR-high, NGFR-high, mixed, or ambiguous. Quantification of replicates 2 and 3 are in Supp. Fig. 11D–E.