Gene regulatory network topology governs resistance and treatment escape in glioma stem-like cells

Abstract

Poor prognosis and drug resistance in glioblastoma (GBM) can result from cellular heterogeneity and treatment-induced shifts in phenotypic states of tumor cells, including dedifferentiation into glioma stem-like cells (GSCs). This rare tumorigenic cell subpopulation resists temozolomide, undergoes proneural-to-mesenchymal transition (PMT) to evade therapy, and drives recurrence. Through inference of transcriptional regulatory networks (TRNs) of patient-derived GSCs (PD-GSCs) at single-cell resolution, we demonstrate how the topology of transcription factor interaction networks drives distinct trajectories of cell state transitions in PD-GSCs resistant or susceptible to cytotoxic drug treatment. By experimentally testing predictions based on TRN simulations, we show that drug treatment drives surviving PD-GSCs along a trajectory of intermediate states, exposing vulnerability to potentiated killing by siRNA or a second drug targeting treatment-induced transcriptional programs governing non-genetic cell plasticity. Our findings demonstrate an approach to uncover TRN topology and use it to rationally predict combinatorial treatments that disrupts acquired resistance in GBM.

Fig. 1.. Pitavastatin causes shift in molecular subtype expressed by PD-GSCs.

(A) Pitavastatin IC₅₀ values for each of 45 PD-GSCs as determined using dose titration assays (below). Labeled PD-GSCs represent a subset deemed as a responders (blue) and non-responders (red) to pitavastatin. Below are drug-dose response and time-course response curves for SN520 (pitavastatin-responsive) and SN503 (pitavastatin-non-responsive) PD-GSC populations. (B) Experimental workflow for longitudinal monitoring of PD-GSC response to pitavastatin treatment. Colored horizontal arrows indicate duration of pitavastatin (magenta), vehicle-control (DMSO, light blue), or untreated control (dark grey). (C) GSVA enrichment scores for each molecular subtype (CL – classical, PN – proneural, MES – mesenchymal) analyzed for all bulk samples collected. (D) UMAP plots of Harmony-integrated scRNA-seq data sets and corresponding individual plots for each PD-GSC phenotype treated with DMSO or pitavastatin (PSTAT) and untreated controls (CTRL) representing D0 time point. (E) Wasserstein distance of transport distances between each consecutive time point for each PD-GSC under each treatment condition (vehicle- or pitavastatin-treatment).

Fig. 2.. Single-cell characterization of PD-GSC response to pitavastatin.

UMAP plots of scRNA-seq profiles, annotated according to treatment conditions (untreated control, vehicle – DMSO, and pitavastatin – PSTAT), for (A) SN520 and (B) SN503. Scatter plots show proportions of each subtype in each PD-GSC population across treatment for (C) SN520 and (D) SN503. (E – F) Flow cytometry analysis of PN and MES markers CD133 (PN) and CD44 (MES) across pitavastatin-treated cells for SN520 and SN503, respectively. Values (grey) indicate percentages of cell populations in each quadrant. Proportions of cells positive for each subtype marker are quantified in the adjacent barplots underneath. (G – H) Heatmap of inferCNV scores for SN520 and SN503, respectively. Cells (rows) are grouped based on treatment conditions (same color annotation as in (A) and (B)). Genes (columns) are arranged according to their chromosomal positions. (I) Dose-response curves of naïve SN520 PD-GSCs (light blue) and SN520 PD-GSCs that survived an initial pitavastatin-treatment (treated – dark blue). Adjacent plot shows corresponding AUC values from dose-response curves generated from subsequent PD-GSC cultures derived from original pitavastatin- or vehicle-control-treatment for SN520 (left) and SN503 (right). Paired t-test results showed a sustained (significant) increase in AUC values of the PSTAT-treated SN520 PD-GSCs relative to their vehicle-control counterparts but not for SN503.

Fig. 3.. Differential expression and pathway enrichment analysis reveals underlying processes driving pitavastatin responses.

(A) Heatmap of the top upregulated DEGs, based on FDR p-values, across the 14 Louvain cell clusters (cl) identified in vehicle-control- and pitavastatin-treated SN520 PD-GSCs. Adjacent UMAP plot with treatment annotation (same as Fig 2A) included for reference. (B) Corresponding UMAP plots of scRNA-seq profiles annotated according to Louvain cell cluster (left) and treatment condition (right) as reference. (C) Cell proportions for each Louvain cluster that belong to each treatment condition for SN520. Significant enrichment of treatment condition within Louvain cluster indicated by asterisk (FDR ≤ 0.05) or double dagger (FDR ≤ 1e-05) (D) Cell proportions for each Louvain cluster that belong to each treatment condition for SN503. Significant enrichment notation identical to that used in (D). (E) Dotplot of hallmark gene sets enriched across SN503 and SN520 PD-GSCs, grouped with respect to either drug-treatment duration or Louvain clustering. Dot size represents the ratio of number of upregulated genes associated with a PD-GSC grouping to the number of genes associated with a specific hallmark gene set. Dot colors indicate significance of enrichment (FDR value). (F) Total number of up- and down-regulated DEGs, relative to untreated control (D0) cells, at each treatment time point for SN503 (red) and SN520 (blue).

Fig. 4.. MINER3 transcriptional regulatory network inference reveals mechanisms of cell-state changes.

(A) Heatmaps of normalized regulon activities across SN520 (top) and SN503 (bottom) PD-GSCs. Regulons (rows) are organized into transcriptional programs (Pr) while single cells (columns) are organized into transcriptional states (St). Left-adjacent color bars indicate what regulons belong to a particular transcriptional program. Left-adjacent color bar indicates transcriptional programs. Top color bars indicate treatment condition (color annotation identical to Fig. 1E) and corresponding transcriptional state for a single cell. (B) Stacked barplot show proportion of cells within each transcriptional state from each treatment condition for SN520 (top) and SN503 (bottom). (C) Boxplot/violin plots of distributions of regulon activity for select programs across treatment conditions for SN520 and SN503. Regulon activity values were capped between the lower 2.5% and 97.5% range of values. Labels indicate program IDs and select hallmark gene sets (95) enriched within each program. The box represents the inter-quantile range (IQR – 25^th and 75^th percentile) and median activity value while the whiskers represent 1.5x IQR. Asterisks indicate statistically significant differences between regulon activity distributions. Single asterisks (*) denote activity distribution of untreated controls (CTRL) is significantly lower than distribution being compared (FDR << 1e-3). Double asterisks (**) denote distribution of untreated controls is significantly higher than either vehicle-treated (DMSO) or pitavstatin-treated (PSTAT) distributions being compared (FDR << 1e-3). (D) Flow diagram outlining approach to derive core TF-TF network from MINER3 results. Final core TF-TF networks derived for (E) SN520 and (F) SN503.

Fig. 5.. Distinct trajectories define SN520 and SN503 pitavastatin response.

(A) UMAP plots of vehicle- and pitavastatin-treated cells for SN520 (left column) and SN503 (right column). Annotation highlights treatment conditions (top row), molecular subtype (2^nd row), pseudotime (3^rd row) and RNA velocity (4^th row). (B) Critical transition index (I_c) of SN520 (blue) and SN503 (red) cells treated with vehicle (DMSO - light) or pitavastatin (PSTAT – dark). (C) LOESS regression of TF expression behavior sorted according to peak expression along pseudotime (Monocle3). Density plots depict distribution of sample time points along pseudotime trajectory. Heatmap shows expression of TFs rank sorted by time of peak expression along pseudotime (color bar beneath heatmap). (D) Select set of LOESS regression of mean program activities with respect to pseudotime. Regulons are clustered based on their dynamic activity profiles with respect to pseudotime. Dashed grey line represents the average shape of the curves for each cluster. Labels indicate which transcriptional programs were grouped into each cluster. Select hallmark gene sets (95) enriched within programs are labeled as well. (E) Boxplots/violin plots of expression of genes associated with indicated pathways/processes (95) on respective treatment days. Relative gene expression values were capped at the lower 2.5% and 97.5% range of values. Labels indicate select hallmark gene sets enriched within subpopulation of cells (treatment time point). Asterisks indicate statistically greater expression in pitavstatin-treated cells (PSTAT) relative to untreated control (CTRL) counterparts (Wilcoxon rank test, FDR << 1e-5). The box represents the inter-quantile range (IQR – 25^th and 75^th percentile), median activity value while the whiskers highlight 1.5x IQR.

Fig. 6.. Dynamic simulations of core TF regulatory network supports phenotypic plasticity of GSCs.

Simulated transcriptional states (black circles) projected along first two principal components. Contour lines represent distribution of PCA scores of TF expression states (core TFs only) for (A) SN520 and (B) SN503 cells. One thousand simulated states were generated using core TF network topologies and corresponding D0 scRNA-seq data for initial conditions (i.c.) as RACIPE inputs. (C) Three plots summarizing results from 1 million RACIPE simulations (independent of (A)) using the core TF-TF network derived from scSYGNAL-520 and randomized initial conditions to explore plausible steady states supported by the network topology. Dendrogram of four distinct simulated steady states. Scatter plot of simulated states projected along first two PCs. Horizontal barplot of rank-ordered TFs based on their importance in distinguishing the four simulated states. Here, importance is defined by the mean decrease in classification accuracy following TF removal from the model, per random forest analysis. (D) Heatmap of expression for SN520 core TFs. Cells (columns) were hierarchically clustered to define experimental states (ES₅₂₀-i), providing a basis of comparison for simulated states (SS₅₂₀-i). Adjacent boxplots of three TFs having high importance in random forest classification. Boxplots (top row) of TF expression distributions for experimental states. Boxplots (bottom row) of simulated TF expression distributions (normalized). (E – F) Corresponding simulation results for SN503. (G) SN520 cell viability following 4-day treatment with either simultaneous treatment with pitavastatin and siRNA (light grey bars) or sequential pitavastatin then siRNA-mediated KD of TFs (dark gray bars). Viabilities are relative to non-template control (NTC)-treated cells. (H) Corresponding bar plots of relative viability for SN503. Asterisks (G – H) indicate significant decrease relative to corresponding NTC treatment (FDR p-values ≤ 0.1).

Fig. 7.. Dynamics of regulon behavior reveal additional targets that guide rational secondary drug selection.

Distribution of activity of select tubulin-associated regulons in single cells across treatments for (A) SN520 and (B) SN503. Asterisks indicate treatments having significantly higher activities relative to the untreated control (CTRL (D0)) (Wilcoxon rank test, * FDR ≤ 1e-20, ** FDR ≤ 1e-150). (C) Experimental design for sequential pitavastatin/vinflunine treatment on multiple PD-GSCs. (D) Dose-response curves for SN520 and SN503 cells treated with pitavastatin alone (PIT, dark gray), or pre-treated with vehicle (DMSO, light blue)/pitavastatin (2μM, pink), followed by 24hr vinflunine treatment (1.5e-9, 4.6e-9, 13.7e-9, 41.2e-9, 123.5e-9, 370.4e-9, 1.10e-6, 3.30e-6, 10.0e-6 30.0e-6 M). Results from 48hr vinflunine treatment included in fig. S16. Adjacent barplots show relative viabilities following various treatments (black dots underneath barplots) including monotherapy with pitavastatin (PIT), or pre-treatment with DMSO (pre-DMSO)/pitavastatin (pre-PIT) followed by vinflunine (VIN). Asterisks/double crosses indicate treatments resulting in significantly lower relative viability than pitavastatin monotherapy (* 1.1μM, FDR ≤ 0.1; ‡ 3.3 μM FDR ≤ 0.1). Color annotation identical to dose-response curves. Error bars represent ±2x standard deviation (N = 3). (E) Depiction of how core TF-TF networks underlying drug-response drive cell state transitions in responder and non-responder PD-GSCs along a Waddington-like phenotypic landscape. Treatment with a primary drug to which cells are sensitive (1° drug^S) activates a highly interconnected network in a responder PD-GSC, driving PMT across surviving cells resulting in acquired resistance to “multiple drugs^R”. Intervention with a second drug (2° drug^S) that targets vulnerabilities in transient states potentiates killing and disrupts PMT. By contrast, the non-responder PD-GSC consists of cell sub-populations (center well) resistant to the primary drug (1° drug^R). Here, treatment with 1° drug^R activates a sparse network that drives surviving cells into multiple distinct drug-resistant states potentially sensitive to secondary interventions.