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 nongenetic cell plasticity. Our findings demonstrate an approach to uncover TRN topology and use it to rationally predict combinatorial treatments that disrupt 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 responders (blue) and nonresponders (red) to pitavastatin. Below are drug-dose response and time-course response curves for SN520 (pitavastatin-responsive) and SN503 (pitavastatin-nonresponsive) 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 gray). (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 datasets 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. Scatterplots show proportions of each subtype in each PD-GSC population across treatment for (C) SN520 and (D) SN503. (E and F) Flow cytometry analysis of PN and MES markers CD133 (PN) and CD44 (MES) across pitavastatin-treated cells for SN520 and SN503, respectively. Values (gray) indicate percentages of cell populations in each quadrant. Proportions of cells positive for each subtype marker are quantified in the adjacent bar plots underneath. (G and 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 area under the curve (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 up-regulated 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 ≤ 1 × 10⁻⁵) (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) Dot plot 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 up-regulated 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 TRN 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 (Tr. Pr.), while single cells (columns) are organized into transcriptional states (Tr. 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 bar plot 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 interquantile range (IQR; 25th and 75th percentile) and median activity value, while the whiskers represent 1.5× IQR. Asterisks indicate statistically significant differences between regulon activity distributions. Single asterisks (*) denote that activity distribution of untreated controls (CTRL) is significantly lower than distribution being compared (FDR << 1 × 10⁻³). Double asterisks (**) denote that distribution of untreated controls is significantly higher than either vehicle-treated (DMSO) or pitavastatin-treated (PSTAT) distributions being compared (FDR << 1 × 10⁻³). (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 (second row), pseudotime (third row), and RNA velocity (fourth 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 gray 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 pitavastatin-treated cells (PSTAT) relative to untreated control (CTRL) counterparts (Wilcoxon rank test, FDR << 1 × 10⁻⁵). The box represents the IQR (25th and 75th percentile) and median activity value, while the whiskers highlight 1.5× 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. Scatterplot of simulated states projected along first two PCs. Horizontal bar plot 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 and F) Corresponding simulation results for SN503. (G) SN520 cell viability following 4-day treatment with either simultaneous treatment with pitavastatin and siRNA (light gray bars) or sequential pitavastatin and then siRNA-mediated knockdown of TFs (dark gray bars). Viabilities are relative to nontemplate control (NTC)–treated cells. (H) Corresponding bar plots of relative viability for SN503. Asterisks [(G) and (H)] indicate significant decrease relative to corresponding NTC treatment (FDR P ≤ 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 ≤ 1 × 10⁻²⁰, **FDR ≤ 1 × 10⁻¹⁵⁰). (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 pretreated with vehicle (DMSO, light blue)/pitavastatin (2 μM, pink), followed by 24-hour vinflunine treatment (1.5 × 10⁻⁹, 4.6 × 10⁻⁹, 13.7 × 10⁻⁹, 41.2 × 10⁻⁹, 123.5 × 10⁻⁹, 370.4 × 10⁻⁹, 1.10 × 10⁻⁶, 3.30 × 10⁻⁶, 10.0 × 10⁻⁶, 30.0 × 10⁻⁶ M). Results from 48-hour vinflunine treatment are included in fig. S16. Adjacent bar plots show relative viabilities following various treatments (black dots underneath bar plots) including monotherapy with pitavastatin (PIT) or pretreatment 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 ±2× SD (N = 3). (E) Depiction of how core TF-TF networks underlying drug response drive cell-state transitions in responder and nonresponder 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 nonresponder PD-GSC consists of cell subpopulations (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.