Longitudinal single-cell RNA sequencing of patient-derived primary cells reveals drug-induced infidelity in stem cell hierarchy

Abstract

Chemo-resistance is one of the major causes of cancer-related deaths. Here we used single-cell transcriptomics to investigate divergent modes of chemo-resistance in tumor cells. We observed that higher degree of phenotypic intra-tumor heterogeneity (ITH) favors selection of pre-existing drug-resistant cells, whereas phenotypically homogeneous cells engage covert epigenetic mechanisms to trans-differentiate under drug-selection. This adaptation was driven by selection-induced gain of H3K27ac marks on bivalently poised resistance-associated chromatin, and therefore not expressed in the treatment-naïve setting. Mechanistic interrogation of this phenomenon revealed that drug-induced adaptation was acquired upon the loss of stem factor SOX2, and a concomitant gain of SOX9. Strikingly we observed an enrichment of SOX9 at drug-induced H3K27ac sites, suggesting that tumor evolution could be driven by stem cell-switch-mediated epigenetic plasticity. Importantly, JQ1 mediated inhibition of BRD4 could reverse drug-induced adaptation. These results provide mechanistic insights into the modes of therapy-induced cellular plasticity and underscore the use of epigenetic inhibitors in targeting tumor evolution.

Fig. 1

Modelling of drug-induced tumor evolution in vitro, in vivo and “in patients”. a, b Cisplatin treatment of PDPCs that are representative of phenotypically heterogeneous (HN137Pri) or homogeneous (HN120Pri) populations. Drug-induced selection of ECAD+ epithelial cells (green), and the elimination of Vim+ mesenchymal (red) cells, in HN137Pri model (a). De novo emergence of Vim + mesenchymal (red) cells within the phenotypically ECAD + epithelial (green) HN120Pri parental population (b) (n = 24 replicates per experimental condition). c, d Immunofluorescence-based characterization of naive and drug-resistant cells from HN120 and HN137 populations. Cells were stained for the expression of epithelial ECAD (in green) and mesenchymal VIM (in red) markers (c), as well as basal KRT5 (in green) and luminal KRT18 (in red) markers (d). e, f In vivo modelling of cisplatin-resistance in subcutaneous Pdx models of HN120Pri and HN137Pri cells (also see Supplementary Fig. 3a-g). Flow-cytometeric analysis of ECAD and VIMENTIN positive (VIM+) cells reveals a reduction of ECAD + population and the de novo gain of VIM + cells in the HN120 model (e). Loss of VIM + cells and an enrichment of ECAD + population is observed post-cisplatin selection in the HN137 model (f). The gating strategy for viable singlets is exemplified in Supplementary Figure 3a. g “In patient” validation of epithelial (ECAD+) to mesenchymal (VIM+) cell-state switch in HN120Met vs. HN120Pri patient tumor, and the selection of ECAD + epithelial cells with concomitant loss of VIM + mesenchymal cells in cisplatin treated HN137 patient recurrent (and thus resistant) tumor in the clinic. n = 24 replicates (a, b) and n = 3 replicates (c–f) per experimental condition). Scale bar = 100 μM (a, b, g), 50 μM (c–d)

Fig. 2

Single cell RNA-seq of naive, drug-resistant and drug-holiday cells from OSCCs. a Schematic representation of the generation of drug-resistant models and scRNA-seq workflow. b RaceID for Hierarchical k-means clustering of 1302 scRNA-seq libraries based on gene expression profiling identifies five major clusters. c t-SNE visualization of clusters identified from OSCC-PDPC scRNA-seq libraries. d Color-coding of t-SNE plot based on the identities of individual OSCC-PDPCs and their drug-resistant/holiday models. e Relative expression of genes associated with basal (KRT5/KRT14), luminal (KRT8/KRT18), epithelial (EPCAM/Ecad), and mesenchymal (VIM/SPARC) markers on tSNE plot suggesting that the clusters represent distinch phenotypic cell-states

Fig. 3

Divergent developmental trajectories of drug-resistance in OSCCs. a, b Distribution of naive primary cells and their drug-resistant/holiday models on tSNE plots for HN120 (a) and HN137 (b). Note the switching of HN120 cells between distinct phenotypic clusters as they evolve from epithelial HN120Pri (cyan) to mesenchymal resistant HN120PCR (purple) and drug-holiday HN120PCRDH (dark-grey) cell states (a). Also note the heterogeneity within the HN137Pri population (spread of ‘red’ cell between two clusters), and the selective retention of the epithelial cell state in the HN137PCR resistant (dark green) and HN137PCRDH holiday model (dark red) (b). c Visualization of HN120, and HN137 naive metastatic cells, as well as their drug-resistant counterparts on tSNE. d Clustering of HN120 primary and metastatic drug-holiday model with intrinsically resistant HN148 cells. e PCA (basal/luminal genes) based hierarchical clustering of OSCC-PDPCs and their drug-resistant/holiday models with PAGODA, note the retention of basal/epithelial properties in HN137PCR, while gain of mesenchymal gene-signature in HN120PCR

Fig. 4

Drug-induced gain or loss of SOX2 determines the phenotypic identify of resistant cells. a Clustering of scRNA-seq libraries based on the expression of stemness genes followed by projection on original RaceID (depicted in Fig. 2d) tSNE plot. Cluster containing epithelial-like HN120Pri and HN137Pri cells (based on Fig. 2e) is marked with a red circle. Note the four sub-clusters within the epithelial-like cells based on the expression of stem cell markers. b Expression of SCC stem-like cell gene, SOX2, in a cluster consisting of epithelial-like HN120Pri and HN137Pri cells. c Beeswarm plot of SOX2 expression in different cell types from 1302 scRNA-seq libraries. d Immunofluorescence-based analysis of SOX2 (green) expression in HN120 and HN137 primary, drug-resistant, and metastatic models. Insets representing the magnified sections from respective immune-micrographs. e Validation of in vitro PDPC models in matched patients. Existence of SOX2 (green) and SOX9 (red) positive cells in HN120 primary patient tumor. Note the gain of SOX2+ (green) subpopulation in cisplatin treated HN137 patient tumor. Green arrows depict the sub-population of bona fide SOX2+ cells in primary tumors. f Limiting dilution assays (LDA) to assess the tumor-initiating potential of SOX2 expressing (HN120Pri, HN137Pri, and HN137PCR) vs. non-expressing HN120PCR cells (n = 4 mice per limiting dilution). g Proposed model for the stem cell-state switch from SOX2→SOX9 in HN120PCR and retention of SOX2 expressing phenotypic state in HN137PCR. h, i RaceID2 projection of the lineage tree of sensitive and resistant lines. HN120Pri and their cisplatin-resistant HN120PCR cells form two distinct lineage trees (h), while a common lineage-tree was observed in HN137Pri and their drug-resistant HN137PCR model (i). Scale bar = 50 μM (d, e)

Fig. 5

Pioneer factor SOX9 drives drug-induced cellular reprogramming and marks the acquired mesenchymal cell state. a Identification of the master regulators of mesenchymal cell fate, as analyzed by Vimentin expression, by an siRNA screen for ~2000 transcription factors in HN120Met. SOX9 was identified as a top candidate hit, the knockdown of which resulted in a marked loss of Vimentin expression. b Beeswarm plot of SOX9 expression in different cell types from 1302 scRNA-seq libraries. c, d Dependency of SOX2 and SOX9 on cisplatin resistance in HN120PCR/HN137PCR models. Graphs depict dose-dependent viability of resistant models treated with cisplatin. Loss-of SOX9 but not SOX2 reversed cisplatin-resistance in HN120PCR (c), while siRNA-mediated knockdown of SOX2, but not SOX9 conferred a significant loss of cisplatin-resistance in HN137PCR cells (n = 3, mean + s.e.m. of biological and technical replicates). e Sphere-forming assays with 500 single cells from HN120Pri, HN120PCR, HN137Pri, and HN137PCR cells. Shown are the effect of genetic loss/gain-of-function of SOX2 and SOX9 on sphere forming ability of treatment naïve and resistant cells (n = 3). f Immunofluorescence-based expression analysis of lineage markers (ECAD in green and Vimentin in red) and stem cell factors (SOX2 in green and SOX9 in red) in FFPE-section of matched OSCC patient tumors: treatment naive and locally recurrent. Inset demonstrating the de novo emergence of Vimentin expressing cells; note the SOX2 to SOX9 switch in these tumors. g Virtual sorting of TCGA-HNSCC in SOX2^High/SOX9^Low, SOX2^High/SOX9^High, SOX2^Low/SOX9^High, and SOX2^Low/SOX9^Low sub-sets. h Differential expression of EMT markers (in red rectangle) in sub-sets based on the “sorted” SOX2 and SOX9 populations. i Five-year survival analysis of the SOX2/SOX9 subsets identified in (g and h). Note patient with SOX2^Low/SOX9^High signatures display poorer survival (Cox regression P value < 0.05). Scale bar = 1000 μM (e), 50 μM (f)

Fig. 6

Poised mesenchymal promoters in treatment naive epithelial cells. a Chip-seq analysis for H3K4me3 in HN120 naïve, resistant and drug-holiday models. H3K4me3 marks on promoters of cisplatin-induced genes associated with cellular-reprogramming (EMT) in HN120 cells. Clusters CP-1 and CP-2 represent poised promoters, whereas CP-3 is associated with closed promoters. b Median expression of cisplatin-dependent upregulated genes from HN120 cells in clusters 1–3 (from Fig. 6a). c Quantitation of H3K27ac signals on promoters of cisplatin-dependent up-regulated genes in clusters 1–3 (from Fig. 6a) b, c mean + s.d., Mann–Whitney U test, ***P < 0.0005, *P < 0.005, #non-significant). d ChIP-seq tracks of H3K4me3 and H3K27ac on VIM, IL6, and GAS6 promoters in HN120Pri naive and drug-resistant/holiday models. Note the gain of H3K27ac marks in PCR/PCRDH cells (green asterisks), compared to naive (red asterisks), on EMT-associated promoters (y-axis = 0–75 (VIM, GAS6), 0–25 (IL6)). e Venn-diagram representing differential H3K27ac peaks in drug-sensitive and drug-resistant HN120 cells. HOMER based discovery of de novo transcription factor motifs in differential H3K27ac peaks. Note the loss of SOX2 motif in the drug-resistant HN120PCR cells, and the enrichment for SOX9 binding sites in the de novo acetylated promoters/enhancers. f Proposed model for SOX2 and SOX9 interplay during stem cell switch. Our data suggests that SOX2 expression maintains the epithelia cell-state while the loss of SOX2 and a concomitant gain of SOX9 results in the activation of an EMT program in HN120 drug-resistant and metastatic tumor cells

Fig. 7

Different modes of drug-resistance uncover potential for combinatorial targeting strategies using epigenetic inhibitors. a Synthetic-lethal screens with cisplatin using siRNAs against chromatin remodelers (see Fig. S6 for complementary screen using small molecule inhibitors). IC50 curves for cell viability screens in naive (b) and drug-resistant (c) HN120 cells with increasing concentration of either cisplatin or JQ1 alone, or a combination of JQ1 with 1 µM cisplatin (mean + s.e.m. n = 3, Two-tail Student’s t test, ***P < 0.0001, *P < 0.001). d, e Proposed model for cellular-reprogramming (adaptive) and ‘clonal selection’ modes of drug-resistance and metastasis. d The phenotypic homogeneous (HN120) model depicts the adaptive mode of evolution where epigenetically poised epithelial cell populations undergo mesenchymal cell-state transition. Different stress (cisplatin or metastasis) conditions leads to similar phenotypic outcome as a result of a SOX9-mediated coordinated activation of H3K4me3- poised promoters, which can be reversed by JQ1 mediated BRD4 inhibition. e Tumor cells exhibiting pre-existing phenotypic heterogeneity (HN137) allows “division-of-labor” where distinct clonal populations of pre-existing ‘designated survivors’ can get selected depending on the context specific nature of different selection pressures