Adaptive chromatin remodeling drives glioblastoma stem cell plasticity and drug tolerance

Sustained Inhibition of Kinase Signaling Enriches for Slow-Cycling GSCs

GSCs maintained in serum-free neurosphere culture conditions share features with neural stem cells (NSCs) and effectively initiate tumors in xenotransplantation assays (Lathia et al., 2015). To investigate proliferative programs in GSCs, we compared single-cell transcriptomes of primary GBM cells to in vitro GSCs (Patel et al., 2014). In primary tumors, only a small fraction of cells displays proliferative markers (2–20% Ki67⁺) (Louis et al., 2016) or express cell cycle signatures (Patel et al., 2014). When we compared developmental and cell cycle signatures, we found only a fraction of stem-like GBM tumor cells display proliferative signatures. In contrast, such signatures are evident in a large majority of in vitro GSCs (Figure 1A). While different GSC lines exhibit variable proliferation (Figure S1A) (Wakimoto et al., 2009), this potentially represents a critical distinction between in vitro and in vivo models.

We therefore considered if proliferative in vitro GSCs could be induced to a non-proliferative state. We treated a panel of GSC lines with small molecule inhibitors of oncogenic signaling pathways (Figure 1B, S1A–S1B, and Table S1). GSC8 and CW1691, both of which harbor focal amplification of PDGFRA, are highly sensitive to dasatinib (IC₅₀ ~ 10 nM, Figure 1B), a PDGFR/Src inhibitor. By contrast, four other models derived from PDGFRA-amplified tumors – GSC87, GSC114, GSC125, CW2907 – are largely insensitive to dasatinib.

In agreement with PDGFRα being the relevant target of dasatinib in the sensitive models, GSC8 and CW1691 are also sensitive to crenolanib, a PDGFR inhibitor that does not appreciably inhibit Src (Figure S1C). Furthermore, overexpression of a drug-resistant PDGFRA mutant (D842V) (Dewaele et al., 2008) partially rescues GSC proliferation in the presence of dasatinib or imatinib, another PDGFR inhibitor (Figure S1D–S1E), while shRNA knockdown of PDGFRA decreased GSC proliferation (Figure S1F). Dasatinib treatment reduced phosphorylation of PDGFRα (pY849), Akt (pT308) and Erk1/2 (pT202/pY204, pT185/pY187) (Figure 1C and S1G). While off-target effects of dasatinib cannot be ruled out, these findings support RTK inhibition as the predominant effector of drug treatment in these models. GSC8 is also sensitive to the MEK inhibitor PD0325901 (Figure S1B), suggesting that the effects of PDGFRa inhibition may operate via MAPK signaling.

Despite the effects of dasatinib, sustained exposure consistently yielded a persistent subpopulation that tolerates higher drug concentrations (Figure 1B). Specifically, treatment of naïve GSC8 (1 µM dasatinib, ~100 times IC₅₀) leads to an acute drop in the fraction of cycling, viable cells followed by a slow recovery (Figure 1D and Figure S2A–S2D). EdU pulse labeling followed by 6-day dasatinib treatment confirmed that initially cycling cells (EdU⁺) retained label and transitioned to a Ki67⁻ state (Figure 1E and Figure S2E). Subsequent washout of dasatinib resulted in rapid loss of EdU⁺ cells. These data suggest that while many cells acutely die upon exposure to dasatinib, a subset of proliferative cells undergoes cell cycle arrest, but remain viable and slowly expand despite the presence of drug. We refer to these drug persistent cells as GSC8^Per. Notably, GSC8^Per lines could also be derived with inhibitors of downstream kinases, including MEK and CDK4/6 (Figure S1B, S1H, and S2F). We also derived persisters from an EGFR-dependent line, GSC26, that tolerate higher concentrations of EGFR inhibitors (Figure S1I).

We isolated populations of GSC8^Per cultured in dasatinib for at least 8 weeks. These cells are relatively insensitive to dasatinib and crenolanib (Figure 1F), and still lack PDGFRα pY849 (Figure 1C), suggesting no mutation of PDGFRα conferring resistance. We considered if other genetic mechanisms might underlie drug-tolerance in GSC8 persisters. It has been reported that drug-tolerance to EGFR inhibitors may be mediated by reduction of EGFR⁺ extrachromosomal DNA and EGFR expression levels (Nathanson et al., 2014). In GSC8^Per, however, total levels of PDGFRα remained relatively constant or increased (Figure S1G and S1J), and low-coverage whole genome sequencing revealed no significant copy number changes at PDGFRA or other loci (Figure S1K). Most critically, the GSC8^Per state is reversible. Withdrawal of dasatinib permits full recovery of growth (Figure 1D), as well as re-sensitization to acute drug-induced arrest even after >4 months of chronic dasatinib treatment (Figure 1G). The rapidity and reversibility of acquired resistance support an epigenetic rather than genetic tolerance mechanism, as demonstrated in several cancer cell lines (Easwaran et al., 2014; Knoechel et al., 2014; Koppikar et al., 2013; Sharma et al., 2010).

Persister GSCs Express Primitive Neurodevelopmental and Quiescence Signatures

To investigate the regulatory circuits that sustain persister GSCs, we profiled gene expression in GSC8^Per, in naïve GSCs, and in short-term treated GSCs (GSC8^12d) at maximal arrest (Figure 1H). Clustering of the 4,084 genes with highest variability across these states revealed gene sets with coherent expression changes (Figure 2A, Table S2). Genes depleted in GSC8^12d (clusters 1–3) and GSC8^Per (clusters 1 and 3) were related to cell cycle and proliferation, consistent with reduced proliferation (Figure 2A, 2C, and S3A). By contrast, genes upregulated in GSC persisters (clusters 4–6) were enriched for signatures derived from quiescent NSCs (qNSCs) (Codega et al., 2014; Llorens-Bobadilla et al., 2015; Martynoga et al., 2013) and quiescent, stem-like medulloblastoma cells (Figure 2C and Figure S3B) (Vanner et al., 2014). Profiling of GSC26 and CW1691 treated with gefitinib or dasatinib, respectively, revealed similar changes in expression for the respective GSC8-derived gene clusters (Figure 2B), and enrichment of quiescent gene signatures (Figure 2C and Figure S3C).

Gene sets upregulated in persisters contain stemness factors previously linked to GSCs (Gangemi et al., 2009; Ikushima et al., 2009; Mehta et al., 2011; Suvà et al., 2014), as well as regulators implicated in stem-like tumor cells in vivo (Patel et al., 2014) (Figure 2A, clusters 4–6, Figure 2D and S3D). GSC8^Per also upregulate classical stemness markers, such as PROM1 (CD133) and SSEA-1 (CD15) (Singh et al., 2004; Son et al., 2009) (Figure 2A, 2E, and S3E). Many of these TFs and markers are already expressed in naïve GSCs but are upregulated on drug treatment (Figure 2D–2E and S3D). Notably, the persisters also upregulate multiple KDMs (Figure 2A, 2D, and S3F–S3I), including KDM3, KDM5, and KDM6 family members, some of which have been previously implicated in slow-cycling drug-resistant cell lines (Roesch et al., 2013; Sharma et al., 2010) and single-cell studies of primary GBM (Patel et al., 2014). Altogether, these expression changes are consistent with the notion that GSC persisters transition to a slow-cycling state and upregulate transcriptional programs shared with qNSCs.

Notch Signaling Underlies GSC^Per Switch

The strong upregulation of prominent Notch signaling genes prompted us to investigate this pathway (Figure 2A and 2D). Various notch pathway genes were upregulated in GSC8^Per, including Notch ligands, receptors, and canonical downstream targets (Figure 3A and Figure S4A–S4B). Consistent with these changes, GSC8 persisters have high levels of Notch1 intracellular domain (N1ICD) (Figure 3B), the proteolytic cleavage product of activated Notch1 receptor that mediates gene induction. The dasatinib-insensitive PDGFRA-amplified models, GSC87 and CW2907, already express Notch pathway genes at higher levels, suggesting that they preexist in a Notch-high state (Figure 3C).

In considering whether Notch signaling is required for the slow-cycling state, we utilized a small molecule inhibitor of γ-secretase, a protease necessary for Notch activation (Kovall and Blacklow, 2010). Treatment with γ-secretase inhibitor (Compound E) reduced GSC8^Per growth, but had little effect on naïve GSC8 (Figure 3D). Notch signaling can also be suppressed genetically using a dominant negative form of mastermind-like (MAML) proteins (dnMAML), critical co-factors for Notch-mediated gene activation (Kovall and Blacklow, 2010; Maillard et al., 2004). Overexpression of dnMAML markedly reduced GSC8^Per proliferation, but had little effect on naïve GSC8 (Figure 3E). Conversely, overexpression of N1ICD was sufficient to induce a reversible growth reduction in GSC8 naïve (Figure 3F and S4D–S4F). These data suggest that Notch activation allows GSCs to transition to a slow-cycling, RTK-independent state.

Primary GBM Tumors Contain Notch-Positive Cells Depleted for Proliferation Markers

To address the in vivo relevance of our findings, we co-stained primary GBM specimens from PDGFRA- and EGFR-amplified tumors for the Notch marker N1ICD and the proliferation marker Ki67. A subset of N1ICD⁺ cells is evident in variable proportions (4.6 to 13%) in all specimens examined. Importantly, the majority of N1ICD⁺ cells (90–100%) were not Ki67⁺ (Figure 4A), suggesting that tumor cells with active Notch signaling are slow-cycling.

We also integrated expression signatures from our in vitro GSC study with single-cell RNA-seq data for primary GBM tumors. We acquired RNA-seq profiles for 300 individual cells from one EGFR-amplified tumor (MGH66), and integrated published single-cell data from two additional EGFR-amplified tumors (Patel et al., 2014). We considered the following 5 gene signatures: (i) a naïve GSC signature comprising genes in clusters 1–3 (Figure 2A); (ii) a persister GSC signature comprising genes in clusters 4–6 (Figure 2A); (iii) a cell cycle signature; (iv) a Notch-target gene signature derived by mapping N1ICD binding sites (see below); and (v) an in vivo qNSC signature (Codega et al., 2014).

We scored each individual tumor cell for each of these 5 signatures. Signatures for naïve GSCs and cell cycle were commonly expressed in a small subset of tumor cells (Figure 4B), consistent with relatively low Ki67 positivity. Signatures for N1ICD-target genes and in vivo qNSCs scored in a distinct subset of persister-like tumor cells (Figure 4B–4C). Altogether, these data indicate that primary GBMs contain a subset of cells that shares features with persister GSCs, including Notch-positivity, stem-like gene expression, and depletion of proliferation markers.

Notch Activity Targets Developmental Enhancers

To investigate how Notch signaling sustains persister GSCs, we mapped N1ICD and its co-factor RBPJ by ChIP-seq across the three GSC8 states. Consistent with increased levels of N1ICD and RBPJ, we detected a ~6-fold increase in the number of N1ICD binding sites in GSC8^12d and GSC8^Per, relative to GSC naïve (2,384 in GSC8^Per vs. 370 in GSC naïve), as well as increased signal intensities over these sites (Figure 5A-5B). As expected, N1ICD binding sites largely overlapped RBPJ binding sites (1750/2384 in GSC8^Per), and the underlying DNA sequences are enriched for RBPJ recognition motifs (Figure 5C). Consistent with prior studies in lymphoid cells (Wang et al., 2011), N1ICD binds primarily to distal sequences containing the active enhancer mark H3K27ac (Figure 5D–5E and S5A) (ENCODE Project Consortium et al., 2012).

To identify putative Notch target genes, we assigned N1ICD peaks detected in GSC8 naïve, GSC8^12d, or GSC8^Per to the nearest gene promoter. We focused on 1,696 target genes identified in GSC8^12d and/or GSC8^Per, but not in GSC8 naïve (Table S3). Consistent with a primary role in transcriptional activation, these candidate N1ICD targets were predominantly upregulated in GSC8^12d and GSC8^Per, relative to GSC8 naïve (Figure 5E–5F), and enriched for Notch pathway and neurodevelopmental genes as well as various canonical Notch targets (Figure 5D and S5B–S5C). Illustrative examples are the HEY1 and HES5 loci (Figure 5D): HEY1 is associated with poor survival in GBM patients (Hulleman et al., 2009) and HES5 is implicated in qNSCs (Imayoshi et al., 2010). In addition to canonical targets, N1ICD bound H3K27ac-marked distal elements in the vicinity of many neurodevelopmental master regulators (e.g., SOX TFs, OLIG1/2) and several pro-survival genes (Figure S5D). These N1ICD and RBPJ binding patterns support a direct mechanistic role for Notch signaling in the activation of regulatory elements and neurodevelopmental genes in slow-cycling GSC persisters.

Dynamic Histone Acetylation and Demethylation Accompany GSC^Per Switch

The dynamic acetylation of N1ICD-bound elements prompted us to investigate global chromatin alterations associated with the persister transition. We mapped H3K27ac genome-wide in naïve GSC8, GSC8^12d, and GSC8^Per, and called enriched intervals, identifying an average of 12,938 sites in each GSC state (naïve: 12,412; GSC8^12d: 12,940; GSC8^Per: 13,462). Approximately 63% of these sites are >5 kb away from any annotated transcriptional start site, thus representing candidate distal regulatory elements or enhancers. These enhancer-like elements are dynamic in activity (~25% common to all states, ~50% specific to one state).

We used K-means clustering to distinguish sets of promoter- and enhancer-like elements with related patterns of activity across the three GSC8 states (Figure 6A). We scanned the DNA sequences underlying the corresponding clusters for over-represented TF motifs (Figure 6A). The resulting predictions of differential TF activity are largely concordant with the TF expression patterns derived by RNA-seq (Figure 2A). GSC8 naïve-specific elements (Figure 6A, cluster II-III) are enriched for MEF2 motifs, consistent with high expression of MEF2C and its dimerization partner ASCL1 in the naïve state (Figure 2A, cluster 3) (Black et al., 1996). ASCL1 has been previously implicated in proliferative stem cell populations (Rheinbay et al., 2013). GSC8^12d-specific elements (Figure 6A, cluster IV-V) are enriched for Foxo motifs, in agreement with upregulation of FOXO3 (Figure 2A: cluster 5, and S3D), a TF with established roles in qNSCs (Paik et al., 2009; Renault et al., 2009). GSC8^Per-specific elements (Figure 6A, cluster VI-VII) are enriched for SOX and NFI motifs, consistent with increased expression of these TFs in GSC8^Per. Thus, concordant differences in regulatory element activity and TF expression distinguish the respective GSC states, and implicate neurodevelopmental regulatory programs in persister GSCs.

Many neurodevelopmental and Notch pathway genes are established targets of Polycomb repressors and the associated histone mark H3K27me3, which restrain their activity in non-expressing cell types (Simon and Kingston, 2013). We noted that the H3K27 methyltransferase EZH2 is downregulated in GSC persisters, while the H3K27 demethylases KDM6A and KDM6B are amongst a large set of KDMs that are upregulated (Figure 2A and Figure S3G–S3H). We therefore mapped and compared the distribution of H3K27me3 across naïve GSC8, GSC8^12d, and GSC8^Per. Many N1ICD-associated loci marked by H3K27me3 in naïve GSC8 undergo coincident H3K27 demethylation, upregulation, and H3K27 acetylation upon transition to the persister state (Figure 6D and Figure S5E). Examples include the HEY1 and HES5 loci, where increased binding of N1ICD and RBPJ is accompanied by gain of H3K27ac and loss of H3K27me3 (Figure 5D).

Broadly, we considered whether increased H3K27ac at enhancer-like elements in GSC8^Per are associated with a reduction in H3K27me3. Indeed, many of these elements are initially marked by H3K27me3 but are strongly depleted for this repressive mark upon activation in GSC8^12d and/or GSC8^Per cells (Figure 6A–6C cluster VII). Consistent with established roles for the respective modifications, genes near such regulated enhancer-like elements have higher expression in the persister cells (Wilcoxon test; P < 10⁻¹⁶).

Comparison of H3K27me3 profiles revealed an asymmetry between the GSC states. Whereas GSC^Per-specific H3K27ac-marked elements tend to have H3K27me3 in naïve GSCs (Figure 6A–B, cluster VII), naïve GSC-specific H3K27ac elements have diminished H3K27me3 in GSC^Per (Figure 6A–B, clusters II-III). This suggests that the transition from naïve to persister state is accompanied by a global loss of H3K27me3 from cis-regulatory elements. To gain more systematic insight, we called H3K27me3-enriched intervals genome-wide in each of the three GSC cell states. Remarkably, H3K27me3 levels across these intervals are markedly reduced in GSC8^12d and GSC8^Per, relative to naïve (Figure 6E). A reduction in H3K27me3 levels also occurred following treatment of GSC8 naïve cells with MEK inhibitor PD0325901 (Figure S6A) and, to a lesser extent, following treatment of GSC26 cells with an EGFR inhibitor (Figure S6B). Altogether, these data suggest that the activation of neurodevelopmental loci and Notch pathway genes is accompanied by pervasive acetylation of cis-regulatory elements, and may be facilitated by a global redistribution of H3K27me3 repressive chromatin.

H3K27me3 Demethylases KDM6A/B are Essential for GSC Persisters

Multiple KDMs are coordinately upregulated in the persisters (Figure 2A), while Notch-high tumors and Notch-high tumor compartments also exhibit significantly higher KDM expression (Brennan et al., 2013; Patel et al., 2014; Verhaak et al., 2010) (Figure 4C–D). Moreover, widespread redistribution of H3K27me3 in cells exiting cycle is suggestive of active histone demethylation. Taken together, these data implicate the H3K27 demethylases KDM6A/B in the transition or maintenance of persister GSCs. We therefore separately knocked out KDM6A and KDM6B in GSC8 naïve and GSC8^Per using CRISPR-Cas9 genome editing (Doench et al., 2014) and single-guide RNAs (sgRNAs) targeting their catalytic domains (Shi et al., 2015) (Figure S7A).

KDM6A and KDM6B knockout minimally affected GSC8 naïve cell growth (Figure 7A, left panel). By contrast, GSC8^Per cells were highly sensitive to knockout of either demethylase (Figure 7A, right panel); growth reduction afforded by KDM6B knockout was most pronounced, consistent with a reported pro-oncogenic role for KDM6B (Hashizume et al., 2014; Ntziachristos et al., 2014). Moreover, GSC8 naïve cells lacking KDM6A or KDM6B were less effective at forming persisters upon dasatinib treatment (Figure 7B).

To further explore this KDM6 dependency, we expanded our studies to other GSC models: GSC4 and GSC87. GSC4 is MYC-amplified, rapidly proliferating, and expresses a cell cycle gene signature akin to naïve GSC8 cells, while GSC87 is slow-cycling and expresses a persister-like gene signature at baseline (Figure S7B). GSC87 is highly CD133⁺, Notch-dependent, insensitive to dasatinib despite focal PDGFRA amplification, and displays similar H3K27ac profiles to GSC8^Per (Figure 1B, S4C, and S7B–S7D). Knockout of either KDM6A or KDM6B had negligible effects on GSC4, but KDM6B knockout significantly impaired persister-like GSC87 (Figure 7C). These data support general roles for KDM6 demethylases in the transition and maintenance of slow-cycling, persister-like GSCs.

To evaluate the contributions of KDM6 enzymatic activity, we treated GSCs with the small molecule inhibitor ‘GSKJ4’ (Kruidenier et al., 2013) (Figure 7D–7E). GSKJ4 treatment (1.5 µM) had negligible effects on the preexisting H3K27me3 landscape of naïve GSC8 but rescued the H3K27me3 loss seen upon dasatinib treatment (Figure 7E and S6C). GSKJ4 treatment did not significantly increase H3K4me3 in these models, supporting its selectivity for KDM6A/B at the doses used (Heinemann et al., 2014) (Figure S6D–S6E). These results support the notion that KDM6 contributes to the global H3K27me3 alterations in dasatinib-treated GSCs.

Consistent with our knockout data, GSKJ4 compromised proliferation of GSC8^12d and GSC8^Per at doses tolerated by naïve GSC8 (Figure 7D). Likewise, GSC8 persisters derived from sustained MEK inhibition are also preferentially sensitive to GSKJ4 (Figure S7E). By contrast, GSC persisters displayed no preferential sensitivity to GSKJ5, an inactive isomer of GSKJ4 (Kruidenier et al., 2013), nor to ‘KDM5-C70’, a recently described KDM5-selective inhibitor (Johansson et al., 2016) (Figure S7E–S7F). Lastly, we found that GSC87, which exhibits persister-like features at baseline, also showed increased sensitivity to GSKJ4, relative to the PDGFRA-dependent line CW1691 (Figure S7G). Thus, preferential KDM6 dependency may be a general characteristic of persister-like GSCs. Altogether, these data suggest that reconfiguration of H3K27me3 landscapes by histone demethylases facilitate the activation of regulatory elements and neurodevelopmental gene expression programs required for the drug tolerant persister GSC state (Figure 7F).

GSC Derivation and Cell Culture

Patient-derived GBM culture lines were obtained from Massachusetts General Hospital (GSC4, GSC8, GSC26, GSC87, GSC114, GSC125) (Wakimoto et. al 2009) and Case Western Reserve University (CW1691, CW2907, CW3246). For CW1691, CW2907 and CW3246, GBM tissues were obtained from excess surgical materials from patients after review from a neuropathologist and used in accordance with an approved protocol by the Institutional Review Board at Cleveland Clinic 2559; GBM cells were derived immediately after dissociation of primary patient tumor (CW1961) or after transient xenograft passage (CW2907). GBM lines were maintained in Neurobasal medium (Life Technologies) supplemented with N2/B27, penicillin/streptomycin (Life Technologies), GlutaMAX (Life Technologies), recombinant human EGF (20 ng/mL, R & D systems), and recombinant human FGF2 (20 ng/mL, R & D systems). GSC114, GSC115, CW1691, CW2907, and CW3246 were additionally supplemented with recombinant human PDGF-AA and PDGF-BB (20 ng/mL each, Shenandoah Biotechnology). CW1691 was grown as adherent monolayers using laminin (5 µg/mL Engelbreth-Holm-Swarm laminin, Sigma).

GSC Persisters

For drug persister studies, fresh inhibitor and media were replenished every 4–6 days to GSC cultures. For downstream studies (e.g., ChIP, immunoblot, gene expression), viable cells were enriched using Lympholyte (Cedarlane) to viability levels comparable to untreated cultures.

Primary Patient Specimens

Single cell RNA-seq data for MGH26, MGH30, and MGH31 was taken from prior published studies (Patel et. al. 2014). Slides for immunohistochemistry for these tumors were cut from archived blocks kept in the MGH pathology tumor bank under Partners Institutional Review Board Protocol 2011P002334. For MGH64 and MGH66, samples were acquired fresh at time of surgical resection at the Massachusetts General Hospital. Patients were consented preoperatively according to the Institutional Review Board Protocol 1999P008145.

Cell Growth Assays

Freshly dissociated single cell suspensions were plated (96-well) in triplicate or quadruplicate at 1,000 to 10,000 cells/well for testing. For 4 day growth assays, CellTiter-Glo (Promega) was added to wells and end point luminescence was measured (BioTek Synergy HTX Platereader). For 12–15 day assays, 1 × compound in media was added at day 4–5 and all wells were replated with fresh media and compound at day 8–10 before viability was measured at day 12–15. For each inhibitor, n=3–5 replicates were used for each concentration, and each inhibitor was repeated in two to four independent experiments.

CRISPR-Cas9 and shRNA Experiments

CRISPR sgRNA sequences were designed according to Doench et. al. (Doench et al., 2014), selected to target the demethylase catalytic domain (Shi et al., 2015), and subcloned into lentiCRISPR v1. sgRNA targeting sequences are included in Table S4. Lentiviruses were produced using standard protocols. Briefly, CRISPR plasmids were cotransfected with GAG/POL and VSVG plasmids into 293T packaging cells using FuGENE HD (Promega) to produce virus. Viral supernatant was collected 72 h after transfection and concentrated using Lenti-X Concentrator (Clontech) following the manufacturer’s instructions. Approximately 1.25×10⁶ adherently attached GSC cells (5 µg/mL Engelbreth-Holm-Swarm laminin, Sigma) were infected with concentrated virus for 24 h, and were selected 48 h after infection with puromycin (GSC4, GSC8: 1 µg/mL, GSC87: 2 µg/mL, Life Technologies) for 4 days. After 4–5 days of recovery, cells were dissociated and plated (96-well) in triplicate or quadruplicate at 2,500 to 5,000 cells/well. Growth was normalized to cells plated on day zero using an ATP standard curve measured at each time point. Each experiment was repeated in two to four independent experiments. Efficiency of genome editing was assessed by PCR amplification (PCR SuperMix, Life Technologies) of 50 ng of genomic DNA using primer sets surrounding the sgRNA-targeted region (~800 to 1,000 bp), followed by subjection of the resultant PCR products to SURVEYOR analysis according to the manufacturer’s protocol (Integrated DNA Technologies, 706020). PCR primers for SURVEYOR analysis are listed in Table S4. For shRNA knockdown experiments, the following lentiviral shRNAs were obtained from the RNAi Consortium in the pLKO.1 vector: PDGFRA (TRCN0000195132, TRCN0000196928, TRCN0000426196, TRCN0000419988) and GFP. Virus generation and growth testing for shRNA experiments followed essentially identical protocols outlined above except the growth assay was begun immediately after puromycin selection due to the immediacy of the phenotype. shRNA infection was performed in two biological replicates and the cell growth assay was performed with four replicates.

Overexpression Studies

The PDGFRA M260I ORF was PCR-amplified from p-DONR223-PDGFRA M260I using primers containing flanking Gateway attB sites and a stop codon, and then subcloned into p-DONR221. Site directed mutagenesis was performed using QuickChange Lightning Site-Directed Mutagenesis Kit according to the manufacturer’s instructions to generate wild-type PDGFRA (I260M) and PDGFRA D842V, which were subsequently subcloned into p-Inducer 20 (pIND20). Notch1 ICD ORF was PCR-amplified from Notch1 intracellular domain-pcw107 (Addgene, plasmid # 64621) using primers containing flanking Gateway attB sites and was subcloned into p-DONR221 and pIND20. dnMAML-GFP ORF was PCR-amplified from MigR1-dnMAML (Maillard et al., 2004) using primers containing flanking Gateway attB sites and was subcloned into p-DONR221 and pIND20. Virus production and infection of cells was carried out as described above, and GSCs were selected with G418 (500 µg/mL) for at least 1 week. pIND20 PDGFRA wt/D842V infected cells were induced with doxycycline (dox, 1 µg/mL) for at least 3 days before cell growth assays were performed. pIND20 dnMAML and pIND20 GFP infected GSC cells were induced with dox (1 µg/mL) for 4 days before mixed in equal proportion with uninfected GSC cells. Dox induction and mixing was performed independently in triplicate. The ratio of GFP⁺ to GFP⁻ cells was measured by flow cytometry at the indicated timepoints. Experiment was performed in two biological replicates from the same pIND20 dnMAML and pIND20 GFP lines.

Flow Cytometric Analysis

For cell cycle analysis, ~1×10⁶ viable dissociated cells were plated and treated with 1 µM EdU for 2 h. Cells were then washed with PBS, incubated with Zombie Aqua viability dye (Biolegend) or Live/Dead fixable far red dead cell stain kit (Life Technologies) for 20 min, washed with PBS, and then fixed with ice-cold 70% ethanol. Staining for EdU was performed using the Click-iT EdU Flow Cytometry Assay kit (Life Technologies) according to the manufacturer’s protocol. DNA content was visualized using FxCycle Violet (Life Technologies). Ki67 (clone B56, BD Biosciences) was stained for 30 min at ambient temperature. For surface marker analysis, cells were washed with PBS and stained with antibodies against CD15 (VIMC6, Miltenyi Biotec), CD133/1 (AC133, Miltenyi Biotec), or PDGFRα (16A1, BioLegend) for 30 min at 4 °C . Gates were determined using an unst ained control, where ~2% of cells were considered positive. All experiments were performed with at least two biological replicates.

Immunohistochemistry and Spectral Imaging of Primary GBM Specimens

Immunohistochemical staining of formalin-fixed, paraffin-embedded sections of four primary GBM specimens was done on a Leica Bond III automated workstation. Dual staining for activated NOTCH1 (N1ICD) and Ki67 was performed following antigen retrieval using Bond Epitope Retrieval 2 solution for 40 min. Slides were first incubated with antibody against N1ICD (D3B8, Cell Signaling Technology, 4147) at a 1:50 dilution for 60 min at ambient temperature, followed by incubation with a rabbit-specific secondary antibody linked to horseradish peroxidase (Bond Polymer Refine Detection kit). N1ICD staining was developed with diaminobenzidine (Leica Detection kit). Slides were then incubated with a mouse monoclonal antibody specific for Ki67 (8D5) at a 1:8,000 dilution for 30 min. Binding of anti-Ki67 was detected using the mouse-specific Bond Polymer Refine Red Detection kit, which detects staining using Fast Red. Slides were then counterstained with hematoxylin, dehydrated, and coverslipped.

To generate spectral libraries, single-stained tissue sections were imaged using the Mantra multispectral imaging platform (PerkinElmer). The spectrally resolved individual profiles between 420–720 nm of 3,3′-diaminobenzidine (DAB; N1ICD), Fast Red (Ki67) and hematoxylin counterstain were used to deconvolute staining patterns in double-stained tissue sections. Five representative areas of each stained tissue section (n = 2 per patient sample) were imaged at 20´ magnification and deconvoluted using the Inform 2.1 software package (PerkinElmer).

Immunoblotting

Cells were lysed on ice using RIPA buffer (50 mM Tris-HCl pH 7.5, 150 mM NaCl, 0.5% sodium deoxycholate, 1% NP-40, 0.1% SDS) supplemented with fresh protease (HALT, Thermofisher; PMSF, 1 mM) and phosphatase inhibitors (Thermofisher). After 30 min, benzonase was added with MgCl₂ (1 mM) to clarify the cell lysate. Immunobloting was performed according to standard procedures.

Tumor Dissociation

MGH66 patient sample was confirmed to be glioblastoma by frozen section diagnosis at the time of resection. Tissue was minced into <1mm pieces with a scalpel and enzymatic ally dissociated using the Brain Tumor Dissociation Kit (Miltenyi Biotec). Digested tissue was strained through a 100 micron cell strainer, washed in PBS and then layered carefully onto a 5mL density gradient (Lympholyte-H, Cedar Lane labs), which was centrifuged at 2,000 rpm for 10 min at room temperature. Dead cells, debris, and RBCs were pelleted and live cells collected from the interface and used for downstream applications including staining, flow cytometry, and single cell analysis. Viability was measured using trypan blue exclusion, which confirmed >90% cell viability.

Fluorescence-activated cell sorting of primary tumors

Tumor cells were suspended in 1% bovine serum albumin in Hanks buffered saline solution (BSA / HBSS) for blocking, and stained first with CD45-Vioblue direct antibody conjugate (Miltenyi Biotec) for 30 min at 4 °C . After washing in cold PBS, cells were resuspended in 1 mL of BSA/HBSS containing 1 µM calcein AM (Life Technologies) and 0.33 µM TO-PRO-3 iodide (Life Technologies) and stained for 30 minutes prior to cell sorting. Cells were run on a FACSAria Fusion Special Order System (Becton Dickinson) using 488nm (calcein AM, 530/30 filter), 640nm (TO-PRO-3, 670/14 filter), and 405nm (Vioblue, 450/50 filter) lasers. Singlet gating was performed using strict forward scatter height versus area criteria to eliminate doublets. Viable cells were identified by staining positive with calcein AM but negative for TO-PRO-3. Single cells were sorted into 96-well plates containing cold buffer TCL buffer (Qiagen) containing 1% β-mercaptoethanol, snap frozen on dry ice, and then stored at –80 °C prior to whole t ranscriptome amplification, library preparation, and sequencing.

Single-cell whole transcriptome amplification, library construction, and sequencing

Libraries from isolated single cells were generated based on the Smart-seq2 protocol (Picelli et al., 2014) with the following modifications. Agencourt RNAClean XP beads (Beckman Coulter) were used to purify single cell RNA prior to oligo-dT primed reverse transcription with Maxima reverse transcriptase and locked TSO oligonucleotide. Amplification of cDNA libraries was done with 20 cycle PCR amplification using KAPA HiFi HotStart ReadyMix (KAPA Biosystems) with subsequent Agencourt AMPure XP bead purification as described. Library preparation by tagmentation was done using the Nextera XT Library Prep kit (Illumina) with custom barcode adapters (sequences available upon request). Libraries from 384 cells from MGH66 with unique barcodes were combined and sequenced using a NextSeq 500 sequencer (Illumina).

Real-Time Quantitative RT-PCR and RNA-seq Library Preparation

Whole RNA was extracted from 1–3×10⁶ cells using the QIAGEN RNeasy kit according to the manufacturer’s protocol. For qRT-PCR, total RNA was reverse-transcribed into cDNA (High Capacitiy cDNA Reverse Transcriptase Kit, Applied Biosystems) and qRT-PCR amplification was performed in triplicate using fast SYBR Green Master Mix (Life Technologies) with specific PCR primers for genes of interest and 18S as an endogenous control. qRT-PCR experiments were performed with two biological replicates. Relative quantification for each target was performed using the comparative Ct method (Applied Biosystems). Sequences for qRT-PCR primers used are listed in Table S4. For RNA-seq library preparation, Poly(A)⁺ RNA was enriched using magnetic oligo(dT)-beads (Life Technologies) and then ligated to RNA adaptors for sequencing. RNA-seq was performed with three biological replicates per GSC line and/or condition.

ChIP-seq

Briefly, formaldehyde-fixed cells were lysed and sheared (Branson S220) on wet ice. The sheared chromatin was cleared and incubated overnight at 4 °C with the following antibodies: H3K4me3 (Millipore, 07-473, Lot 2207275), H3K27me3 (Millipore, 07-449, lot 2382150), H3K27ac (Active Motif, 39133 lot 25812006), RBPJ (Abcam, ab25949, lot GR169397-1), and cleaved Notch1 (CST, 4147, lot #4). Antibody-chromatin complexes were immunoprecipitated with protein G magnetic Dynal beads (Life Technologies), washed, eluted, reverse crosslinked, and treated with RNAse A followed by proteinase K. ChIP DNA was purified using Ampure XP beads (Beckmann Coulter) and then used to prepare sequencing libraries for sequencing with the Next-Seq Illumina genome analyzer. Experiments were performed once.

Levels phosphorylated PDGFRα, Akt, and Erk1/2

For each assay n = 2 biological replicates involving separate drug treatments.

Cell growth assays (Cell-titer Glo)

N=3–4 replicates per cell line. N=2–4 independent experiments per compound, and 3–5 replicates per concentration (see figure legends for more detail). For overexpression studies n= 2. Relative CTG (Cell-tier Glo) values were normalized to T₀ (y-axis) over a time course (x-axis) following treatment. Error was calculated as either s.d or s.e.m across all biological replicates. IC₅₀ values: Calculated as standard, from at least three biological replicates.

CRISPR-Cas9 and shRNA experiments

Each CRISPR-Cas9 sgRNA experiment was repeated in two to four independent experiments involving separate infections; in each experiment, cells were plated in n = 3–5 replicates for the cell growth assay. shRNA infection was performed in two biological replicates and the cell growth assay was performed with n = 4 replicates. Growth was normalized to cells plated on day zero using an ATP standard curve measured at each time point.

Overexpression studies

Dox induction and mixing was performed independently in triplicate. The ratio of GFP⁺ to GFP⁻ cells was measured by flow cytometry at the indicated timepoints, and error bars were calculated as s.e.m. across all three biological replicates (indicated in the figure legends). Experiment was performed in two biological replicates from the same pIND20 dnMAML and pIND20 GFP lines.

Flow Cytometric Analysis

PDGFRα staining was performed with two biological replicates, while staining for CD133 and CD15 were performed with three biological replicates. Gates were determined using an unstained control, where ~2% of cells were considered positive. Error bars were calculated as s.e.m. across three replicates and significance was determined using a two-tailed Student’s T-Test. Cell cycle analysis was performed with at least three biological replicates per time point, and gates were determined for each timepoint by the relative proportions of DAPI and EdU labelling.

Immunohistochemistry and Spectral Imaging of Primary GBM Specimens

Single slides from four tumors were stained with antibodies against Notch1 and Ki67 (see above methods) in two separate replicates. From each slide, five representative regions were quantified using spectral imaging, yielding >10,000 cells quantified per sample over two replicates.

Single-cell RNA-seq data analysis

Libraries from 384 cells from a single tumor (MGH66) were combined and sequenced. Paired-end, 38-base reads were mapped using Bowtie (Langmead et al., 2009) with parameters “-q --phred33-quals -n 1 -e 99999999 -l 25 -I 1 -X 2000 -a -m 15 -S -p 6” to the UCSC hg19 human transcriptome, which permits alignment of sequences with single base changes to account for point mutations. RSEM v1.2.3 (Li and Dewey, 2011) in paired end mode using parameters “--estimate-rspd --paired end -sam -p 6” was used to quantify expression in the form of normalized transcripts per million (TPM) for each gene in each single cell. Data for tumors MGH26, MGH28, and MGH30 were prepared as previously published in (Patel et al., 2014).

To remove low quality data for each tumor, we excluded all cells with less than 3,000 detectable genes and all genes with a mean TPM less than 10. The TPM values were subjected to base 2 logarithmic transformation after adding 1 to each value to avoid singularities. To account for global mean expression differences across cells, for each cell a size factor was calculated as described in (Anders and Huber, 2010) with the following modifications: i) computations were performed in log-space and adjusted accordingly; ii) the mean over all genes as opposed to the median was used ultimately. The cell-specific size factors were subtracted from the transformed TPM values. To assess the cell cycle state we used the cell cycle signature derived in our previous study (see (Patel et al., 2014) for details). To score single-cells for individual signatures (e.g. cell cycle or persister signatures) we computed mean TPM values for the corresponding signature gene sets and transformed resulting values into standard scores (z-scores). Gene signatures used are listed in Table S3.

Population RNA-seq analysis

Paired-end reads were aligned to UCSC transcriptome (hg19) using Bowtie (Langmead et al., 2009) (version 0.12.7) with the following parameters: --chunkmbs 512 -q --phred33-quals -n 0 -l 25 -I 1 -X 2000 -p 6 -a -m 15. Gene expression, quantified as transcript per million (TPM), was estimated using RSEM (Li and Dewey, 2011) with the following parameters: --fragment-length-max 1000 --estimate-rspd --paired-end.

Data processing was performed using R version ≥3.1.0 (R Core Team, 2014; http://www.Rproject.org/) and BioConductor (http://www.bioconductor.org). All heatmaps were generated using the R package ggplot2 (Wickham 2009). To identify variable genes across the time course of drug treatment, triplicate RNA-seq data sets were generated for each time point and processed as outlined above. The resulting TPM values were employed for differential gene expression analysis using the R package EBSeq (Leng et al., 2015). To identify variable genes the following pairwise comparisons were performed: GSC8 naïve versus GSC8^12d, GSC8 naïve versus GSC8^Per, GSC8^12d versus GSC8^Per. We considered all genes with a posterior probability of differential expression (PPDE) greater than 0.5 in at least one of the comparisons as variable genes. A table of all genes and significantly differentially expressed genes is provided in Table S2. GSEA version 2.1 (Subramanian et al., 2005) was carried out using signal-to-noise on log2 + 1 transformed TPM values as the metric; genes with mean TPM less than 10 under all given conditions were excluded from analysis. Gene signatures used for GSEA are listed in Table S3. For the comparison with GTEx expression data (GTEx Consortium, 2013), RPKM values were computed.

ChIP-seq data analysis

All experiments were performed once. Reads were aligned to hg19 using BWA (Li and Durbin, 2009) and identical ChIP-seq sequence reads were collapsed to avoid PCR duplicates. In order to avoid possible saturation biases, reads were downsampled to approximately similar numbers (~20 million reads). Peaks were called using HOMER v4.6 (Heinz et al., 2010) using matched inputs with the following parameters: H3K4me3, –histone –tagThreshold 30; H3K27ac, –histone –tagThreshold 50; H3K27me3, –histone –size 3000 –minDist 2500 –F 1.5 –L 0 –FDR 0.1; Notch1 ICD, –factor; and RBPJ, –factor. TF motif enrichment analysis was performed using HOMER v4.6 on 1 kb windows centered on previously called H3K27ac, Notch1 ICD, and RBPJ peaks (parameters: –size given –mask) and p-values presented were provided by the software. Quantification of the H3K27me3 ChIP-seq signals, using fore- and background normalization, was essentially performed as outlined in (Zhu et al., 2013) with few modifications. To exclude genomic regions with potential copy number aberrations, we quantified read counts derived from whole cell extracts ChIP-seq within 5 kb genomic windows. Read counts were performed using the R package GenomicRanges (Lawrence et al., 2013). To compute normalization constants, all genomic windows with read counts equal to 0 or greater than 3 standard deviations from the mean were excluded from the analysis. For a given H3K27me3 ChIP-seq data set, the background constant was estimated as the median of the read count distribution within 5 kb windows overlapping the 5% most highly expressed genes, as active genes are considered to be devoid of repressive chromatin modifications. The foreground signal was estimated as the median of the read count distribution within 5 kb windows centered around each called peak and subtracting the background signal. To quantify H3K27me3 signal within 5 kb regions of interest, the background signal was subtracted from the corresponding read count and the resulting value divided by the foreground signal. Values smaller or larger than 0 or 1 were mapped to 0 or 1, respectively. The fore- and background normalization of H3K27ac was performed analogously but with the following modifications. To account for the narrower signal, all operations were performed using 1 kb windows. For a given H3K27ac ChIP-seq data set, the background signal was estimated as the median of the read count distribution within all 1 kb windows across the genome but not overlapping with a called peak. To account for greater dynamic range, the foreground signal was estimated as the 0.95-quantile of the read count distribution within 1 kb windows centered around each called peak and subtracting the background constant. To estimate the H3K27ac signal within a 1 kb region of interest, the background constant was subtracted from the corresponding read count and the resulting value was divided by the foreground constant. Values smaller or larger than 0 or 1 were mapped to 0 or 1, respectively. In order to quantify ChIP-seq intensities within genomic regions of different sizes (e.g. genes), the obtained normalization constants were scaled accordingly. Association of peaks with neighboring genes was performed using the R package ChIPpeakAnno (Zhu et al., 2010).

ChIP-seq Binding Profiles

ChIP-seq profile plots were generated using ngs.plot (Shen et al., 2014). To assess H3K27me3 levels across broad domains, H3K27me3 metaprofiles were generated over the union of peaks, called in each of the corresponding GSC8 persister conditions, and having a minimum size of 10 kb.