NF1 regulates mesenchymal glioblastoma plasticity and aggressiveness through the AP-1 transcription factor FOSL1

The molecular basis underlying glioblastoma (GBM) heterogeneity and plasticity is not fully understood. Using transcriptomic data of human patient-derived brain tumor stem cell lines (BTSCs), classified based on GBM-intrinsic signatures, we identify the AP-1 transcription factor FOSL1 as a key regulator of the mesenchymal (MES) subtype. We provide a mechanistic basis to the role of the neurofibromatosis type 1 gene (NF1), a negative regulator of the RAS/MAPK pathway, in GBM mesenchymal transformation through the modulation of FOSL1 expression. Depletion of FOSL1 in NF1-mutant human BTSCs and Kras-mutant mouse neural stem cells results in loss of the mesenchymal gene signature and reduction in stem cell properties and in vivo tumorigenic potential. Our data demonstrate that FOSL1 controls GBM plasticity and aggressiveness in response to NF1 alterations.


Introduction
Gliomas are the most common primary brain tumor in adults. Given the strong association of the isocitrate dehydrogenase 1 and 2 (IDH1/2) genes mutations with glioma patients survival, the 2016 WHO classification, which integrates both histological and molecular features, has introduced the distinction of IDH-wildtype (IDH-wt) or IDH-mutant (IDH-mut) in diffuse gliomas (Louis et al., 2016). IDH-wt glioblastoma (GBM) represents the most frequent and aggressive form of gliomas, characterized by high molecular and cellular inter-and intra-tumoral heterogeneity.
Large-scale sequencing approaches have evidenced how concurrent perturbations of cell cycle regulators, growth and survival pathways, mediated by RAS/MAPK and PI3K/AKT signaling, play a significant role in driving adult GBMs (Brennan et al., 2013;Cancer Genome Atlas Research Network, 2008;Verhaak et al., 2010). Moreover, various studies have classified GBM in different subtypes, using transcriptional profiling, being now the proneural (PN), classical (CL), and mesenchymal (MES) the most widely accepted (Phillips et al., 2006;Verhaak et al., 2010;Wang et al., 2017).
Patients with the MES subtype tend to have worse survival rates compared to other subtypes, both in the primary and recurrent tumor settings . The most frequent genetic alterations -neurofibromatosis type 1 gene (NF1) copy number loss or mutation -and important regulators of the MES subtype, such as STAT3, CEBPB, and TAZ, have been identified (Bhat et al., 2011;Carro et al., 2010;Verhaak et al., 2010). Nevertheless, the mechanisms of regulation of MES GBMs are still not fully understood. For example, whether the MES transcriptional signature is controlled through tumor cell-intrinsic mechanisms or influenced by the tumor microenvironment (TME) is still an unsolved question. In fact, the critical contribution of the TME adds another layer of complexity to MES GBMs. Tumors from this subtype are highly infiltrated by non-neoplastic cells as compared to PN and CL subtypes . Additionally, MES tumors express high levels of angiogenic markers and exhibit high levels of necrosis (Cooper et al., 2012).
Even though each subtype is associated with specific genetic alterations, there is a considerable plasticity among them: different subtypes coexist in the same tumors and shifts in subtypes can occur over time Sottoriva et al., 2013). This plasticity may be explained by acquisition of new genetic and epigenetic abnormalities, stem-like reprogramming, or clonal variation (Fedele et al., 2019). It is also not fully understood whether the distinct subtypes evolve from a common glioma precursor . For instance, PN and CL tumors often switch phenotype to MES upon recurrence, and treatment also increases the mesenchymal gene signature, suggesting that MES transition, or epithelial to mesenchymal (EMT)-like, in GBM is associated with tumor progression and therapy resistance (Bhat et al., 2013;Halliday et al., 2014;Phillips et al., 2006). Yet, the frequency and relevance of this EMT-like phenomenon in glioma progression remains unclear. EMT has also been associated with stemness in other cancers (Mani et al., 2008;Ye et al., 2015). Glioma stem cells (GSCs) share features with normal neural stem cells (NSCs) such as self-renewal and ability to differentiate into distinct cellular lineages (astrocytes, oligodendrocytes, and neurons) but are thought to be responsible for tumor relapse, given their ability to repopulate tumors and their resistance to treatment (Bao et al., 2006;Chen et al., 2012). GSCs heterogeneity is also being increasingly observed (Bhat et al., 2013;Mack et al., 2019;Richards et al., 2021), but whether genotype-to-phenotype connections exist remain to be clarified.
Here we report that NF1 loss, by increasing RAS/MAPK activity, modulates FOSL1 expression, which in turn plays a central function in the regulation of MES GBM. Using a surrogate mouse model of MES GBM and patient-derived MES brain tumor stem cells (BTSCs), we show that FOSL1 is responsible for sustaining cell growth in vitro and in vivo, and for the maintenance of stem-like properties. We propose that FOSL1 is an important regulator of GBM stemness, MES features and plasticity, controlling an EMT-like process with therapeutically relevant implications.

FOSL1 is a key regulator of the MES subtype
To study the tumor cell-intrinsic signaling pathways that modulate the GBM expression subtypes, we assembled a collection of transcriptomic data (both expression arrays and RNA-sequencing) of 144 samples derived from 116 independent BTSC lines (see Materials and methods for details). Samples were then classified according to the previously reported 50-gene glioma-intrinsic transcriptional subtype signatures and the single-sample gene set enrichment analysis (ssGSEA)-based equivalent distribution resampling classification strategy . Principal component analysis (PCA) showed a large overlap of the transcription profile among BTSCs classified either as CL/PN while most of the MES appeared as separate groups ( Figure 1A and Supplementary file 1). This separation is consistent with early evidence in GSCs (Bhat et al., 2013) and holds 92% of concordance in the identification of a recent two transcriptional subgroups classification of single-GSCs defined as developmental (DEV) and injury response (INJ) (Richards et al., 2021). Differential gene expression analysis comparing mesenchymal versus non-mesenchymal BTSCs confirmed the clear separation among the two groups, with only a minor fraction of cell lines showing a mixed Source data 1. Source data of Figure 1A, B, D-F.    expression profile ( Figure 1B and Supplementary file 2), further supporting that GSCs exist along a major transcriptional gradient between two cellular states (Bhat et al., 2013;Richards et al., 2021).
To reveal the signaling pathways underlying the differences between MES and non-MES BTSCs, we then applied a network-based approach based on the Algorithm for the Reconstruction of Accurate Cellular Networks (ARACNe) (Basso et al., 2005;Carro et al., 2010), which identifies a list of TFs with their predicted targets, defined as regulons. The regulon for each TF is constituted by all the genes whose expression data exhibit significant mutual information with that of a given TF and are thus expected to be regulated by that TF (Castro et al., 2016;Fletcher et al., 2013). Enrichment of a relevant gene signature in each of the regulons can point to the TFs acting as master regulators (MRs) of the response or phenotype (Carro et al., 2010;Fletcher et al., 2013). Master regulator analysis (MRA) identified a series of TFs, among which FOSL1, VDR, OLIG2, SP100, ELF4, SOX11, BNC2, ASCL1, SALL2, and POU3F3 were the top 10 most statistically significant (Benjamini-Hochberg p<0.0001) ( Figure 1C and Supplementary file 3). FOSL1, VDR, SP100, ELF4, and BNC2 were significantly upregulated in the MES BTSCs, while OLIG2, SOX11, ASCL1, SALL2, and POU3F3 were upregulated in the non-MES BTSCs ( Figure 1D and We then analyzed the CGGA and TCGA pan-glioma datasets (Ceccarelli et al., 2016;Zhao et al., 2017) and observed that FOSL1 expression is elevated in the IDH-wt glioma molecular subgroup ( Figure 1E and Supplementary file 4) with a significant upregulation in the MES subtype in bulk tumors, and it is also enriched in MES-like cells (Neftel et al., 2019) at the single-cell level (Figure 1-figure supplement 2A-C). Importantly, high expression levels were associated with worse prognosis in IDH-wt tumors ( Figure 1F), thus suggesting that FOSL1 could represent not only a key regulator of the glioma-intrinsic MES signature, but also a putative key player in MES glioma pathogenesis.
NF1 modulates the MES signature and FOSL1 expression NF1 alterations and activation of the RAS/MAPK signaling have been previously associated with the MES GBM subtype (Brennan et al., 2013;Verhaak et al., 2010;Wang et al., 2016;Wang et al., 2017). However, whether NF1 plays a broader functional role in the regulation of the MES gene signature (MGS) in IDH-wt gliomas still remains to be established.
We initially grouped, according to the previously described GBM subtype-specific gene signatures, a subset of IDH-wt glioma samples of the TCGA dataset for which RNA-seq data were available (n = 229) (see Materials and methods for details). By analyzing the frequency of NF1 alterations (either point mutations or biallelic gene loss), we confirmed a significant enrichment of NF1 alterations in MES versus non-MES tumors (Fisher's exact test p=0.0106) (Figure 2A, B). Importantly, we detected higher level of FOSL1 mRNA in the cohort of IDH-wt gliomas with NF1 alterations (Student's t test p=0.018) ( Figure 2C), as well as a significant negative correlation between FOSL1 and NF1 mRNA levels (Pearson R = À0.44, p=7.8e-12) ( Figure 2D and Supplementary file 4).
To test whether a NF1-MAPK signaling is involved in the regulation of FOSL1 and the MES subtype, we manipulated NF1 expression in patient-derived GBM tumorspheres of either MES or non-MES subtypes. To recapitulate the activity of the full-length NF1 protein, we transduced the cells with the NF1 GTPase-activating domain (NF1-GRD), spanning the whole predicted Ras GTPaseactivating (GAP) domain (McCormick, 1990). NF1-GRD expression in the MES cell line BTSC 233 led to (i) inhibition of RAS activity as confirmed by analysis of pERK expression upon EGF or serum stimulation (   Figure 2 continued on next page also observed a significant decrease of two well-characterized mesenchymal features, namely CHI3L1 expression ( Figure 2F) as well as the ability of MES GBM cells to differentiate into osteocytes, a feature shared with mesenchymal stem cells (Ricci-Vitiani et al., 2008;Tso et al., 2006;Figure 2-figure supplement 2D). Moreover, NF1-GRD expression led to a significant reduction of the FOSL1 regulon and the MGSs, with a concurrent increase of the OLIG2 regulon and the non-MES gene signatures (non-MGSs) ( Figure 2G).
Conversely, NF1 knockdown with three independent shRNAs (shNF1_1, shNF1_4, and shNF1_5) in two non-MES lines (BTSC 3021 and BTSC 3047) ( The observed NF1-mediated gene expression changes might be potentially driven by an effect on FOSL1 or other previously described mesenchymal TFs (such as BHLHB2, CEBPB, FOSL2, RUNX1, STAT3, and TAZ;) (Bhat et al., 2011;Carro et al., 2010). Interestingly, only FOSL1, and to some extent CEBPB, was consistently downregulated upon NF1-GRD expression (Figure 2-figure supplement 2H) and upregulated following NF1 knockdown (Figure 2-figure supplement 2I). To then test whether FOSL1 was playing a direct role in the NF1-mediated regulation of mesenchymal genes expression, we overexpressed FOSL1 in the MES GBM lines transduced with the NF1-GRD ( Figure 2-figure supplement 2J). qRT-PCR analysis showed that FOSL1 was able to rescue the NF1-GRD-mediated downregulation of mesenchymal genes, such as ITGA3, ITGA5, SERPINE1, and TNC ( Figure 3J). Lastly, exposure of NF1 silenced cells to the MEK inhibitor GDC-0623, led to a reduction of FOSL1 upregulation, both at the protein and the mRNA levels, as well as to a downregulation of the mesenchymal genes ITGA3 and SERPINE1 (Figure 2-figure supplement 3A, B).
Overall these evidences implicate the NF1-MAPK signaling in the regulation of the MGSs through the modulation of FOSL1 expression.

Fosl1 deletion induces a shift from a MES to a PN gene signature
To further explore the NF1-MAPK-FOSL1 axis in MES GBM, we used a combination of the RCAS-Tva system with the CRISPR/Cas9 technology, recently developed in our laboratory (Oldrini et al., 2018), to induce Nf1 loss or Kras mutation. Mouse NSCs from hGFAP-Tva; hGFAP-Cre; Trp53 lox ; ROSA26-LSL-Cas9 pups were isolated and infected with viruses produced by DF1 packaging cells transduced with RCAS vectors targeting the expression of Nf1 through shRNA and sgRNA (shNf1 and sgNf1) or overexpressing a mutant form of Kras (Kras G12V ). Loss of NF1 expression was confirmed by western blot, and FRA-1 was upregulated in the two models of Nf1 loss compared to Ctrl. (J) qRT-PCR analysis of MES genes expression upon NF1-GRD and FOSL1 co-expression in BTSC 232 and BTSC 233 cells. qRT-PCR data in (E), (H), and (J) are presented as mean ± SD (n = 3, technical replicates), normalized to 18S rRNA expression; Student's t test, *p 0.05, **p 0.01, ***p 0.001, ns = not significant. The online version of this article includes the following source data and figure supplement(s) for figure 2: Source data 1. Source data of Figure 2F. Source data 2. Source data of Figure 2A, C-E, G-J.             parental cells and further upregulated in cells infected with Kras G12V ( Figure 3A). Consistent with activation of the Ras signaling, as a result of both Nf1 loss and Kras mutation, the MEK/ERK pathway was more active in infected cells compared to parental cells ( Figure 3A). Higher levels of activation of the MEK/ERK pathway were more pronounced in the Kras mutant cells and were associated with a stronger induction of mesenchymal genes such as Plau, Plaur, Timp1, and Cd44 ( Figure 3B). (B) mRNA expression of Fosl1 and MES genes (Plau, Plaur, Timp1, and Cd44) in infected p53-null NSCs compared to parental cells (not infected). Data from a representative of two experiments are presented as mean ± SD (n = 3), normalized to Gapdh expression. Student's t test, relative to parental cells: ns = not significant, *p 0.05, **p 0.01, ***p 0.001. (C) FRA-1 expression detected by western blot in p53-null Kras G12V NSCs upon transduction with sgRNAs targeting Fosl1, after selection with 1 mg/mL puromycin; vinculin was used as loading control. (D) Gene set enrichment analysis (GSEA) results of p53-null Kras G12V sgFosl1_1 versus sgCtrl NSCs. (E, F) mRNA expression of MES (E) and PN genes (F) in sgCtrl and sgFosl1_1 p53-null Kras G12V NSCs. Data from a representative of two experiments are presented as mean ± SD (n = 3, technical replicates), normalized to Gapdh expression. Student's t test, relative to sgCtrl: *p 0.05; **p 0.01; ***p 0.001. The online version of this article includes the following source data for figure 3: Source data 1. Source data of Figure 3A. Source data 2. Source data of Figure 3C. Source data 3. Source data of Figure 3B Taking advantage of the Cas9 expression in the generated p53-null NSCs models, Fosl1 was knocked out through sgRNAs. Efficient downregulation of FRA-1 was achieved with two different sgRNAs ( Figure 3C and Figure 2-figure supplement 3E). Cells transduced with sgFosl1_1 and sgFosl1_3 were then subjected to further studies.
As suggested by the data presented here on the human BTSCs datasets and cell lines, FOSL1 appears to be a key regulator of the MES subtype. Consistently, RNA-seq analysis followed by GSEA of p53-null Kras G12V sgFosl1_1 versus sgCtrl revealed a significant loss of the MGSs and increase in the N-MGSs ( Figure 3D). These findings were validated by qRT-PCR with a significant decrease in expression of a panel of MES genes (Plau, Itga7, Timp1, Plaur, Fn1, Cyr61, Actn1, S100a4, Vim, Cd44) ( Figure 3E) and increased expression of PN genes (Olig2, Ncam1, Bcan, Lgr5) in the Fosl1 knock-out (KO) Kras G12V NSCs ( Figure 3F). A similar trend was observed in the Fosl1 KO shNf1 NSCs (Figure 2-figure supplement 3F, G), and the extent of MSG regulation appeared proportional to the extent of MAPK activation by individual perturbations ( Figure 3A).
Altogether, these data indicated that Kras G12V -transduced cells, which show the highest FOSL1 expression and mesenchymal commitment, are a suitable model to functionally study the role of a MAPK-FOSL1 axis in MES GBM.
Fosl1 depletion affects the chromatin accessibility of the mesenchymal transcription program and differentiation genes FOSL1 is a member of the AP-1 TF super family, which may be composed of a diverse set of homoand heterodimers of the individual members of the JUN, FOS, ATF, and MAF protein families. In GBM, AP-1 can act as a pioneer factor for other transcriptional regulators, such as ATF3, to coordinate response to stress in GSCs (Gargiulo et al., 2013). To test the effect of Fosl1 ablation on chromatin regulation, we performed open chromatin profiling using ATAC-seq in the p53-null Kras G12V NSCs model ( Figure 3C). This analysis revealed that Fosl1 loss strongly affects chromatin accessibility of known cis-regulatory elements such as transcription start sites (TSS) and CpG islands (CGI), as gauged by unsupervised clustering of Fosl1 wild-type and KO cells ( Figure 4A). Consistent with a role for FOSL1/FRA-1 in maintaining chromatin accessibility at direct target genes, deletion of Fosl1 caused the selective closing of chromatin associated with the major AP-1 TFs binding sites ( Figure 4B). Upon Fosl1 loss, profiling of the motifs indicated that chromatin associated with AP-1/2 TFs binding were closed and -conversely -a diverse set of general and lineage-specific TFs, including MFZ1, NRF1, RREB1, and others ( Figure 4C), were opened. The genes associated with changes in chromatin accessibility upon Fosl1 loss are involved in several cell fate commitment, differentiation, and morphogenesis programs ( Figure 4D, E). Next, we investigated chromatin remodeling dynamics using limma and identified 9749 regions with significant differential accessibility (absolute log2 fold-change >1, FDR < 0.05). Importantly, Fosl1 loss induced opening of chromatin associated with lineage-specific markers, along with closing of chromatin at the loci of genes, associated with mesenchymal GBM identity in human tumors and BTSC lines ( Figure 4F-H). Taken all together, this evidence further indicates that FOSL1/FRA-1 might modulate the mesenchymal transcriptional program by regulating the chromatin accessibility of MES genes.

Fosl1 deletion reduces stemness and tumor growth
Ras activating mutations have been widely used to study gliomagenesis, in combination with other alterations as Akt mutation, loss of Ink4a/Arf or Trp53 (Friedmann-Morvinski et al., 2012;Holland et al., 2000;Koschmann et al., 2016;Muñoz et al., 2013;Uhrbom et al., 2002). Thus, we then explored the possibility that Fosl1 could modulate the tumorigenic potential of the p53-null Kras mutant cells.
Cell viability was significantly decreased in Fosl1 KO cell lines as compared to sgCtrl ( Figure 5A). Concomitantly, we observed a significant decreased percentage of cells in S-phase (mean values: sgCtrl = 42.6%; sgFosl1_1 = 21.6%, Student's t test p 0.001; sgFosl1_3 = 20.4%, Student's t test p=0.003), an increase in percentage of cells in G2/M (mean values: sgCtrl = 11.7%, sgFosl1_1 = 28.4%, Student's t test p 0.001; sgFosl1_3 = 23.4%, Student's t test p=0.012) ( Figure 5B), and a Clustering of the Fosl1-WT (sgCtrl, n = 4) and Fosl1-depleted (sgFosl1_1 and sgFosl1_3, n = 8) samples is based upon the bias corrected deviations in chromatin accessibility (see Materials and methods). (B) tSNE visualization of cellular similarity between Fosl1-depleted and control cells based on chromatin accessibility. Samples are color-coded according to the cell type (black, red, and green for sgCtrl, sgFosl1_1, and sgFosl1_3 cells, respectively), or by directional z-scores. (C) Volcano plot illustrating the mean difference in bias-corrected accessibility deviations between Fosl1-deficient and control cells against the FDR-corrected p-value for that difference. The top differential motifs are highlighted in violet and red, indicating decreased and increased accessibility, respectively. (D, E) Top enriched Gene Ontology (GO) biological processes pathways for the regions with decreased (D) and increased (E) chromatin accessibility upon Fosl1 loss. The nodes represent the functional categories from the respective databases, color-coded by the significance of enrichment (FDR < 0.05). The node size indicates the number of query genes represented among the ontology term, and the edges highlight the relative relationships among these categories.  KO p53-null Kras G12V neural stem cells (NSCs) measured by MTT assay; absorbance values were normalized to day 1. Data from a representative of three independent experiments are presented as mean ± SD (n = 10, technical replicates). Two-way ANOVA, relative to sgCtrl for both sgFosl1_1 and sgFosl1_3: ***p 0.001. (B) Quantification of cell cycle populations of control and Fosl1 KO p53-null Kras G12V NSCs by flow cytometry analysis of PI staining. Data from a representative of two independent experiments are presented as mean ± SD (n = 3, technical replicates). Student's t test, relative to sgCtrl: *p 0.05; **p 0.01; ***p 0.001. (C) Representative limiting dilution experiment on p53-null Kras G12V sgCtrl and sgFosl1_1 NSCs, calculated with extreme limiting dilution assay (ELDA) analysis; bar plot inlet shows the estimated stem cell frequency with the confidence interval; chi-square p<0.0001. (D) Heatmap of expression of stem cell (yellow) and lineage-specific (neuronal -purple, astrocytic -green, and oligodendrocytic -orange) genes, comparing sgCtrl and sgFosl1_1 p53-null Kras G12V NSCs. Two biological replicates are shown. (E) Quantification of pixel area (fold-change relative to sgCtrl) of CD44, GFAP, and OLIG2 relative to DAPI pixel area per field of view in control and Fosl1 KO p53-null Kras G12V NSCs. Data from a representative of two independent experiments; Student's t test, relative to sgCtrl: ***p 0.001. (F) Kaplan-Meier survival Another aspect that contributes to GBM aggressiveness is its heterogeneity, attributable in part to the presence of GSCs. By using limiting dilution assays, we found that Fosl1 is required for the maintenance of stem cell activity, with a stem cell frequency estimate of sgFosl1_1 = 28.6 compared to sgCtrl = 3 (chi-square p<2.2e-16) ( Figure 5C). Moreover, RNA-seq analysis showed that sgFosl1_1 cells downregulated the expression of stem genes (Elf4,Klf4,Itgb1,Nes,Sall4,L1cam,Melk,Cd44,Myc,Fut4,Cxcr4,Prom1) while upregulating the expression of lineage-specific genes: neuronal (Map2, Ncam1, Tubb3, Slc1a2, Rbfox3, Dcx), astrocytic (Aldh1l1, Gfap, S100b, Slc1a3), and oligodendrocytic (Olig2, Sox10, Cnp, Mbp, Cspg4) ( Figure 5D). The different expression of some of the stem/differentiation markers was confirmed also by immunofluorescence analysis. While Fosl1 KO cells presented low expression of the stem cell marker CD44, differentiation markers as GFAP and OLIG2 were significantly higher when compared to sgCtrl cells ( Figure 5E and We then sought to test whether (i) p53-null Kras G12V NSCs were tumorigenic and (ii) Fosl1 played any role in their tumorigenic potential. Intracranial injections of p53-null Kras G12V NSCs in nu/nu mice led to the development of high-grade tumors with a median survival of 37 days in control cells (n = 9). In contrast, sgFosl1_1-injected mice (n = 6) had a significant increase in median survival (54.5 days, log-rank p=0.0263) ( Figure 5F). Consistent with our in vitro experiments ( Figure 3D-F), Fosl1depleted tumors failed to support mesenchymal program ( Figure 5G-I). Western blot, immunohistochemical, and qPCR analysis show the reduction of MES markers (VIM, CD44, and S100A4) and the expression of the PN marker OLIG2 to depend on Fosl1 expression ( Figure 5G-J).
Overall, our data in p53-null Kras mutant NSCs support the conclusion that, besides controlling cell proliferation, Fosl1 plays a critical role in the maintenance of the stem cell activity and tumorigenicity.

Fosl1 amplifies mesenchymal gene expression
To further assess the role of Fosl1 as a key player in the control of the MGS, we used a mouse model of inducible Fosl1 overexpression containing the alleles Kras LSLG12V ; Trp53 lox ; ROSA26 LSLrtTA-IRES-EGFP ; Col1a1  (here referred to as Fosl1 tetON ). Similar to the loss-of-function approach here used, this allelic combination allows the expression of Kras G12V and deletion of p53 after Cre recombination. Moreover, the expression of the reverse tetracycline transactivator (rtTA) allows, upon induction with doxycycline (Dox), the ectopic expression of Fosl1 (Flag tagged), under the control of the Col1a1 locus and a tetracycline-responsive element (TRE or Tet-O) (Belteki, 2005;Hasenfuss et al., 2014).
NSCs derived from Fosl1 WT and Fosl1 tetON mice were infected in vitro with a lentiviral vector expressing the Cre recombinase and efficient infection was confirmed by fluorescence microscopy as the cells expressing the rtTA should express GFP (data not shown). FRA-1 overexpression, as well as Flag-tag expression, was then tested by western blot after 72 hr of Dox induction ( Figure 6A). When Fosl1 tetON NSCs were analyzed by qRT-PCR for the expression of MES/PN markers, a significant upregulation of most MES genes and downregulation of PN genes was found in the cells overexpressing Fosl1 ( Figure 6B, C), thereby complementing our findings in Fosl1 KO cells. curves of nu/nu mice injected with p53-null Kras G12V sgCtrl (n = 9) and sgFosl1_1 (n = 6) NSCs. Log-rank p=0.0263. (G) Western blot analysis using the indicated antibodies of four sgCtrl and four sgFosl1_1 tumors (showing low or no detectable expression of FRA-1); vinculin was used as loading control. (H) Representative images of IHCs using the indicated antibodies. Scale bars represent 100 mm. (I) mRNA expression of MES genes in the samples sgCtrl-T4 (higher FRA-1 expression) and sgFosl1_1-T3 and -T4 (no detectable FRA-1 expression). (J) mRNA expression of PN genes in samples as in (H). Data from a representative of two experiments are presented as mean ± SD (n = 3, technical replicates), normalized to Gapdh expression. Student's t test for sgFosl1_1 tumors, relative to sgCtrl-T4: ns = not significant, *p 0.05, **p 0.01, ***p 0.001. The online version of this article includes the following source data and figure supplement(s) for figure 5: Source data 1. Source data of Figure 5G. Source data 2. Source data of Figure 5A-C, E, F, I, J.  Source data 1. Source data of Figure 6A. Source data 2. Source data of Figure 6E. Source data 3. Source data of Figure 6H. Source data 4. Source data of Figure 6B-D, F, G, I,J.  To investigate if the MES phenotype induced with Fosl1 overexpression would have any effect in vivo, p53-null Kras G12V Fosl1 tetON NSCs were intracranially injected into syngeneic C57BL/6J wildtype mice. Injected mice were randomized and subjected to Dox diet (food pellets and drinking water) or kept as controls with regular food and drinking water with 1% sucrose. Fosl1 overexpressing mice (+Dox) developed larger tumors that were more infiltrative and aggressive than controls (-Dox), which mostly grew as superficial tumor masses instead ( Figure 6D). This phenotype appears to be independent of tumor cells proliferation as gauged by Ki-67 staining and does not affect overall survival ( Figure 6-figure supplement 1A, B).
Tumorspheres were derived from -Dox and +Dox tumor-bearing mice, and Fosl1 expression was manipulated in vitro through addition or withdrawal of Dox from the culture medium. In the case of tumorspheres derived from a -Dox tumor, when Dox was added for 19 days, high levels of FRA-1 expression were detected by western blot (Figure 6E). At the mRNA level, Dox treatment also greatly increased Fosl1 expression, as well as some of the MES genes ( Figure 6F), while the expression of PN genes was downregulated ( Figure 6G). Conversely, when Dox was removed from +Doxderived tumorspheres for 19 days, the expression of FRA-1 decreased ( Figure 6H, I), along with the expression of MES genes ( Figure 6I), while PN genes were upregulated ( Figure 6J). These results confirm the essential role of Fosl1 in the regulation of the MGS in p53-null Kras G12V tumor cells and the plasticity between the PN and MES subtypes.

FOSL1 controls growth, stemness, and mesenchymal gene expression in patient-derived BTSCs
To prove the relevance of our findings in the context of human tumors, we analyzed BTSC lines characterized as non-MES (h676, h543, and BTSC 268) or MES (BTSC 349, BTSC 380, and BTSC 233) (this study and Ozawa et al., 2014). By western blot, we found that consistent with what was observed either in human BTSCs ( Figure 1D) or mouse NSCs ( Figure 3A), MES cell lines expressed high levels of FRA-1 and activation of the MEK/ERK pathway ( Figure 7A).
We then tested whether FOSL1/FRA-1 modulates the MGS via direct target regulation. To this end, we first identified high-confidence FOSL1/FRA-1 binding sites in chromatin immunoprecipitation-seq (ChIP-seq) previously generated in the KRAS mutant HCT116 colorectal cancer cell line (see Materials and methods), and then we determined the counts per million reads (CPM) of the enhancer histone mark H3K27Ac in a set of MES (n = 10) and non-MES BTSCs (n = 10) (Mack et al., 2019), selected based on the highest and lowest FOSL1 expression, respectively. PCA showed a marked separation of the two groups of BTSCs ( Figure 7F). Differential enrichment analysis by DESeq2 revealed 11748 regions statistically significant (FDR < 0.005) for H3K27Ac at FOSL1/FRA-1 binding sites in either MES or non-MES BTSCs ( Figure 7G). Next, we compared H3K27Ac distribution over FOSL1/FRA-1 binding sites to that of the non-MES MR OLIG2. This analysis showed that  Figure 7J).
Altogether, our data support that FOSL1/FRA-1 regulates MES gene expression and aggressiveness in human gliomas via direct transcriptional regulation, downstream of the NF1-MAPK-FOSL1 signaling.

Discussion
The most broadly accepted transcriptional classification of GBM was originally based on gene expression profiles of bulk tumors (Verhaak et al., 2010), which did not discriminate the contribution of tumor cells and TME to the transcriptional signatures. It is now becoming evident that both cell-intrinsic and -extrinsic cues can contribute to the specification of the MES subtype (Bhat et al., 2013;Hara et al., 2021;Neftel et al., 2019;Schmitt et al., 2021;Wang et al., 2017). Bhat and colleagues had shown that while some of the MES GBMs maintained the mesenchymal characteristics when expanded in vitro as BTSCs, some others lost the MGS after few passages while exhibiting a higher non-MGSs (Bhat et al., 2013). These data, together with the evidence that xenografts into immunocompromised mice of BTSCs derived from MES GBMs were also unable to fully restore the MES phenotype (Bhat et al., 2013), suggested that the presence of an intact TME potentially contributed to the maintenance of a MGS. In support of this, Schmitt and colleagues have recently shown that innate immune cells divert GBM cells to a proneural-to-mesenchymal transition (PN-to-MES) that also contributes to therapeutic resistance .
The transcriptional GBM subtypes were lately redefined based on the expression of glioma-intrinsic genes, thus excluding the genes expressed by cells of the TME (Richards et al., 2021;Wang et al., 2017). Our MRA on the BTSCs points to the AP-1 family member FOSL1 as one pioneer TF contributing to the cell-intrinsic MGS. Previous tumor bulk analysis identified a related AP-1 family member FOSL2, together with CEBPB, STAT3, and TAZ, as important regulators of the MES GBM subtype (Bhat et al., 2011;Carro et al., 2010). While FOSL1 was also listed as a putative MES MR (Carro et al., 2010), its function and mechanism of action have not been further characterized since then. Our experimental data show that FOSL1 is a key regulator of GBM subtype plasticity and MES transition, and define the molecular mechanism through which FOSL1 is regulated. While here Source data 1. Source data of Figure 7A. Source data 2. Source data of Figure 7B. Source data 3. Source data of Figure 7C-G, J.     we have focused on the TFs contributing to MES specifications, previous studies had highlighted the role of other TFs, some of which were also identified in our MRA, such as OLIG2, SALL2, and ASCL1, as important molecules for non-MES GBM cells . Moreover, using a similar MRA, Wu and colleagues have recently described also SOX10 as another TF that contributes to the identity of non-MES GBM cells. Strikingly, loss of SOX10 resulted in MES transition associated with changes in chromatin accessibility in regions that are specifically enriched for FRA-1 binding motifs (Wu et al., 2020). Lastly, using an unbiased CRISPR/Cas9 genome-wide screening, Richards and colleagues had shown that few of the top TFs identified here, such as FOSL1, OLIG2, and ASCL1, are genes essential specifically either for MES GSCs (FOSL1) or for non-MES GSCs (OLIG2 and ASCL1) (Richards et al., 2021). This evidence further strengthen the relevance of the MRA that we have performed in the identification of important regulators of GBM subtype-specific cell biology.
Although consistently defined, GBM subtypes do not represent static entities. The plasticity between subtypes happens at several levels. Besides the referred MES-to-PN change in cultured GSCs compared to the parental tumor (Bhat et al., 2013), a PN-to-MES shift often occurs upon treatment and recurrence. Several independent studies comparing matched pairs of primary and recurrent tumors demonstrated a tendency to shift towards a MES phenotype, associated with a worse patient survival, likely as a result of treatment-induced changes in the tumor and/or the microenvironment (Phillips et al., 2006;Varn et al., 2021;Wang et al., 2016;Wang et al., 2017). Moreover, distinct subtypes/cellular states can coexist within the same tumor (Neftel et al., 2019;Patel et al., 2014;Richards et al., 2021;Sottoriva et al., 2013;Varn et al., 2021;Wang et al., 2019) and targeting these multiple cellular components could result in more effective treatments .
PN-to-MES transition is often considered an EMT-like phenomenon, associated with tumor progression (Fedele et al., 2019). The role of FOSL1 in EMT has been studied in other tumor types. In breast cancer cells, FOSL1 expression correlates with mesenchymal features and drives cancer stem cells  and the regulation of EMT seems to happen through the direct binding of FRA-1 to promoters of EMT genes such as Tgfb1, Zeb1, and Zeb2 (Bakiri et al., 2015). In colorectal cancer cells, FOSL1 was also shown to promote cancer aggressiveness through EMT by direct transcription regulation of EMT-related genes (Diesch et al., 2014;Liu et al., 2015).
It is well established that NF1 inactivation is a major genetic event associated with the MES subtype (Verhaak et al., 2010;Wang et al., 2017). However, this is probably a late event in MES gliomagenesis as all tumors possibly arise from a PN precursor and just later in disease progression acquire NF1 alterations that are directly associated with a transition to a MES subtype . Moreover, NF1 deficiency has been linked to macrophage/microglia infiltration in the MES subtype . The fact that the enriched macrophage/microglia microenvironment is also able to modulate a MES phenotype suggests that there might be a twoway interaction between tumor cells and TME. The mechanisms of NF1-regulated chemotaxis and whether this relationship between the TME and MGS in GBM is causal remain elusive.
Here, we provide evidence that manipulation of NF1 expression levels in patient-derived BTSCs has a direct consequence on the tumor-intrinsic MGS activation and that such activation can at least in part be mediated by the modulation of FOSL1. Among the previously validated MRs, only CEBPB appears also to be finely tuned by NF1 inactivation. This suggests that among the TFs previously characterized (such as FOSL2, STAT3, BHLHB2, and RUNX1), FOSL1 and CEBPB might play a dominant role in the NF1-mediated MES transition that occurs in a glioma cell-intrinsic manner. However, whether FOSL1 contributes also to the cross-talk between the TME and the cell-intrinsic MGS still has to be established.
Furthermore, we show that FOSL1 is a crucial player in glioma pathogenesis, particularly in a MAPK-driven MES GBM context (Figure 8). Most likely, the existence of a NF1-MAPK-FOSL1 axis goes beyond GBM pathogenesis since FOSL1 appears to be upregulated in concomitance with NF1 mutations in multiple tumor types (Figure 8-figure supplement 1). Our findings broaden its previously described role in KRAS-driven epithelial tumors, such as lung and pancreatic ductal adenocarcinoma (Vallejo et al., 2017). NF1 inactivation results in Ras activation, which stimulates downstream pathways as MAPK and PI3K/Akt/mTOR. RAS/MEK/ERK signaling in turn regulates FRA-1 protein stability (Casalino et al., 2003;Verde et al., 2007). FOSL1 mRNA expression is then most likely induced by binding of the SRF/Elk1 complex to the serum-responsive element (SRE) on FOSL1 promoter (Esnault et al., 2017). At the same time, FRA-1 can then directly bind to its own promoter to activate its own expression (Diesch et al., 2014;Lau et al., 2016) and those of MES genes. This further generates a feedback loop that induces MGS, increases proliferation and stemness, sustaining tumor growth. FRA-1 requires, for its transcriptional activity, heterodimerization with the AP-1 TFs JUN, JUNB, or JUND (Eferl and Wagner, 2003). Which of the JUN family members participate in the MES gene regulation and whether FOSL1/FRA-1 activates MES gene expression and simultaneously represses non-MES genes requires further investigation. Of note, pancancer analysis of anatomically distinct solid tumors suggested that c-JUN/JUNB and FOSL1/2 are bona fide canonical AP-1 TF configurations in mesenchymal states of lung, kidney, and stomach cancers . Intriguingly, in support of a direct role in the repression of non-MES genes in GBM cells, it has been hypothesized, though not formally demonstrated, that FOSL1/FRA-1 could act as a transcriptional repressor of a core set of neurodevelopmental TFs, including OLIG2 and SALL2 (Fiscon et al., 2018).
In conclusion, FOSL1/FRA-1 is a key regulator of the MES subtype of GBM, significantly contributing to its stem cell features, which could open new therapeutic options. Although FOSL1/FRA-1 pharmacological inhibition is difficult to achieve due to its enzymatic activity, a gene therapy approach targeting FOSL1/FRA-1 expression through CRISPR/Cas9 or PROTAC, for instance, could constitute attractive alternatives to treat mesenchymal GBM patients.
To analyze the RNA-seq samples, we used the Nextpresso pipeline (Graña et al., 2018). For the Affymetrix microarrays, raw data were downloaded from the GEO repository (https://www.ncbi.nlm. nih.gov/geo/) and subsequently the 'affy' R package (Gautier et al., 2004) was used for robust multi-array average normalization followed by quantile normalization. For genes with several probe sets, the median of all probes had been chosen and only common genes among all the datasets (n = 9889) were used for further analysis. To avoid issues with the use of different transcriptomic platforms, each dataset was then scaled (mean = 0, SD = 1) before assembling the combined final dataset. Transcriptional subtypes were obtained using the 'ssgsea.GBM.classification' R package , using 1000 permutations. For differential gene expression studies, we selected the 133 BTSCs lines that had concordant ssGSEA results, with MES BTSCs classified both as MES and INJ and non-MES BTSCs classified both as PN/CL and DEV. Differential gene expression (MES vs. non-MES BTSCs) was performed using the 'limma' R package (Ritchie et al., 2015), taking into account the possible batch differences due the different datasets assembled. The MRA was performed using the 'RTN' R package (Castro et al., 2016). Normalized BTSC expression data were used as input to build a transcriptional network (TN) for 887 TFs present in the dataset. TF annotations were obtained from the human TF atlas version 1.0.1 (http://humantfs.ccbr. utoronto.ca/) (Lambert et al., 2018). p-Values for network edges were computed from a pooled null distribution using 1000 permutations. Edges with an adjusted-p-value<0.05 were kept for data processing inequality (DPI) filtering. In the TN, each target can be connected to multiple TFs and regulation can occur as a result of both direct and indirect interactions. DPI filtering removes the weakest interaction in any triangle of two TFs and a target gene, therefore preserving the dominant TF-target pairs and resulting in a filtered TN that highlights the most significant interactions (Fletcher et al., 2013). Post-DPI filtering, the MRA computes the overlap between the transcriptional regulatory unities (regulons) and the input signature genes using the hypergeometric distribution (with multiple hypothesis testing corrections). To identify MRs, the differential gene expression between MES and non-MES was used as a phenotype.
TCGA and CGGA data analysis RSEM normalized RNA-seq data for the TCGA GBMLGG and CGGA datasets were downloaded from the Broad Institute Firebrowse (http://gdac.broadinstitute.org) and the Chinese Glioma Genome Atlas (updated November 2019) (http://www.cgga.org.cn/), respectively. NF1 copy number alterations and point mutations for the TCGA GBMLGG were obtained at the cBioPortal (https:// www.cbioportal.org). FOSL1 expression and NF1 mutational status for the TCGA datasets in Figure 8-figure supplement 1 were obtained from the Timer2.0 web portal (http://timer.cistrome. org/) (Li et al., 2020). Transcriptional subtypes were inferred using the 'ssgsea.GBM.classification' R package as indicated above. Glioma molecular subtypes information was downloaded from the Glio-Vis web portal (http://gliovis.bioinfo.cnio.es) (Bowman et al., 2017). Survival analysis was performed using the 'survival' R package. Stratification of the patients has been done by comparing the patients with the 30% FOSL1 high expression to the 30% FOSL1 low expression.

Gene expression array and GSEA
For gene expression profiling of the BTSC lines of the Freiburg dataset, total RNA was prepared using the RNeasy kit (QIAGEN #74104) or the AllPrep DNA/RNA/Protein mini kit (QIAGEN #80004) and quantified using 2100 Bioanalyzer (Agilent). One-and-a-half microgram of total RNA for each sample was sent to the genomic facility of the German Cancer Research Center (DKFZ) in Heidelberg (Germany), where hybridization and data normalization were performed. Hybridization was carried out on Illumina HumanHT-12v4 expression BeadChip. GSEA was performed using the GSEA software (http://www.broadinstitute.org/gsea/index.jsp) (Subramanian et al., 2005). Gene signatures are listed in Supplementary file 5.

ATAC-seq
ATAC-seq was performed on 60,000 p53-null Kras G12V mouse NSCs transduced with either sgFosl1_1, sgFosl1_3, or sgCtrl. Briefly, the cells were pelleted in PBS and tagmentation was performed in either 50 mL of master mix containing 25 mL 2xTD buffer, 2.5 mL transposase, and 22.5 mL nuclease-free water (Nextera DNA Library Prep, Illumina #FC-121-1030) or in 50 mL of tagmentation mix containing 25 mL of TAPS-DMF buffer (80 mM TAPS, 40 mM MgCl 2 , 50% vol/vol DMF), 0.625 mL in-house-produced hyperactive Tn5 enzyme and 24.4 mL nuclease-free water (adapted from Hennig et al., 2018). The tagmentation reactions were incubated for 1 hr at 37˚C with moderate agitation (500-800 rpm). After the incubation, 5 mL of Proteinase K (Invitrogen #AM2548) were added to the samples to stop the transposition. The cells were subsequently lysed by adding 50 mL of AL buffer (QIAGEN #19075) and incubating for 10 min at 56˚C. The DNA was extracted by means of 1.8Â vol/vol AMPure XP beads (Beckman Coulter #A63881) and eluted in 18 mL. To determine the optimal number of PCR cycles required for library amplification, 2 mL of each sample were taken as template for qPCR using KAPA HiFi HotStart ReadyMix (Roche #7958927001) and 1xEvaGreen Dye (Biotium #31000). The whole probe amplification was performed in 50 mL qPCR volume with 8-12 mL of template DNA. Primers were previously described (Buenrostro et al., 2013). Each library was individually quantified utilizing the Qubit 3.0 Fluorometer (Life Technologies). The appropriate ATAC-seq libraries laddering pattern was determined with TapeStation High Sensitivity D1000 ScreenTapes (Agilent #5067-5584). The libraries were sequenced on the Illumina NextSeq500 using the High Output V2 150 cycles chemistry kit in a 2 Â 75 bp mode.

ATAC-seq analysis
The reads were adaptor-trimmed using the trim-galore v0.6.2 -nextera (https://www.bioinformatics.babraham.ac.uk/projects/trim_galore/). The mapping was conducted using the bowtie2 v2.3.5 (Langmead and Salzberg, 2012) default parameters. The differential chromatin accessibility analysis was performed in SeqMonk taking mouse CpG islands (mCGI) and TSSs ± 500 bp as the probe set (GRCm38 assembly). Counts were normalized by means of the 'Read Count Quantification' function with additional correction for total count (CPM), log transformation, and size factor normalization. The differential accessibility between the Fosl1-WT and Fosl1-depleted cells was determined utilizing the limma pipeline (Ritchie et al., 2015).
The 'chromVAR' R package (Schep et al., 2017) was used to analyze the chromatin accessibility data, perform corrections for known technical biases, and identify motifs with differential deviation in chromatin accessibility between the samples. The enrichGO function from the 'clusterProfiler' R package  was used to visualize the relevant pathways in Figure 4D, E. The bamCoverage function of the deepTools2 tool (Ramírez et al., 2016) was used to generate BigWig from aligned files for subsequent visualization with the 'karyoploteR' R package (Gel and Serra, 2017).

Vectors, virus production, and infection
Flag-tagged NF1-GRD (aminoacids 1131-1534) was amplified by PCR from human cortical tissue (epilepsy patient) and first cloned in the pDRIVE vector (QIAGEN #231124). Primers are listed in Supplementary file 6. The NF1-GRD sequence was then excised by restriction digestion using PmeI and SpeI enzymes and subcloned in the modified pCHMWS lentiviral vector (kind gift from V. Baekelandt, University of Leuven, Belgium) sites by removing the fLUC region. The correct sequence was verified by sequencing. For NF1 silencing, shNF1 from pLKO (Sigma, TRCN0000238778) (shNF1_1) was subcloned in pGIPZ lentiviral vector (Open Biosystems). The corresponding short hairpin sequence was synthetized (GATC) and amplified by PCR using XhoI and EcoRI sites containing primers. The PCR product was digested using XhoI and EcoRI and subcloned into the pGIPZ vector previously digested with XhoI and PmeI followed by digestion with EcoRI. The two vector fragments were ligated with NF1 short hairpin fragment. The correct insertion and sequence was validated by sequencing. In addition, experiments were performed using shNF1-pGIPZ clone V2LHS_76027 (shNF1_4) and V2LHS_260806 (shNF1_5).
Gp2-293 packaging cell line (Clontech #631458) was grown in DMEM (Sigma-Aldrich #D5796) with 10% FBS. Lentiviruses generated in this cell line were produced using calcium-phosphate precipitate transfection and co-transfected with second-generation packaging vectors (pMD2G and psPAX2). High-titer virus was collected at 36 and 60 hr following transfection.
All cells were infected with lenti-or retroviruses by four cycles of spin infection (200 Â g for 2 hr) in the presence of 8 mg/mL polybrene (Sigma-Aldrich #H9268). Transduced cells were selected after 48 hr from the last infection with 1 mg/mL puromycin (Sigma-Aldrich #P8833).
Mice were anesthetized with 4% isofluorane and then injected with a stereotactic apparatus (Stoelting) as previously described (Hambardzumyan et al., 2009). After intracranial injection, all mice were routinely checked and sacrificed when developed symptoms of disease (lethargy, poor grooming, weight loss, and macrocephaly).

Immunohistochemistry
Tissue samples were fixed in 10% formalin, paraffin-embedded, and cut in 3 mm sections, which were mounted in Superfrost Plus microscope slides (Thermo Scientific #J1810AMNZ) and dried. Tissues were deparaffinized in xylene and re-hydrated through graded concentrations of ethanol in water, ending in a final rinse in water.

Reverse transcription quantitative PCR
RNA from NSCs and frozen tissue was isolated with TRIzol reagent (Invitrogen #15596-026) according to the manufacturer's instructions. For reverse transcription PCR (RT-PCR), 500 ng of total RNA was reverse transcribed using the High Capacity cDNA Reverse Transcription Kit (Applied Biosystems #4368814). Quantitative PCR was performed using the SYBR Select Master Mix (Applied Biosystems #4472908) according to the manufacturer's instructions. qPCRs were run and the melting curves of the amplified products were used to determine the specificity of the amplification. The threshold cycle number for the genes analyzed was normalized to GAPDH. Mouse and human primer sequences are listed in Supplementary file 6.
RNA from BTSC cells was prepared using the RNeasy kit or the AllPrep DNA/RNA Protein Mini Kit and used for first-strand cDNA synthesis using random primers and SuperscriptIII reverse transcriptase (Life Technologies #18080-085). Primer sequences used in qRT-PCR with SYBR Green are listed in Supplementary file 6. Quantitative real-time PCR (qRT-PCR) STAT3 and CEBPB were performed using pre-validated TaqMan assays (Applied Biosystems): STAT3: Hs01047580, CEBPB: Hs00270923 and 18S rRNA: Hs99999901.

MTT assay
Cells were seeded in 96-well culture plates (1000 cells per well, 10 wells per cell line) and grown for 7 days. At each timepoint (days 1, 3, 5, and 7), cell viability was determined by MTT assay. Briefly, 10 mL of 5 mg/mL MTT (Sigma-Aldrich #M5655) was added to each well and cells were incubated for 4 hr before lysing with a formazan solubilization solution (10% SDS in 0.01 M HCl). Colorimetric intensity was quantified using a plate reader at 590 nm. Values were obtained after subtraction of matched blanks (medium only).
Cell cycle analysis: propidium iodide (PI) staining Cells were harvested and washed twice with PBS prior to fixation with 70% cold ethanol, added drop-wise to the cell pellet while vortexing. Fixed cells were then washed, first with 1% BSA in PBS, then with PBS only and stained overnight with 50 mg/mL PI (Sigma-Aldrich #P4170) and 100 mg/mL RNase A (Roche #10109142001) in PBS. Samples were acquired in a FACSCanto II cytometer (BD Biosciences), and data were analyzed using FlowJo software.

BrdU and EdU incorporation assays
Cells were pulse-labeled with 10 mM BrdU (Sigma-Aldrich #B9285) for 2 hr, harvested and washed twice with PBS prior to fixation with 70% cold ethanol, and added drop-wise to the cell pellet while vortexing. DNA denaturation was performed by incubating samples for 10 min on ice with 0.1 M HCl with 0.5% Tween-20, and samples were then resuspended in water and boiled at 100˚C for 10 min. Anti-BrdU-FITC antibody (BD Biosciences #556028) was incubated according to the manufacturer's protocol. After washing with PBSTB (PBS with 0.5% Tween-20% and 1% BSA), samples were resuspended in 25 mg/mL PI and 100 mg/mL RNase A diluted in PBS. Samples were acquired in a FACSCanto II cytometer (BD Biosciences), and data were analyzed using FlowJo software.
EdU incorporation was assessed using the EdU-Click594 Cell Proliferation Imaging Kit (Baseclick GmbH) according to the manufacturer's instructions. 96 hr after transduction, 2.0 Â 10 4 BTSC 233 cells were seeded on laminin-coated glass coverslips in a 24-well cell culture plate. Pictures were acquired using an Axiovert Microscope (Zeiss).
Fluorescence signal was quantified as the ratio of green/red pixel area relative to DAPI pixel area per field of view in a total of 36 fields per condition analyzed.

Neurosphere formation assay and limiting dilution analysis
Neurospheres were dissociated and passed through a 40 mm mesh filter to eliminate non-single cells. Decreasing cell densities were plated in ultra-low attachment 96-well plates (Corning #CLS3474), and fresh medium was added every 3-4 days. The number of positive wells for the presence of spheres was counted 2 weeks after plating. Limiting dilution analysis was performed using ELDA R package (http://bioinf.wehi.edu.au/software/elda/).

RNA-sequencing and analysis on mouse NSCs
For the p53-null Kras G12V NSCs, 1 mg of total RNA from the samples was used. cDNA libraries were prepared using the 'QuantSeq 3 'mRNA-Seq Library Prep Kit (FWD) for Illumina' (Lexogen #015) by following the manufacturer's instructions. Library generation is initiated by reverse transcription with oligo(dT) priming, and a second-strand synthesis is performed from random primers by a DNA polymerase. Primers from both steps contain Illumina-compatible sequences. Adapter-ligated libraries were completed by PCR, applied to an Illumina flow cell for cluster generation, and sequenced on an Illumina HiSeq 2500 by following the manufacturer's protocols. Sequencing read alignment and quantification and differential gene expression analysis was performed in the Bluebee Genomics Platform, a cloud-based service provider (https://www.illumina.com/company/about-us/mergersacquisitions/bluebee.html). Briefly, reads were first trimmed with bbduk from BBTools (BBMap -Bushnell B, https://sourceforge.net/projects/bbmap/) to remove adapter sequences and polyA tails. Trimmed reads were aligned to the GRCm38/mm10 genome assembly with STAR v 2.5 (Dobin et al., 2013). Read counting was performed with HTSeq (Anders et al., 2015). The list of stem/differentiation markers was compiled by combining a previously described gene list (Sandberg et al., 2013) with other markers (Bazzoli et al., 2012).
For the p53-null shNf1 NSCs, total RNA samples (500 ng) were converted into sequencing libraries with the 'NEBNext Ultra II Directional RNA Library Prep Kit for Illumina' (NEB #E7760), as recommended by the manufacturer. Briefly, polyA+ fraction is purified and randomly fragmented, converted to double-stranded cDNA, and processed through subsequent enzymatic treatments of end-repair, dA-tailing, and ligation to adapters. Adapter-ligated library is completed by PCR with Illumina PE primers. The resulting purified cDNA libraries were applied to an Illumina flow cell for cluster generation and sequenced on an Illumina NextSeq 550 by following the manufacturer's protocols. We then used the Nextpresso pipeline (Graña et al., 2018) for alignment and quantification.

Osteogenesis differentiation assay
The osteogenesis differentiation assay was performed using the StemPro Osteogenesis Differentiation Kit (Life Technologies #A1007201) according to the manufacturer's instructions. Briefly, 5 Â 10 3 cells/cm 2 were seeded on laminin-coated glass coverslips in a 24-well cell culture plate. Cells were incubated in MSC Growth Medium at 37˚C, 5% CO 2 for 21 days, replacing the medium every 4 days. Cells were then fixed with 4% formaldehyde, stained with Alizarin Red S solution (pH 4.2), and mounted on microscope slides. Pictures were acquired using an Axiovert Microscope (Zeiss).

Active Ras pull-down assay
Active Ras pull-down assay was performed using Active Ras pull-down assay kit (Thermo Fisher Scientific #16117) according to the manufacturer's instructions. Briefly, cells were grown in 10 cm plates at 80-90% confluency in the presence or absence of growth factors (EGF, FGF, and LIF) and lysed with the provided buffer. Lysates were incubated with either GDP or GTP for 30 min followed by precipitation with GST-Raf1-RBD. Western blot was performed with the provided anti-RAS antibody (1:250).

Statistical analysis
All statistical analyses were performed using R programming language (3.6.3). Statistical differences between groups were assessed by one-way ANOVA, two-way ANOVA, or unpaired two-tailed Student's t tests unless otherwise specified.
Kaplan-Meier survival curves were produced with GraphPad Prism, and p-values were generated using the log-rank statistics.
Results are presented as mean ± standard deviation (SD), and statistical significance was defined as p 0.05 for a 95% confidence interval.