# Cellular senescence in malignant cells promotes tumor progression in mouse and patient Glioblastoma

Rana Salam, Alexa Saliou, Franck Bielle, Mathilde Bertrand, Christophe Antoniewski, Catherine Carpentier, Agusti Alentorn, Laurent Capelle, Marc Sanson, Emmanuelle Huillard, et al.

► To cite this version:

Rana Salam, Alexa Saliou, Franck Bielle, Mathilde Bertrand, Christophe Antoniewski, et al. Cellular senescence in malignant cells promotes tumor progression in mouse and patient Glioblastoma. Nature Communications, 2023, 14 (1), pp.441. 10.1038/s41467-023-36124-9. hal-03969354v2

HAL Id: hal-03969354

https://hal.sorbonne-universite.fr/hal-03969354v2

Submitted on 12 Sep 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

# Cellular senescence in malignant cells promotes tumor progression in mouse and patient Glioblastoma

Received: 29 April 2022

Accepted: 16 January 2023

Published online: 27 January 2023

<figure><img src="image_3.png" /></figure>

Check for updates

Rana Salam ©<fnref n="1" />,⁸, Alexa Saliou<fnref n="1" />,⁸, Franck Bielle©<fnref n="1" />,²,³, Mathilde Bertrand⁴, Christophe Antoniewski⁵, Catherine Carpentier<fnref n="1" />, Agusti Alentorn©<fnref n="1" />,⁶, Laurent Capelle⁷, Marc Sanson<fnref n="1" />,³,⁶, Emmanuelle Huillard<fnref n="1" />, Léa Bellenger⁵, Justine Guégan©⁴ & Isabelle Le Roux©<fnref n="1" />✉

Glioblastoma (GBM) is the most common primary malignant brain tumor in adults, yet it remains refractory to systemic therapy. Elimination of senescent cells has emerged as a promising new treatment approach against cancer. Here, we investigated the contribution of senescent cells to GBM progression. Senescent cells are identified in patient and mouse GBMs. Partial removal of p16INK4a-expressing malignant senescent cells, which make up less than 7 % of the tumor, modifies the tumor ecosystem and improves the survival of GBM-bearing female mice. By combining single cell and bulk RNA sequencing, immunohistochemistry and genetic knockdowns, we identify the NRF2 transcription factor as a determinant of the senescent phenotype. Remarkably, our mouse senescent transcriptional signature and underlying mechanisms of senescence are conserved in patient GBMs, in whom higher senescence scores correlate with shorter survival times. These findings suggest that senolytic drug therapy may be a beneficial adjuvant therapy for patients with GBM.

Diffuse gliomas are the most common primary malignant brain tumors in adults<fnref n="1" />. Glioblastoma (GBM; IDH-wild type glioma, grade 4) is the most aggressive glioma and despite intensive conventional therapy which includes surgery, radiation, and both concurrent and adjuvant temozolomide (TMZ) chemotherapy, GBM remains treatment-resistant and disease progression is fatal, with a median survival below 15 months². Distinct factors may account for current treatment failures, including tumor invasiveness, an immunosuppressive microenvironment and intratumoral heterogeneity. Novel approaches are therefore required to find effective therapeutic strategies.

Cellular senescence is a permanent cell cycle arrest mediated by p53/p21<sup>CIP1</sup> and/or p16<sup>INK4A</sup>/Rb pathways and is defined by a combination of features including a senescence-associated secretory phenotype (SASP), anti-apoptotic program, and increased lysosomal content, the latter allowing histochemical detection of senescence associated-β-galactosidase activity (SA-β-gal)³. In cancer, cellular senescence is triggered by multiple stresses such as DNA damage, oncogene activation, therapeutic agents, or elevated reactive oxygen species (ROS). SASP is defined by the secretion of a plethora of factors including cytokines, chemokines, growth factors, extracellular matrix (ECM) components, and proteases, which together can stimulate

<footnote n="1">¹Paris Brain Institute (ICM), Hôpital Pitié-Salpêtrière, Inserm U 1127, CNRS UMR 7225, Sorbonne Université, Genetics and Development of Brain Tumors Team, Paris, France. ²AP-HP, Hôpital de la Pitié-Salpêtrière—Charles Foix, Département de Neuropathologie, Paris, France. ³Paris Brain Institute (ICM), Hôpital Pitié-Salpêtrière, Inserm U 1127, CNRS UMR 7225, Sorbonne Université, Onconeurotek Tumor Bank, Paris, France. ⁴Paris Brain Institute (ICM), Hôpital Pitié-Salpêtrière, Inserm U 1127, CNRS UMR 7225, Sorbonne Université, Data Analysis Core Platform, Paris, France. ⁵Sorbonne Université, CNRS FR3631, Inserm USO37, Institut de Biologie Paris Seine (IBPS), ARTbio Bioinformatics Analysis Facility, Paris, Institut Français de Bioinformatique (IFB), Paris, France. ⁶AP-HP, Hôpital de la Pitié-Salpêtrière—Charles Foix, Service de Neurologie 2-Mazarin, Paris, France. ⁷AP-HP, Hôpital de la Pitié-Salpêtrière—Charles Foix, Service de Neurochirurgie, Paris, France. ⁸These authors contributed equally: Rana Salam, Alexa Saliou. ✉e-mail: isabelle.leroux@icm-institute.org</footnote>
angiogenesis, modulate the composition of the ECM and promote an epithelial-to-mesenchymal transition⁴,⁵. Depending on the context, senescence exerts two opposite effects during tumorigenesis⁶. In some contexts, senescent cells prevent the proliferation of pre-malignant cancer cells, as SASP factors stimulate the immune clearance of onco-gene or therapy-induced senescent tumor cells⁷⁻⁹. Conversely, in persistently senescent cells, the SASP can either directly induce tumor growth¹⁰ or contribute to immune suppression, thus allowing tumor progression¹¹,¹². Many studies have assessed the function of senescence in developing tissues and age-related diseases by the in vivo removal of senescent cells, either using chemical or genetic senolytics¹³⁻¹⁷. A common genetic approach employs *p16*<sup>INK4a</sup> regulatory sequences to drive the inducible expression of INK-ATTAC or p16-3MR, which selectively eliminate senescent cells expressing *p16*<sup>INK4a</sup>, leading to apoptosis¹⁶,¹⁷. This senolytic strategy efficiently reduces the adverse effects of therapy-induced senescent cells including cancer recurrence in a mouse breast cancer model¹⁸.

A few in vivo studies have begun to examine the role of cellular senescence in gliomas. Using mouse patient-derived xenograft (PDX) models, it was shown that IL6, a universal SASP component as well as a cytokine expressed by immune cells, promotes the growth of patient glioma stem cells (GSCs) and contributes to glioma malignancy¹⁹. Conversely, loss of PTEN-PRMT5 signaling induces senescent GSCs to slow down tumorigenesis²⁰. Furthermore, loss of one allele of p53 reduces H-RasV12 oncogene-induced senescence in an orthotopic GBM model, as evidenced by reduced SA-β-gal staining, and decreases mouse survival time²¹. Finally, a recent study revealed the dual effect of therapy-induced senescence (TIS) following BMI1 inhibitor treatment of diffuse intrinsic pontine glioma tumor (a pediatric high-grade glioma), which initially attenuates tumor cell self-renewal and growth, but later leads to SASP-mediated tumor recurrence²². This study confirmed the detrimental function of persistent senescent cells in glial tumors and suggested that senescent cells could represent an actionable target to mitigate the process of gliomagenesis⁶.

Recent single-cell RNA sequencing studies classified the intra-tumoral heterogeneity of malignant GBM cells²³⁻²⁷, which can be subdivided into four main cellular plastic states: oligodendrocyte progenitor cell-like (OPC-like), neural progenitor cell-like (NPC-like), astrocyte-like (AC-like) and mesenchymal-like (MES-like) states²³. The relative abundance of these cellular states within the tumor defines three GBM transcriptomic subtypes, with proneural (PN-GBM) and classical (CL-GBM) GBMs associated with neurodevelopmental programs and mesenchymal GBM (MES-GBM) associated with injury response programs²⁴⁻²⁹. OPC-like and NPC-like states are enriched in PN-GBM whereas AC-like and MES-like states are enriched in CL-GBM and MES-GBM, respectively²³. Notably, stemness programs are heterogeneous even within a single tumor, and PN- and MES-GSCs could contribute to the genetic heterogeneity observed in patient GBM²⁴,²⁶,³⁰. Each transcriptional GBM subtype is associated with distinct molecular alterations and patient outcomes. MES-GBM is correlated with enhanced activation of anti-inflammatory (or tumor-promoting) macrophages²⁹,³¹⁻³³. Mutations in *NF1*, *TP53*, and *PTEN* genes, and increased NF-κB signaling are prevalent in this GBM subtype²⁸. Interestingly, *PTEN* loss induces cellular senescence and activates NF-κB signaling, which initiates and maintains the SASP³⁴⁻³⁶. Together these findings support the idea that cellular senescence could contribute to the intra-tumoral heterogeneity of GBM.

In this study, we investigated whether cellular senescence participates in GBM tumor progression using patient-resected GBM tissues and a mouse GBM model³⁷. We identify senescent cells in patient and mouse GBMs. Partial removal of senescent cells expressing high levels of *p16*<sup>INK4a</sup> using a ganciclovir-inducible p16-3MR transgenic line¹⁷ improves the survival of GBM-bearing mice. To identify the cells expressing high levels of *p16*<sup>INK4a</sup>, and to characterize the action of these cells on the tumor ecosystem, we combined single-cell and bulk

RNA sequencing (RNAseq) analysis at early and late time points after the senolytic treatment. This approach leads to the identification of the NRF2 transcription factor and its selected targets as a signal triggering the pro-tumoral activity of *p16*<sup>INK4a</sup> expressing senescent cells. Using these data, we define an unbiased senescence signature that we successfully used to interrogate GBM patient data sets, revealing that higher senescence scores correlate with shorter survival times.

## Results

### Identification of senescent cells in patient and mouse gliomas

We first searched for senescent cells in 28 freshly resected diffuse gliomas from patients by performing SA-β-gal staining coupled with immunohistochemistry (IHC) (14 GBMs, 5 astrocytomas, 9 oligodendrogliomas; Supplementary Fig. 1a). Ninety-eight percent of SA-β-gal positive (SA-β-gal+) cells were negative for the cell cycle marker Ki67 (Ki67−) and 72% of SA-β-gal+ cells were positive for the cell cycle inhibitor *p16*<sup>INK4A</sup>, strongly suggesting that the majority of SA-β-gal+ cells were senescent (Fig. 1a; SA-β-gal+ Ki67−: 98.35 ± 1.96%; SA-β-gal+ *p16*<sup>INK4A+</sup>: 72.51 ± 10.52%; *n* = 4 *p16*<sup>INK4A</sup>-non deleted tumors). In gliomas harboring a mutation of p53, few SA-β-gal+ cells expressed high levels of mutant p53 (Fig. 1a and Supplementary 1a) showing that some SA-β-gal+ cells were malignant cells. To identify further senescent cells, we used cell type-specific markers. Some SA-β-gal+ cells co-expressed GFAP, which could either represent parenchymal astrocytes or malignant cells, OLIG2, an OPC marker, or IBA1, a microglia/macrophage marker (Fig. 1a). To establish a quantitative measure of senescent cell burden, we quantified the percentage of the tumor area containing SA-β-gal+ cells, and used these measures to stratify tumors into three senescent categories: (1) >1% (*n* = 10/28) but below 7%; (2) ≤1–0.1% (*n* = 13/28); and (3) ≤0.1% (*n* = 5/28) senescent cells (Supplementary Fig. 1a and b). Notably, no diffuse glioma types nor the sex of the patient were associated with a particular senescent category.

We next investigated whether specific molecular alterations were associated with each of the senescent tumor categories, as defined by SA-β-Gal cell percentages. Homozygous deletion of *CDKN2A*, encoding for *p16*<sup>INK4A</sup>, is carried by 54% of patient GBMs (cbioportal.org). As *p16*<sup>INK4A</sup> is a mediator of senescence, we annotated *p16*<sup>INK4a</sup> status of each tumor, as well as examined other common molecular alterations, including p53, PTEN, NF1, and EGFR mutations. Notably, we did not find any association of a specific molecular alteration with a single senescent category (Supplementary Fig. 1a).

We then studied senescence in an immuno-competent GBM mouse model, employing a modified version of a model developed by Friedmann-Morvinski et al.³⁷. This model recapitulates the molecular alterations identified in MES-GBM: the loss of *Pten* and *p53* and the inactivation of *Nf1* triggered by the ectopic expression of H-RasV12 (Fig. 1b). Six-to-eight week-old GlastreERT2/+;Pten<sup>fl/fl</sup> female mice were intracranially injected with a lentivirus encoding H-RasV12-IRES-eGFP and shp53, into the subventricular zone (SVZ). Mice were sacrificed when they reached disease endpoints, hereafter referred to as late timepoint (Fig. 1b). These tumors displayed a heterogenous histo-pathology similar to that described in patient GBMs³⁵ (Supplementary Fig. 1c). By qPCR analysis, elevated expression levels of *Ink4/ARF* (encoding *p16*<sup>INK4a</sup>, *p15*<sup>INK4b</sup>, *p19Arf*) and *p21*, both of which encode senescence-mediating proteins, were detected in the tumor (GFP+) cells compared with the surrounding parenchyma (GFP−) cells (Fig. 1c). In contrast, *p53* mRNA levels were similarly low within GFP+ tumor cells and adjacent GFP− cells, in agreement with the presence of shp53 in the lentivirus (Fig. 1c). We observed, similarly to patient tumors, that 95% of SA-β-gal+ cells were Ki67− (94.56 ± 2.27%, *n* = 6) and 75% were *p19Arf*+ (75.17 ± 4.37%; *n* = 5). Again, these results strongly suggest that the majority of SA-β-gal+ cells were senescent. Notably, *p16*<sup>INK4a</sup> protein could not be examined in mouse tissues due to the lack of a suitable
<figure><img src="image_4.png" /><figcaption>Tumor 15</figcaption></figure>

<caption>Tumor 3</caption>

<caption>b</caption>

adult mice (6-8 weeks)

mouse genetics: GLAST<sup>CreERT2/+</sup>; Pten<sup>fl/fl</sup>

lentivirus (lv): H-RasV12-GFP-shp53

<figure><img src="image_5.png" /><figcaption>c</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>d</figcaption></figure>

<figure><img src="image_7.png" /></figure>

antibody. We also identified SA-β-gal+ LAMINB1- cells, illustrating the characteristic loss of the nuclear envelope integrity of senescent cells. Within the murine tumor, senescent cells were of distinct types, and included either malignant (GFP+) tumor cells, glial cells (GFAP+, OLIG2+), or microglia/macrophage (IBA1+) (Fig. 1d). We did not detect any senescent endothelial cells (CD31+; Fig. 1d). In general, the senescent cells were sparsely distributed in the tumor, and mostly located in

proliferative areas or adjacent to necrotic regions (Supplemen-
tary Fig. 1d).

Altogether these data reveal that cellular senescence is associated
with primary gliomagenesis, including in the mouse GBM model,
which recapitulates the histopathology, senescence features, and cell
identities of patient GBMs. We thus further used this model to address
the function of senescent cells during primary gliomagenesis.
**Fig. 1 | Identification of senescent cells in patient and mouse gliomas.**

a Representative SA-β-gal staining (blue) coupled with immunohistochemistry (IHC, brown) on two non-fixed patients GBM cryosections samples (Ki67 and GFAP: *n* = 16; *p16*<sup>INK4A</sup>; *n* = 12; IBA1: *n* = 10; p53: *n* = 7 and OLIG2: *n* = 6 patient GBMs). b Left: genetics of the mouse mesenchymal GBM model (mouse and injected lentivirus (lv)). The timeline represents the induction of tumorigenesis with tamoxifen intraperitoneal (i.p.) injections (TMX, 100 mg/kg/day for 5 days). Brains are harvested when mice reach endpoints. Right: the representative stereomicroscopic image of a mouse brain with a GFP+ GBM. c Relative transcript levels shown as ratios of normalized values of mouse GBM (GFP+, *n* = 5) over surrounding

parenchyma (GFP-, *n* = 4). Data are represented as the mean ± SEM. Statistical significance was determined by the Wilcoxon–Mann–Whitney test (*p* < 0.05; ns, not significant). d Representative SA-β-gal staining (blue) coupled with IHC (brown) on mouse GBM cryosections. (Ki67, p19, IBA1 and GFP: *n* = 8; GFAP: *n* = 6; OLIG2 and CD31: *n* = 5; LMNB1: *n* = 4 independent mouse GBMs). Arrowheads in a, d point to SA-β-gal+ cells co-labeled for the markers OLIG2, GFP, GFAP, p19<sup>AR</sup>, or IBA1. For Ki67, LMNB1, and CD31 IHC, the arrowheads point to mono-labeled SA-β-gal+ cells. Scale bars, a and d: 20 µm. H hematoxylin, HE hematoxylin, and eosin, i.p. intraperitoneal, TMX tamoxifen. Raw data are provided as a Source Data file.

## Senescent cell's partial removal increases the survival of GBM-bearing mice

We introduced the *p16-3MR* transgene in the mouse GBM model, which allowed the selective removal of senescent cells expressing *p16*<sup>INK4a</sup> with ganciclovir (GCV) injections¹⁷. Levels of the *p16-3MR* transgene were elevated in the tumor (GFP+) cells compared with the surrounding parenchyma (GFP−) cells, similar to *p16*<sup>INK4a</sup> expression, suggesting that in our model *p16-3MR* expression followed the same trend as the endogenous expression of *p16*<sup>INK4a</sup> (Supplementary Fig. 2a). Remarkably, the median survival of GBM-bearing mice harboring *p16-3MR* that were treated with GCV (p16-3MR+GCV) increased significantly compared with WT mice treated with GCV (WT+GCV) or p16-3MR mice treated with vehicle (p16-3MR+vhc) (Fig. 2a–c). Similarly, the survival of GBM-bearing mice treated with the senolytic drug ABT263 (Navitoclax, an inhibitor of the anti-apoptotic proteins BCL2 and BCL-XL<sup>13</sup>) increased significantly compared with control mice (WT+vhc) (Fig. 2b and d). Together these results strongly suggest that senescent cells act as a pro-tumoral mechanism during primary gliomagenesis.

To confirm the tumor-promoting function of senescent cells, we further studied GBM mice carrying the p16-3MR transgene. First, we analyzed whether senescence hallmarks decreased in p16-3MR+GCV GBMs compared with controls at the late time point (i.e., disease endpoint). We quantified the percentage of the tumor area (defined by GFP expression) encompassing SA-β-gal cells and found that it decreased 2.2-fold (from 2.26% to 1.02%) in p16-3MR+GCV tumors compared with WT+GCV GBMs (Fig. 2e and f). On average, about 2% of the tumor area was comprised of SA-β-gal cells in WT+GCV GBMs, which corresponds to senescent category one as we defined using patient gliomas (Supplementary Fig. 1a).

We next performed bulk RNA sequencing (RNAseq) of the tumors with or without senescent cells. In agreement with the inter-tumoral heterogeneity of patient GBMs, heat maps of the bulk RNAseq data revealed inter-tumoral heterogeneity of mouse GBMs independent of the treatment (Supplementary Fig. 2b and d). Gene set enrichment analysis (GSEA, Supplementary Fig. 2f) of p16-3MR+GCV GBMs compared with WT+GCV GBMs revealed an upregulation of cell cycle components (E2F targets), a downregulation of pathways involved in cancer (Notch signaling, mTORC1 signaling, epithelial–mesenchymal transition, angiogenesis), and modulation of the immune system (TNFA signaling via NFKB, Interferon responses, IL2-Stat5 signaling). In addition, bulk RNAseq analysis revealed a slight decrease in *p16*<sup>INK4a</sup> levels in p16-3MR+GCV GBMs compared with control GBMs (Fig. 2g). By qPCR analysis, significantly decreased expression levels of *p16*<sup>INK4a</sup>, *p15*<sup>INK4b</sup>, *p19*<sup>ARf</sup> were detected in treated GBMs compared with control GBMs (Fig. 2h). Finally, GSEA revealed a significant downregulation of senescence pathways (Fig. 2i, Supplementary Fig. 2g; Supplementary Data 1). SASP genes whose expression was significantly decreased in p16-3MR+GCV compared with WT+GCV GBMs included *Fn1*, *Plau*, *Timp1*, *Ereg*, and *Bmp2*<sup>39–41</sup>, the qPCR analysis further validated *Ereg* decrease (Fig. 2j, Supplementary Fig. 2c, e, and h). These SASP genes encode growth factors and extracellular matrix components or remodelers.

Collectively our data show that at the late timepoint, when mice were sacrificed due to tumor burden, there was an increased survival of

GBM-bearing mice associated with the partial removal of *p16*<sup>INK4a</sup> senescent cells, therefore pointing to the tumor-promoting action of senescent cells during gliomagenesis.

## Identification of p16<sup>INK4a Hi</sup> cells in a subset of malignant cells

To unveil the identity of senescent cells expressing *p16*<sup>INK4a</sup> and targeted by the p16-3MR transgene with GCV<sup>17</sup>, we performed droplet-based single cell RNAseq (scRNAseq) on FACs sorted cells from WT and p16-3MR GBM cells collected 7 days after the last GCV injection, hereafter named early timepoint (Fig. 3a and b; Supplementary Data 1). At this stage, WT+GCV GBMs (*n* = 2) exhibit increased tumor growth compared with p16-3MR+GCV GBMs (*n* = 2) (Supplementary Fig. 3a–c). Uniform manifold approximation and projection (UMAP) clustering at 0.5 resolution revealed 22 clusters with distinct gene expression signatures in each sample in the two conditions (Fig. 3c and Supplementary Fig. 3d). Non-malignant cells were identified based on the expression of the pan-leukocyte marker *CD45* (*Pptrc*), and malignant cells were identified by their expression of the 3' long terminal repeat (3'LTR) of the injected lentivirus and by copy number variations (CNV) (Fig. 3d and Supplementary Fig. 3e). Cells in each of the 22 clusters expressed variable levels of *p16*<sup>INK4a</sup> (*Cdkn2a*) however, only malignant tumor cells expressed high levels of *p16*<sup>INK4a</sup>. The p16-3MR mice were used in different cellular contexts. Injection of GCV always decreased significantly *p16*<sup>INK4a</sup> levels. However, this decrease never exceeded *p16*<sup>INK4a</sup> basal expression levels corresponding to the levels observed in the organ in absence of senescent cells<sup>17,18</sup>. We hypothesized that cells expressing *p16*<sup>INK4a</sup> at a level ≥4 (hereafter, we refer to *p16*<sup>INK4a Hi</sup> cells) were the cells targeted by p16-3MR with GCV (Fig. 3d). This threshold was chosen as *p16*<sup>INK4a Hi</sup> cells represent 3% (412/13563) of the tumor cells, a percentage that is in agreement with the area of SA-β-gal staining in the tumors (Supplementary Fig. 1a and Fig. 2f).

We then focused our analyses on the malignant cell compartment. The *p16*<sup>INK4a Hi</sup> cells were mostly present in cluster 0 which comprises the highest cell number in WT+GCV GBMs (2910 out of 13 563 cells; Fig. 3d). Further UMAP clustering of malignant cells at 0.6 resolution identified 17 clusters in the two conditions (Fig. 3e). GSEA using the mouse gene lists published by Weng et al.<sup>42</sup> allowed the malignant cell clusters to be assigned cellular identities, which predominantly included cycling cells, pri-oligodendrocyte progenitor cell-like (pri-OPC-like), committed OPC-like (COP-like), myelinating oligodendrocyte (mOL), astrocyte-like (AC), neural progenitor-like (NP-like), and hypoxic cells (HC) (Fig. 3e and f, Supplementary Fig. 3f; Supplementary Data 1). The labeling of the clusters was also in agreement with GSEA using human gene lists published by Bhaduri et al.<sup>26</sup> (Supplementary Fig. 3g). Some clusters exhibited mixed cell identities. The astrocyte cluster shared gene signatures of astrocytes, endothelial cells, and ependymal cells whereas the pri-OPC-like 1 (pOPC1) and pri-OPC-like 2 (pOPC2) clusters shared gene signatures of pri-OPC-like cells, astrocytes and COP cells (Fig. 3f). Of note, the enrichment score of each subpopulation differed very little between p16-3MR+GCV and WT+GCV GBMs, except for the pOPC1-3 clusters (Fig. 3f).

The *p16*<sup>INK4a Hi</sup> cells were mainly grouped in the astrocyte cluster and to a lesser extent in the NP-like cluster (Fig. 3g). The levels of *p16*<sup>INK4a</sup> decreased significantly in the astrocyte cluster in p16-3MR
<figure><img src="image_8.png" /><figcaption>c</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>b</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>d</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>e</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>HE</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>SA-β-gal</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>f</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>g</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>h</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>i</figcaption></figure>

<figure><img src="image_18.png" /><figcaption>j</figcaption></figure>

<figure><img src="image_19.png" /></figure>

+GCV GBMs compared with WT+GCV GBMs. No other malignant and microenvironment (CD45+) clusters showed a significant difference in *p16*<sup>ink4a</sup> levels between the two conditions (Fig. 3h; Supplementary Data 2). Therefore, this analysis identifies the astrocyte cluster as senescent. On line with a senescent phenotype, the astrocyte cluster shared an inflammatory signature (gene signatures of microglia and tumor-associated macrophages) (Supplementary Fig. 3g). Remarkably,

the percentage of cells in the astrocyte cluster decreased in p16-3MR +GCV GBMs compared with WT+GCV GBMs (astrocyte cluster from 7.75% to 3.21%; Fig. 3i), in agreement with the partial removal of *p16*<sup>ink4a Hi</sup> cells by the p16-3MR transgene in the presence of GCV.

Altogether, scRNAseq analysis identifies senescent malignant *p16*<sup>ink4a Hi</sup> cells belonging to the astrocyte cluster displaying inflammatory phenotype and targeted by the p16-3MR with GCV.
**Fig. 2** | Senescent cell's partial removal increases the survival of GBM-bearing mice. a Timeline of tumorigenesis induction (lv: H-RasV12-GFP-shp53) and removal of senescent cells with GCV, 21 days post lv injection in the p16-3MR transgenic mouse. b Timeline of tumorigenesis induction (lv-luc: H-RasV12-GFP-P2A-Luc2-shp53) and removal of senescent cells with GCV in the p16-3MR mouse or with ABT263 in WT mouse when head to body bioluminescence ratio reached 2. c Kaplan-Meier survival curves (solid lines) of WT (n = 10, median survival 38 days) and p16-3MR (n = 9, median survival 51 days) mice treated with GCV as shown in (a). Kaplan-Meier survival curves (dotted lines) of p16-3MR mice treated with vhc (n = 14, median survival 36 days) or GCV (n = 15, median survival 46 days) as shown in (b). d Kaplan-Meier survival curves of WT mice treated with vhc (n = 11, median survival 34 days) or ABT263 (n = 11, median survival 46 days) as shown in (b). e Representative HE, GFP IHC, and SA-β-Gal staining on adjacent mouse GBM cryosections. Right panels represent higher magnifications of the left panels. Scale bars, left panels: 2.5 mm, right panels: 20 µm. f Quantification of the SA-β-gal area over the tumor (GFP+) area (n = 7 biologically independent animals/group). g Relative transcript levels of p16<sup>lnk4a</sup>, shown as FPKM estimates extracted from bulk

RNAseq analysis (WT+GCV, n = 5; p16-3MR+GCV, n = 9). h Relative transcript levels are shown as ratios of normalized values of p16-3MR+GCV GBMs (n = 6) over WT+GCV GBMs (n = 4). i GSEA graphs from bulk RNAseq data in p16-3MR+GCV GBMs compared with WT+GCV GBMs. SASP gene list from Gorgoulis et al.<sup>3</sup> (Supplementary Data 1). j Relative transcript levels of genes in WT+GCV and p16-3MR+GCV GBMs extracted from bulk RNAseq data (WT+GCV, n = 5; p16-3MR+GCV, n = 9). c and d Statistical significance was determined by Mantel-Cox log-rank test. b Raw p-values from the log-rank tests are included in the figure and significance is indicated by * for p-values below the 5% level after correction by the Benjamini-Hochberg procedure. f-h, j Data are represented as the mean ± SD and statistical significance was determined by the Wilcoxon-Mann-Whitney test (*p* < 0.05; **p* < 0.01, ns, not significant). TMX tamoxifen, vhc vehicle, gav. gavage, GCV ganciclovir, i.p. intraperitoneal, lv lentivirus, lv-luc lentivirus-luciferase, GSEA gene set enrichment analysis, FDR false discovery rate, NES normalized enrichment score, r. enrichment score running enrichment score. Raw data are provided as a Source Data file.

## Partial removal of p16<sup>lnk4a Hi</sup> malignant cells impacts the remaining malignant cells

We next analyzed whether the partial removal of p16<sup>lnk4a Hi</sup> cells impacted the remaining malignant cells in our GBM model. Although the size of the tumors at this early timepoint differed between p16-3MR+GCV and WT+GCV GBMs (Supplementary Fig. 3a-c), the percentage of cycling cells (pOPC3, G2/M1, G2/M2, G1/S1, G2/S2 clusters) remained stable (WT+GCV: 28.78%; p16-3MR+GCV: 28.05%). In addition, three clusters of the oligodendroglial lineage pOPC2, COP, and mOL, increased in cell proportions upon the partial removal of p16<sup>lnk4a Hi</sup> cells (pOPC2 from 6.81% to 13.10%; COP from 1.24% to 4.60%; mOL from 1.31% to 4.30%; Fig. 3i). We further validated these results using bulk RNAseq data of WT+GCV and p16-3MR+GCV GBMs collected at the early timepoint (Supplementary Fig. 3h). GSEA using the Weng et al.<sup>42</sup> gene lists showed no difference in the expression of cycling genes (Supplementary Fig. 3i). In contrast, there was an increase of transcripts associated with COP and mOL gene signatures upon the partial removal of p16<sup>lnk4a Hi</sup> cells (Fig. 3j). The increase in cell numbers (scRNAseq) and in gene signatures (bulk RNAseq) of the oligodendroglial lineage suggests a shift of the malignant cellular states upon senolytic treatment. Indeed, GSEA revealed an increase in OPC-like and NPC-like states and their associated proneural transcriptional subtype following the partial removal of p16<sup>lnk4a Hi</sup> cells. In parallel, GSEA showed a decrease in the MES-like state and the mesenchymal transcriptional subtype (Fig. 3j). Of note, no significant change in transcripts associated with the stemness gene signature defined by Tirosh et al.<sup>43</sup> was revealed between the two conditions (Fig. 3j). Remarkably, all of these phenotypic traits perdured until the late timepoint (Fig. 3j).

Altogether, based on scRNAseq analyses, p16<sup>lnk4a Hi</sup> senescent cells are a small subset of malignant cells. Their partial removal impacts the remaining malignant cells leading to a long-lasting switch to a more oligodendroglial-like phenotype and a decrease in the expression of genes signatory of mesenchymal cell identity.

## Modulation of the immune compartment following the partial removal of p16<sup>lnk4a Hi</sup> cells

The mesenchymal transcriptional GBM subtype is associated with enhanced expression of anti-inflammatory and tumor-promoting macrophages<sup>29,31,44</sup>. We, therefore, examined the immune compartment in the scRNAseq data at the early timepoint following the partial removal of p16<sup>lnk4a Hi</sup> cells. UMAP clustering of CD45+ cells revealed seven clusters in the two experimental conditions (Fig. 3d, Supplementary Fig. 3d and Fig. 4a). Differentially expressed (DE) genes and GSEA allowed the labeling of these clusters into infiltrating bone marrow-derived macrophages (BMDM), resident microglia and T cells<sup>45</sup> (Fig. 4b and c, Supplementary Fig. 4a; Supplementary Data 1). All the BMDM and microglia clusters harbored an anti-inflammatory gene signature. Furthermore, the BMDM-like1 and microglia clusters

shared an antagonist pro-inflammatory gene signature<sup>32</sup> (Fig. 4c; Supplementary Data 1). In addition, the proportion of the immune fraction within the tumor hardly varied between WT+GCV and p16-3MR+GCV GBMs. However, the number of T cells increased (from 9% to 27%), whereas the number of BMDM decreased (from 41% to 30%) upon the partial removal of p16<sup>lnk4a Hi</sup> cells (Fig. 4d). This latter phenotype was confirmed by GSEA on bulk RNAseq data, which showed a significant decrease in transcripts associated with a core BMDM signature at the early and late timepoints in p16-3MR+GCV GBMs compared with controls (Supplementary Fig. 4b). In addition, an estimation of the abundances of the main immune cell types from our bulk RNAseq data using CIBERSORT pointed to a significant decrease of BMDM upon partial removal of p16<sup>lnk4a Hi</sup> cells at the late timepoint (Supplementary Fig. 4c). Further, qPCR analysis on bulk RNA at the late timepoint showed significantly decreased expression levels of *Gda* and *Crip1*, two BMDM markers<sup>32</sup> in treated GBMs compared with controls whereas expression levels of *Tmem119* and *P2ry12*<sup>32</sup> two microglia markers, did not vary between both conditions (Supplementary Fig. 4d).

We next examined whether the activity of immune cell types was altered in GBM tumors partially depleted of senescent cells. GSEA on scRNAseq data at the early timepoint revealed an upregulation of TNFA signaling via the NFKB pathway in the microglia cluster and a downregulation of genes associated with the epithelial-to-mesenchymal transition, inflammatory and hypoxia pathways in the BMDM clusters in p16-3MR GBMs compared with WT+GCV GBMs (Fig. 4e). Close examination of the DE genes in these pathways revealed a significant increase in the expression of genes associated with a pro-inflammatory signature (*Ccl4*, *Tnf*, *Il1a*, *Il1b*, *Csf1*) in the microglia cluster and a significant decrease in the expression of genes related to an anti-inflammatory signature (*Cxcl2*, *Vegfa*, *Tgfbi*, *Spp1*, *Thbs1*, *Hmox1*, *Hif1a*) in the BMDM clusters (Fig. 4f, Supplementary Data 2). In addition, T cell cluster analysis revealed a decrease in the expression of genes regulating the activity of T cells, including the immune checkpoint genes, *Ctla4*, *Lag3*, and *Pdcd1* (encoding PD1) (Fig. 4f). Consistent with these data, GSEA of bulk RNAseq data revealed a decrease in transcripts associated with an anti-inflammatory pathway following the partial removal of p16<sup>lnk4a Hi</sup> cells at the early and late timepoints (Supplementary Fig. 4b).

Collectively, the transcriptomic analysis at single-cell and bulk levels shows that the partial removal of p16<sup>lnk4a Hi</sup> malignant cells modulates the abundance and the activity of tumor-associated macrophages.

## Identification of NRF2 activity and its putative targets in p16<sup>lnk4a Hi</sup> malignant cells

To explore the regulators of senescence in p16<sup>lnk4a Hi</sup> malignant cells, we performed pathway enrichment analysis with the ENCODE and ChEA
<figure><img src="image_20.png" /><figcaption>b</figcaption></figure>

<figure><img src="image_21.png" /><figcaption>c</figcaption></figure>

<figure><img src="image_22.png" /><figcaption>d</figcaption></figure>

<figure><img src="image_23.png" /><figcaption>e</figcaption></figure>

<figure><img src="image_24.png" /><figcaption>g</figcaption></figure>

<figure><img src="image_25.png" /><figcaption>f</figcaption></figure>

<figure><img src="image_26.png" /><figcaption>h</figcaption></figure>

<figure><img src="image_27.png" /><figcaption>i</figcaption></figure>

<figure><img src="image_28.png" /><figcaption>j</figcaption></figure>

<figure><img src="image_29.png" /></figure>

consensus TFs from ChIP-X database using Enrichr¹⁶ on three gene sets enriched in p16<sup>ink4a</sup> Hi senescent cells: (i) differentially downregulated genes in the p16-3MR+GCV vs. the WT+GCV astrocyte cluster from the scRNAseq data at the early timepoint (Early; Fig. 5a and b; Supplementary Data 3); (ii) differentially upregulated genes in p16<sup>ink4a</sup> positive vs. p16<sup>ink4a</sup> negative malignant cells from scRNAseq analysis at the late timepoint (Late (1); Fig. 5c and d, Supplementary Fig. 5a-e;

Supplementary Data 3); (iii) differentially downregulated genes in the p16-3MR+GCV vs. the WT+GCV GBMs from the bulk RNAseq data at the late timepoint (Late (2); Fig. 5e and f; Supplementary Data 3). Remarkably, the *Nuclear Factor Erythroid 2 Like 2 (Nfe2l2)* signaling pathway was enriched in the three gene sets. NRF2 encoded by *Nfe2l2* is an antioxidant defense system that appears to be a plausible candidate to trigger the pro-tumoral action of p16<sup>ink4aHi</sup> senescent cells as it
**Fig. 3 | Identification of p16<sup>ink4a hi</sup> cells in a subset of malignant cells.** a Timeline of the mouse GBM generation for scRNAseq at the early timepoint. b Scheme of the scRNAseq experiment. c UMAP plots of WT+GCV (*n* = 2) and p16-3MR+GCV (*n* = 2) GBM cells at a 0.5 resolution and annotated malignant cells and TME cells. d Violin plots of the expression of CD45, 3'LTR, and p16<sup>ink4a</sup> in WT+GCV GBM cells per cluster. e UMAP plots of WT+GCV (*n* = 2) and p16-3MR+GCV (*n* = 2) GBM malignant cells and annotated cell type at a 0.6 resolution. f GSEA dot plots of DE genes (FDR < 0.05; avlogFC>0.25) in WT+GCV (gray dots) and p16-3MR+GCV (red dots) GBMs of gene lists from Weng et al.<sup>42</sup> (Supplementary Data 1). g Violin plots of the expression of p16<sup>ink4a</sup> in malignant cells per cluster. The red box indicates the cells with an expression of p16<sup>ink4a</sup> ≥ 4 (p16<sup>ink4a hi</sup> cells). h Bar plots representing the significance of p16<sup>ink4a</sup> fold change per cluster in p16-3MR+GCV GBMs compared

with WT+GCV GBMs. The arrowheads point to a decrease (arrowheads down) or increase (arrowheads up) in the fold change. i Bar plots representing the percentage of malignant cells per cluster in WT+GCV and p16-3MR+GCV GBMs. The arrowheads point to clusters whose cell number varies between the two conditions. j GSEA ridge plot of gene lists from Weng et al.<sup>103</sup>, Neftel et al.<sup>23</sup>, and Wang et al.<sup>29</sup> (Supplementary Data 1) between p16-3MR+GCV and WT+GCV GBMs at the early and late time points. Analysis performed from bulk RNAseq data. TMX tamoxifen; GCV ganciclovir; vhc vehicle; i.p. intraperitoneal; Iv-lv lentivirus-luciferase; TME tumor microenvironment, UMAP uniform manifold approximation and projection, LTR long terminal repeat, DE differentially expressed, GSEA gene set enrichment analysis, FDR false discovery rate, NES normalized enrichment score, r. enrichment score running enrichment score. Raw data are provided as a Source Data file.

induces cellular senescence in fibroblasts<sup>47</sup> and confers a selective advantage in cancer cells<sup>48</sup>. Among the identified NRF2 putative targets, three genes were common to all three data sets (*Tgif1*, *Plaur*, *Gja1*) whereas eight genes were shared between two of the three gene lists (*Dap*, *Esd*, *Lmna*, *Areg*, *Igfbp3*, *Cdkn2b*, *Tnc*, and *Peak1*) (Fig. 5g; Supplementary Data 3). As illustrated on the heatmap, the combined expression of *Nfe2l2* and 11 putative target genes were unique to p16<sup>ink4a hi</sup> cells in WT+GCV GBMs (Fig. 5h).

Immunohistochemistry on GBM cryosections collected at the late timepoint revealed that NRF2 was expressed in a few scattered cells (Fig. 5i and Supplementary Fig. 5f). Quantification of the NRF2 expression area in the tumor showed a modest decrease in p16-3MR+GCV compared with WT+GCV GBMs (Fig. 5j). Of note, the expression of NRF2 in cells expressing low levels of p16<sup>ink4a</sup>, most probably CD45+ cells, may have concealed decreased NRF2 expression in senescent cells (Fig. 5h). We further examined the expression of three NRF2 putative target genes whose encoded proteins are associated with senescence, glioma progression or glioma resistance, respectively, namely urokinase plasminogen activator receptor (uPAR) encoded by *Plaur*<sup>49</sup>, Tenascin-C (TNC)<sup>50</sup> or Connexin43 (CX43) encoded by *Gja1*<sup>51</sup>. These proteins were expressed in a few scattered cells throughout the tumor. TNC was expressed in more cells than uPAR and CX43 in line with their transcript expression at the single cell level (Fig. 5h and i, Supplementary Fig. 5f). Quantification of CX43 by IHC and TNC by western blot revealed a significant downregulation of these proteins in p16-3MR+GCV GBMs compared with WT+GCV GBMs, strengthening *Gja1* and *Tnc* as NRF2 target genes in GBM (Fig. 5k and l, Supplementary Fig. 5g). We then assessed whether interactions between NRF2 selected targets and the immune fractions were observed in GBMs. We interrogated for ligand-receptor interactions between cluster 0, enriched in p16<sup>ink4a hi</sup> cells, and the immune clusters in the scRNAseq data at the early timepoint using CellPhoneDB (Fig. 3d, Supplementary Fig. 5h). In silico analysis highlighted possible interactions between TNC and PLAUR, expressed in malignant cells, and integrins receptors expressed in the immune clusters. Remarkably, putative TNC-aVb3 and PLAUR-aVb3 ligand-receptor interactions between malignant cells and T cells were abolished upon partial removal of p16<sup>ink4a hi</sup> cells (Supplementary Fig. 5h).

Altogether these data identify NRF2 activity and its putative target genes in p16<sup>ink4a hi</sup> senescent cells and suggest that this signal could in part trigger the detrimental action of senescent cells during gliomagenesis.

## Knockdown of NRF2 in malignant cells recapitulates most features of the senolytic treatment

NRF2 has pleiotropic actions depending on cellular context. Tumor-suppressing effects of NRF2 are mediated via the maintenance of a functional immune system<sup>48</sup>. For instance, in a mouse lung cancer model, NRF2 activity in immune cells contributes to suppressing tumor progression<sup>52</sup>. To directly test whether NRF2 triggers the tumor-promoting action of malignant senescent cells, we used a knockdown approach, introducing a

microRNA targeting NRF2 (miR-NRF2) into the lentivirus used to induce gliomagenesis. We analyzed the resultant tumors at the late timepoint (Fig. 6a and b, Supplementary, Fig. 6a). Quantification of NRF2 by IHC revealed a significant decrease of the protein in miR-NRF2-GBMs compared with miR-control (ctl)-GBMs (Fig. 6c and d, Supplementary Fig. 6d). Notably, NRF2 is also expressed in CD45+ cells, not targeted by our approach, which persisted in miR-NRF2-GBMs. We performed bulk RNAseq and GSEA of miR-NRF2- and miR-ctl- GBMs at the late timepoint, and found a significant downregulation of canonical NRF2 targets and NRF2 targets from the combined analysis, confirming knockdown of NRF2 using our miR-NRF2 (Fig. 6e; Supplementary Data 3).

We next asked whether the knockdown of NRF2 in malignant cells impacted cellular senescence. The percentage of the tumor area encompassing SA-β-gal cells was similar in miR-NRF2 GBMs compared with miR-ctl GBMs, suggesting that NRF2 knockdown in malignant cells does not induce the death of senescent cells (Fig. 6f and g, Supplementary Fig. 6e). However, GSEA performed on bulk RNAseq data from miRNRF2- and miR-ctl-GBMs revealed a significant downregulation of SASP genes associated with senescence (Fig. 6e). Among these genes, *Mmp1a*, *Mmp3*, *Mmp10*, *Plau*, *Col1a2*, *Timp1*, and *Thbs1* were differentially expressed between the two conditions (Supplementary Fig. 6b). This result strongly suggests that NRF2 regulates directly or indirectly the expression of SASP genes. Further, we examined whether NRF2 knockdown mimicked the phenotype of senolytic treatment. GSEA revealed a significant decrease in the expression of genes associated with a mesenchymal identity, a BMDM signature, and anti-inflammatory pathways, similar to the gene signature changes observed in p16-3MR+GCV GBMs (Fig. 6e; Supplementary Fig. 6c). However, genes associated with an oligodendroglial identity were not modulated upon NRF2 knockdown (Fig. 6e), in contrast to the effect of the partial removal of p16<sup>ink4a hi</sup> senescent cells.

Finally, as NRF2 activity protects against DNA-damaging agents and prevents carcinogenesis<sup>53</sup>, we explored whether NRF2 knockdown impacted the onset of tumorigenesis. Live bioluminescence imaging showed no difference in the onset of tumorigenesis between the two tumor types (Fig. 6h). Most importantly, the presence of miR-NRF2 in malignant cells significantly increased the survival of GBM-bearing mice compared with controls (Fig. 6i), an effect that was more marked than upon the partial removal of p16<sup>ink4a hi</sup> senescent cells (Figs. 2c, d, and 6i). One major difference between the two paradigms was the decrease of transcripts linked to the cell cycle, which occurred only in the knockdown of NRF2 in malignant cells (Fig. 6e and Supplementary 6c). One reason for this difference could be that all malignant cells were targeted by the miR strategy, whereas the p16-3MR+GCV paradigm only partially removed senescent cells.

Collectively our results show that NRF2 knockdown recapitulates most features of senolytic treatment and strongly supports NRF2 as a cellular senescence regulator in malignant cells.
<figure><img src="image_30.png" /></figure>

<figure><img src="image_31.png" /></figure>

<figure><img src="image_32.png" /><figcaption>**Fig. 4 | Modulation of the immune compartment following p16<sup>INK4a Hi</sup> cells partial removal.** a UMAP plots of CD45+ cells in WT+GCV and p16-3MR+GCV GBMs at a 0.5 resolution and annotated cell type. b Violin plots representing the expression of selected DE genes (FDR < 0.05; avlogFC > 0.25) per cluster in WT+GCV GBMs. c GSEA dot plot of DE genes (FDR < 0.05; avlogFC>0.25) in WT+GCV (gray dots) and p16-3MR+GCV (red dots) CD45+ clusters of core-BMDM, core-MG, pro-inflammatory and anti-inflammatory pathways as defined in Bowman et al.⁴⁵ and Darmanis et al.³² (Supplementary Data 1). d Chart pies representing the percentage of CD45+ cells per cluster in WT+GCV and p16-3MR+GCV GBMs. e GSEA graphs representing the enrichment score of Hallmark gene lists in p16-3MR+GCV compared with WT+GCV microglia clusters and pooled BMDM clusters. The</figcaption></figure>

<figure><img src="image_33.png" /></figure>

<figure><img src="image_34.png" /></figure>

<figure><img src="image_35.png" /></figure>

barcode plot indicates the position of the genes in each gene set; red represents positive Pearson's correlation with p16-3MR+GCV expression and blue with WT+GCV expression. f Dot plots of the relative expression of selected genes in WT+GCV and p16-3MR+GCV microglia, pooled BMDM and T cells clusters. Statistical significance of the expression of genes in p16-3MR+GCV compared with WT+GCV clusters was determined by the Wilcoxon-Mann-Whitney test (ns, not significant, \*p < 0.05; \*\*p < 0.01; \*\*\*p < 0.001). UMAP uniform manifold approximation and projection, BMDM bone marrow-derived macrophages, MG microglia, DE differentially expressed, EMT epithelial to mesenchymal transition, GSEA gene set enrichment analysis, FDR false discovery rate, NES normalized enrichment score. Raw data are provided as a Source Data file.

## Mouse senescent signature is conserved in patient GBMs and its enrichment score is predictive of a worse survival

We next explored whether p16<sup>INK4a Hi</sup> senescent cells are conserved in patient GBMs. We first established a senescent signature from scRNAseq data at the early timepoint (Fig. 3). We compared the transcriptome of p16<sup>INK4a Hi</sup> cells in astrocyte and NP-like clusters

(p16<sup>INK4a Hi</sup> group) with the remaining malignant cells in WT+GCV GBMs (Figs. 3g and 7a and b). GSEA in the p16<sup>INK4a Hi</sup> group revealed down-regulation of cell cycle pathways and an oligodendroglial state and increased expression of genes associated with inflammation, NRF2 signaling, MES-like state, mesenchymal transcriptional GBM subtype (Supplementary Fig. 7a and b). We further selected a list of 31
<figure><img src="image_36.png" /></figure>

<figure><img src="image_37.png" /></figure>

<figure><img src="image_38.png" /></figure>

<figure><img src="image_39.png" /></figure>

<figure><img src="image_40.png" /></figure>

<figure><img src="image_41.png" /></figure>

<figure><img src="image_42.png" /></figure>

<figure><img src="image_43.png" /></figure>

<figure><img src="image_44.png" /></figure>

<figure><img src="image_45.png" /></figure>

<figure><img src="image_46.png" /></figure>

<figure><img src="image_47.png" /></figure>

<figure><img src="image_48.png" /></figure>

<figure><img src="image_49.png" /></figure>

<figure><img src="image_50.png" /></figure>

genes to define a GBM senescence signature. Among the 278 differ-
entially upregulated genes (FDR < 0.05) in the p16<sup>ink4a</sup> Hi group, we
selected genes that were expressed in more than 90% of p16<sup>ink4a</sup>Hi cells
and presented a log2-fold change superior to 0.8 between the two
groups. As expected, senescence-associated genes were enriched in
the astrocyte cluster and to a lesser extent in the NP-like cluster
(Fig. 7c). The encoded proteins were associated with diverse cellular

processes compatible with cellular senescence, such as cell cycle arrest
(Cdkn1a, Cdkn2a, Cdkn2b), lysosomal function (Ctsb, Ctsd, Ctsl, Ctsz,
Lamp1, Lamp2), cellular growth (Igfbp2, Igfbp3), extracellular matrix
interaction (Sparc, Tnc, Sdc4, Lgals1, Timp1, Mt1), cytoskeleton inter-
action (Pdlim4, S100a11, Tmsb4x, Sep11) and cancer (Tm4sf1, Ociad2,
Emp3) (Fig. 7c). We then computed a single-sample GSEA (ssGSEA)
senescent Z-score corresponding to the enrichment Z-score of the 31
**Fig. 5 | Identification of NRF2 activity and its putative targets in p16<sup>ink4a</sup> Hi malignant cells.** a Timeline of the mouse GBM generation for scRNAseq at the early timepoint (EARLY). b Barplot corresponding to significantly enriched pathways (ENCODE and ChEA consensus TFs from ChIP-X, Enrichr) in differentially downregulated genes (FDR < 0.05; avlogFC > 0.25) in the p16-3MR+GCV compared with the WT+GCV astrocyte clusters from the scRNAseq data (as shown in a). c Timeline of the mouse GBM generation for scRNAseq at the late timepoint (LATE (1)). d Barplot corresponding to significantly enriched pathways in differentially upregulated genes (FDR < 0.05; logFC > 0.5) in p16<sup>ink4a</sup> positive vs. p16<sup>ink4a</sup> negative malignant cells from the scRNAseq data (as shown in c). e Timeline of the mouse GBM generation for bulk RNAseq at the late timepoint (LATE (2)). f Barplot corresponding to significantly enriched pathways in differentially down-regulated genes (FDR < 0.05; logFC > 0.5) in p16-3MR+GCV compared with WT+GCV GBMs from the bulk RNAseq data (as shown in e). g Venn diagram of NRF2 putative targets between the 3 gene sets as shown in (a, c, and e). h Heatmaps of Nrf2 and its 11 identified putative targets in WT+GCV and p16-3MR GBMs. Cells are classified into five

categories according to p16<sup>ink4a</sup> expression levels. i Representative immunohistochemistry (IHC, brown) counterstained with hematoxylin (H, purple) on mouse GBM cryosections at the late timepoint. H hematoxylin. (NRF2: WT+GCV, *n* = 7; p16-3MR+GCV *n* = 7; uPAR: WT+GCV, *n* = 3; p16-3MR+GCV, *n* = 3; CX43: WT+GCV, *n* = 5; p16-3MR+GCV, *n* = 6; TNC: WT+GCV, *n* = 5; p16-3MR+GCV, *n* = 7 independent mouse GBMs). Scale bar: 20 µm. j Quantification of the NRF2 area (IHC) over the tumor area (WT+GCV, *n* = 7; p16-3MR+GCV, *n* = 7 independent mouse GBMs). k Quantification of the CX43 area (IHC) over the tumor area (WT+GCV, *n* = 6; p16-3MR+GCV, *n* = 6 independent mouse GBMs). l Quantification of the ratio of TNC over β-TUBULIN expression (western blot) (WT+GCV, *n* = 5; p16-3MR+GCV, *n* = 7 independent mouse GBMs). Raw data are shown in Supplementary Fig. 5g, j–l data are presented as the mean ± SD. Statistical significance was determined by the Wilcoxon–Mann–Whitney test (*p* < 0.05). i.p., intraperitoneal, Iv lentivirus, Iv-luc lentivirus-luciferase, TMX tamoxifen, DE differentially expressed GCV ganciclovir, H hematoxylin. Raw data are provided as a Source Data file.

genes of the senescence signature in all malignant cells of WT+GCV and p16-3MR+GCV GBM transcriptomes. As expected, the astrocyte cluster in WT+GCV GBMs contained the largest high senescent Z-score distribution rate (Supplementary Fig. 7c). For unbiased analysis, we defined the high distribution rate as the highest decile. This percentage dropped in all clusters in p16-3MR+GCV GBMs compared with WT+GCV GBMs, as predicted by the partial removal of p16<sup>ink4a Hi</sup> senescent cells (Supplementary Fig. 7c).

We then applied the ssGSEA senescence Z-score to three single-cell data sets of patient GBMs<sup>23,26,54</sup>. We analyzed separately the Neftel dataset according to the sequencing technology (Smartseq2 (SS2) and 10X). The range of the ssGSEA senescence Z-score of cells from patient GBMs was similar to those identified in mouse GBMs (Fig. 7d; Supplementary Fig. 7c). Hence, similar senescent transcriptomic profiles were observed in cells from mouse and patient GBMs. GBMs from the three data sets contained high senescent Z-score distribution rates, with the exception of 2/31 tumors, possibly due to the small number of cells sequenced in these samples (MGH126: 201 cells and MGH151: 151 cells). In summary, this analysis strongly suggests that the mouse senescent signature is conserved in cells from patient GBMs.

To assess whether the ssGSEA senescence score could be used as a prognostic factor for patients with GBM, we interrogated The Cancer Genome Atlas (TCGA) GBM data sets and performed a Cox regression analysis with five variables. Cellular senescence was linked to aging and our mouse senescence signature was based on the expression of p16<sup>ink4a</sup>. Therefore, we used as variables the ssGSEA senescence score, p16<sup>ink4a</sup> copy number alteration status, and the age of the patient, in addition to the Karnofsky score and sex variables commonly used for patients with GBM (Supplementary Data 4). Multivariate Cox regression model showed that regardless of p16<sup>ink4a</sup> status, the age of the patient, the sex, and the Karnofsky score, the enrichment of the senescence score predicted a worse survival (hazard ratio above 1) in patients with GBM (Fig. 7e; Supplementary Data 4).

Finally, we tested whether NRF2 activity could account for some of the tumor-promoting action of cellular senescence in patient GBMs, similar to mouse GBMs. First, SA-β-gal staining coupled with IHC on cryosections revealed the expression of NRF2, TNC, and CX43 in SA-β-gal+ cells in patient GBM samples (Fig. 7f and Supplementary 7d). As described above, we next interrogated ssGSEA NRF2 target scores in TCGA GBM data sets and performed a Cox regression analysis. NRF2 targets corresponded to the 59 genes identified in the combined analysis (Supplementary Data 3; Fig. 5g). The Cox regression model showed that regardless of p16<sup>ink4a</sup> status, the age of the patient, the sex, and the Karnosky score, an enriched NRF2 target gene score predicted worse survival in patients with GBM (Fig. 7g; Supplementary Data 4).

In summary, our data show that cells enriched for the mouse senescent signature are present in patient GBMs and that the

enrichment scores of senescence and of NRF2 targets are correlated with worse survival in patients with GBM.

## Discussion

Depending on the context, cellular senescence plays both beneficial and detrimental roles during tumor progression. Here, we revealed the tumor-promoting action of malignant senescent cells in mouse and patient GBMs. The mouse MES-GBM model used in the present study, even though its genetics differ from the patient GBMs, recapitulated the histopathology, the heterogeneity of cellular states, the infiltration of BMDM specific to the mesenchymal transcriptional GBM subtype, and the senescent features of patient GBMs (see also ref. 44). Partial removal of p16<sup>ink4a Hi</sup> malignant senescent cells modified the tumor ecosystem and improved the survival of GBM-bearing mice. The difference in survival following a senolytic treatment appeared to be relatively modest, nonetheless, this difference was significant and was observed in the p16-3MR paradigm when compared to a first control cohort (WT+GCV). These results were repeated using a second control cohort (p16-3MR+vhc). In addition, we observed a benefit effect of another senolytic, ABT263, in a treated cohort compared to a control one. This finding is remarkable given the fact that senescent cells represented <7% of the tumor and that their removal using the p16-3MR transgene was only partial. These findings suggest that senolytic drug therapy may be a benefit for patients with GBM.

By combining single-cell and bulk RNA sequencing, immunohistochemistry, and genetic knockdowns, our study established a link between senescence and NRF2 activity in the context of GBM. Chronic activation of NRF2 contributes to tumor growth, metastasis, treatment resistance, and poorer prognosis in patients with cancer<sup>48</sup>. Depending on the context, NRF2 promotes or delays fibroblast senescence<sup>47,48</sup>. NRF2 binds to antioxidant-responsive elements (AREs) and controls the expression of a battery of genes regulating metabolism, intracellular redox balance, apoptosis, and autophagy<sup>48</sup>. Under a homeostatic state, cytoplasmic NRF2 binds to KEAP1, which mediates its proteasomal degradation. Impairment of NRF2-KEAP1 binding, either by phosphorylated-p62/SQSTM1 or by elevated ROS permits NRF2 nuclear translocation and subsequent activation of target genes<sup>55–57</sup>. Of note, p62/SQSTM1-mediated degradation of KEAP1 and NRF2 promotes in vitro glioma stem cell survival<sup>58</sup>. NRF2 is hyperactivated preferentially in the MES-GBM subtype<sup>58</sup> however, NRF2/NFE2L2 and genes regulating its activity (KEAP1, SQSTM1, CUL3, RBX1, SKP1, CUL1, BTRC, SYVNI) are rarely altered in GBMs (see cbioportal.org). The putative triggers of senescence in our GBM model and in patient GBMs are also known to regulate NRF2 activity such as hypoxia, ROS, PI3K–AKT pathway (enhanced by the loss of PTEN)<sup>48,59</sup> or Nrf2 transcription such as RAS oncogene (K-RAS)<sup>60</sup>. Future work is required to establish the contribution of these triggers in the process of cellular senescence and NRF2 activity in GBMs. In the present study, we identified NRF2
<figure><img src="image_51.png" /><figcaption>c</figcaption></figure>

<figure><img src="image_52.png" /><figcaption>e</figcaption></figure>

<figure><img src="image_53.png" /><figcaption>f</figcaption></figure>

<figure><img src="image_54.png" /><figcaption>h</figcaption></figure>

<figure><img src="image_55.png" /><figcaption>b</figcaption></figure>

<figure><img src="image_56.png" /><figcaption>d</figcaption></figure>

<figure><img src="image_57.png" /><figcaption>g</figcaption></figure>

<figure><img src="image_58.png" /><figcaption>i</figcaption></figure>

<figure><img src="image_59.png" /></figure>

putative targets that are not canonical NRF2 targets. These targets encode for growth factors (AREG, IGFBP3, TGIF1), ECM components or remodelers (TNC, uPAR, ESD), or cell-cell interactors (CX43) and have been previously identified as SASP factors³⁹⁻⁴¹. TNC and CX43 are of particular interest. Indeed, CX43 participates in the formation of microtubes that interconnect malignant cells, creating a cellular network resistant to treatment⁵¹. Further, the pro-tumoral functions of

TNC have been described, independently of the senescence context⁶¹,⁶². TNC is a component of the glioma ECM that binds to integrin receptors, EGF receptor (EGFR), and SYNDECAN 4 (SDC4) and regulates angiogenesis, proliferation, and cell migration⁵³. Hence, TNC functions could partly be responsible for the tumor-promoting phenotype of the p16 Ink4a H1 malignant senescent cells. As mentioned above, NRF2 has pleiotropic actions depending on the cellular context and is
**Fig. 6** | Knockdown of NRF2 in malignant cells recapitulates most features of the senolytic treatment. a Scheme of the lentiviral vector containing either a miR-NRF2 or a miR-ctl. b Timeline of the mouse GBM generation at the late timepoint. c Representative NRF2 IHC staining (brown) on miR-ctl (n = 4) and miR-NRF2 (n = 4) GBM cryosections. Necrotic areas are outlined in red dashed lines. d Quantification of the NRF2 area (IHC) over the tumor area (miR-ctl n = 4; miR-NRF2 n = 4). e GSEA ridge plot on bulk RNAseq of miR-NRF2-GBMs compared with miR-ctl-GBMs (see Supplementary Data 1 for gene lists). f Representative SA-β-gal (blue) staining on miR-ctl (n = 4) and miR-NRF2 (n = 4) GBM cryosections. Necrotic areas are outlined in red dashed lines. g Quantification of the SA-β-gal area over the tumor area (miR-

ctl n = 4; miR-NRF2 n = 4). h Boxplot representing the onset of tumorigenesis in miR-ctl (n = 10) and miR-NRF2 (n = 10) mice following post-lentiviral injection. The onset of tumorigenesis is defined when the bioluminescence reached 3e10⁶. i Kaplan-Meier survival curves of miR-ctl (n = 10, median survival 38.5 days) and miR-NRF2 mice (n = 10, median survival 55.5 days). Statistical significance was determined by the Mantel-Cox log-rank test (**p < 0.01). Scale bar, c and f 50 µm. d, g, h Statistical significance was determined by Wilcoxon-Mann-Whitney test (*p < 0.05). Iv lentivirus, miR-ctl miR-control, H hematoxylin, GSEA gene set enrichment analysis, sen. senescence, TAM-associated macrophages, BMDM bone marrow-derived macrophages. Raw data are provided as a Source Data file.

expressed in multiple cell types. Of note, in this study, we did not address the function of NRF2 in the GBM microenvironment, notably in CD45+ cells. To sum up, our findings suggest that a senolytic treatment may represent a therapeutic strategy to eliminate NRF2+ malignant senescent cells without targeting NRF2+ cells in the microenvironment.

Single-cell RNAseq analysis of mouse GBMs allowed the comprehensive characterization of the pro-tumorigenic malignant senescent cells. Although our approach focused primarily on p16<sup>INK4a</sup>HI senescent cells in a mouse MES-GBM model, our findings show that the senescence signature we established in this study, is applicable to GBMs regardless of p16<sup>INK4a</sup> status (Fig. 7e). The presence in the senescence signature of three genes (*Cdkn2a*, *Cdkn2b*, *Cdkn1a*) encoding for cyclin-dependent kinase inhibitors warrants cell cycle exit and entry into senescence. Of note, *CDKN1A* (*p21<sup>CIP1</sup>*) is rarely mutated in patient GBMs (0.4%) and *p21<sup>CIP1</sup>* mediates senescence in many tissues³. Furthermore, we presume that the senescent signature defined in the study is specific to detrimental senescence. Indeed, the enrichment of the senescence score predicts a worse survival in patients with GBM and multiple genes in the signature encode for proteins whose activities are associated with tumor aggressiveness and/or worse patient prognosis (CD151<sup>63</sup>; EMP3<sup>64</sup>; IGFBP2<sup>65</sup>; LGALS1<sup>66</sup>; TMSB4X<sup>67</sup>; TNC/SDC4<sup>62</sup>; SPARC<sup>68</sup>; TIMP1<sup>69</sup>). Future studies will determine whether the senescence scoring could be used in the diagnosis of patients with GBM to improve the design of personalized treatment and effective combinatorial strategies.

In this study, we showed that senolytic treatments applied to GBM-bearing mice delay temporarily tumor growth (Fig. 2c and d). Previous work investigated the action of the specific inhibitors of the anti-apoptotic BCL2 and BCL-XL proteins, ABT737 and ABT263, in the context of GBM mouse models. Similarly to our results (Fig. 2d), ABT737 treatment increases the survival of immunodeficient mice grafted with the human glioma cell line U-251MG<sup>70</sup>. ABT263 treatment when combined with drugs decreasing McI1 levels (a BCL2 family protein), decreased tumor volume in a heterotopic (subcutaneous) model of proneural GBM<sup>71,72</sup>. Furthermore, ABT263 treatment provides a survival benefit only when applied in combination with the oncometabolite 2-R-2-hydroxyglutarate (2-HG; produced by IDH1-mutated tumors) in a proneural GBM model<sup>73</sup>. Cellular senescence was not addressed in these studies but in the light of our results using the p16-3MR transgene, some of the cells targeted by this approach may be malignant senescent cells, extending the notion of detrimental senescent cells to distinct GBM subtypes in agreement with the senescence score analysis performed on data from patient GBMs (Fig. 7). All together our findings and these previous studies, raise the possibility to use senotherapy to improve the outcome of a patient with GBM.

Many important issues remain to be solved before envisioning senotherapy in the context of GBM such as the nature of the senolytic, the timing of its administration (neoadjuvant-concomitant adjuvant or adjuvant), and the most effective combination of treatment. First of all, these drugs should cross the brain-blood-barrier. Novel molecules are expected to be discovered in the near future as the field of senotherapies, which includes drugs eliminating senescent cells (senolytic drugs: anti-BCL2 and BCL-XL<sup>13</sup>, dasatinib and quercitin<sup>14</sup>, cardiac glycosides<sup>74-76</sup>, CAR T cells<sup>77</sup>) and drugs inhibiting their function

(senostatic drugs such as Metformin<sup>78</sup>), is under active investigation<sup>79</sup>. Despite encouraging results in mouse GBM models, clinical studies show major side effects of ABT263 such as thrombocytopenia and neutropenia caused by BCL-XL inhibition<sup>80</sup>, limiting the use of BCL-XL inhibitors as safe and effective anticancer agents. The selective BCL2 inhibitor ABT199 (Venetoclax) which spares the platelets, displays variable results on senescent glioma cells, in vitro<sup>81,82</sup>. Administration of a combination of dasatinib and quercitin, a multi-tyrosine kinase inhibitor and a natural flavonoid respectively, eliminates efficiently senescent cells in different pathologies in the mouse including neurodegenerative diseases<sup>83,84</sup>. Results of ongoing clinical trials using these drugs on cohorts of patients with mild cognitive impairment and Alzheimer's disease (NCT04785300i, NCT04685590ii) should encourage or refute this approach for targeting senescence in the central nervous system.

Heterotopic and orthotopic mouse models of cancers showed the efficacy of the one-two punch sequential therapy defined by a pro-senescence therapy, followed by senolytic therapy clearing the induced-senescent cells and preventing the accumulation of detrimental persistent senescent cells<sup>18,77,85-87</sup>. The efficacy of such a strategy implies a homogeneous response to the pro-senescence therapy. GBMs which are characterized by intra-tumoral heterogeneity, are thus not the ideal candidate for the use of this strategy. Nonetheless, radiation and TMZ chemotherapy induce senescence in glioma cells in vitro and sensitize these cells to anti-apoptotic inhibitors with distinct efficacy according to the patient-derived cell lines, suggesting that this strategy may be relevant for some patients with GBM<sup>82,88</sup>. Based on our work which focused on naturally-occurring senescence, we hypothesize that senolytics combined with conventional treatment could have a double action by depleting both resident and therapy-induced senescent cells which may potentialize their effect. Also, we would like to propose another combined therapeutical strategy that would take advantage of the tumor ecosystem modifications following senolytic treatment. For example, promyelinating drugs may amplify the plasticity of malignant cells towards an oligodendroglial-like differentiated phenotype initiated by the senolytic treatment (Fig. 3). Another strategy could be to target the immune cells. Indeed, although a thorough study on the consequence of senolytic treatment on the immune system is required, our study provides evidence of a decreased anti-inflammatory phenotype after treatment (Fig. 4). Therefore, senolytic treatment may prime GBM to respond to immunotherapy. This hypothesis is attractive as the immunotherapies with the anti-PD1 and PD-L1 antibodies did not show an extension of the overall survival in treating patients with recurrent GBM<sup>89,90</sup>. One possible explanation for this failure could be that GBMs contain very few immune effector cells<sup>91</sup>. Further work on immunocompetent GBM models, reproducing the intra-heterogeneity of patient GBMs, is now needed to evaluate the effect of senotherapy on glioma progression and assess their efficacy as companion therapy.

## Methods

### Ethics

Fresh patient GBM samples were selected from the Pitié-Salpêtrière tumor bank Onconeurotek. They were reviewed by our senior
<figure><img src="image_60.png" /><figcaption>**b** p16<sup>INK4a</sup> Hi cells in NP-like and astrocyte clusters vs all malignant cells</figcaption></figure>

<figure><img src="image_61.png" /></figure>

<figure><img src="image_62.png" /></figure>

<figure><img src="image_63.png" /></figure>

<figure><img src="image_64.png" /></figure>

\#Events: 322; Global p-value (Log-Rank): 1.5414e-10

AIC:3138.05; Concordance Index:0.63

<figure><img src="image_65.png" /></figure>

<figure><img src="image_66.png" /></figure>

\#Events: 322; Global p-value (Log-Rank): 7.4205e-09

AIC:3146.26; Concordance Index:0.62

pathologist (F.B.) to validate the histological features and confirm patients' diagnoses. Collection of tumor samples and clinical-pathological information were obtained upon patients' informed consent and ethical board approval, as stated by the Declaration of Helsinki. The ethical approval was obtained from the ethical committee «CPP-Île de France VI-Groupe Hospitalier Pitié Salpêtrière». Molecular characterizations were performed as previously described⁹².

All animal care and treatment protocols complied with European legislation (no. 2010/63/UE) and national (French Ministry of Agriculture) guidelines for the use and ethical treatment of laboratory animals. All experiments on animals were approved by the ethical committee in animal experimentation Charles Darwin no. 5, Paris (approval APAFIS 9131). In the approved document the maximal tumor burden, as defined as end-points in this study, was determined by a loss of 10-15% of the mouse body weight in the 48 h interval and/or the
**Fig. 7 |** Mouse senescent signature is conserved in patient GBM and its enrichment score is predictive of worse survival. a Timeline of the mouse GBM generation for scRNAseq at the early timepoint. b Volcano plot of differentially expressed (DE) genes (−0.5 < log2FC > 0.5; FDR < 0.05) between p16<sup>INK4a</sup>HI cells (gene expression ≥ 4) of astrocyte and NP-like clusters compared with the remaining malignant cells in WT+GCV GBMs. c Heatmap of the 31 senescence signature genes in WT+GCV GBMs. d Top: Violin plots of the single-sample GSEA (ssGSEA) senescent Z-score in all patient GBM cells. Patient GBMs data were extracted from Bhaduri et al.<sup>26</sup>, Johnson et al.<sup>54</sup>, and Neftel et al.<sup>23</sup>. Bottom: Barplots of the percentage of the ssGSEA senescent Z-score distribution rate in all patient GBM cells. High and Low distribution rates correspond to the highest and lowest decile, respectively. e Table representing a Cox regression model using the ssGSEA-senescence score (sen-score), p16<sup>INK4a</sup> copy number alteration (p16-CNA) stratified

into a group without alteration (normal) and a group harboring homozygous
recessive deletion (homdel_rec) in the *INK4a* locus, the sex, the age of the patients
and the Karnofsky score. f Representative SA-β-gal staining (blue) coupled with IHC
(brown) and counterstained with hematoxylin (H) on patient GBM cryosections.
Three patient GBMs were analyzed per antibody. Scale bar: 10 µm. g Table repre-
senting a Cox regression model using the ssGSEA NRF2 targets score (NRF2 targets
score), p16<sup>INK4a</sup> copy number alteration, the sex, the age of the patients, and the
Karnofsky score. e and g Data were extracted from The Cancer Genome Atlas
(TCGA) GBM data sets and statistical significance was determined by a Log-Rank
test. The error bars correspond to 95% Confidence Intervals (CIs). GCV ganciclovir,
TMX tamoxifen, i.p. intraperitoneal, Iv lentivirus, Iv-luc lentivirus-luciferase, OS
overall survival. Raw data are provided as a Source Data file.

degradation of the general condition of the animal such as prostration/
decubitus, loss of alertness, cutaneous lesions.

Mouse breeding

To generate the GBM mouse model, we crossed Glast<sup>creERT2/+</sup> mice<sup>93</sup> with the Pten<sup>fl/fl</sup> mice<sup>94</sup>. Glast<sup>creERT2/+</sup>; Pten<sup>fl/fl</sup> males were bred with either Pten<sup>fl/fl</sup> or Pten<sup>fl/fl</sup>; p16-3MR/+females<sup>17</sup> to generate Glast<sup>creERT2/+</sup>; Pten<sup>fl/fl</sup> and Glast<sup>creERT2/+</sup>; Pten<sup>fl/fl</sup>; p16-3MR/+mice, named WT and p16-3MR mice respectively. All animals used in the study were from mixed genetic backgrounds C57BL/6J and Swiss and they were 6-8-week-old females except for the mice used for scRNAseq at the early timepoint that were 14-week-old females.

Plasmid construction

H-RasV12-shp53-(IRES)-GFP-(2a)-Firefly-luciferase vector was gener-
ated from the H-RasV12-shp53-(IRES)-GFP³⁷ construct using the Gibson
Assembly technique⁹⁵. The terminal IRES-GFP region of the initial
vector and a P2A-luciferase2 sequence were both flanked with a shared
sequence overlap and amplified by PCR. They were then inserted into
the Sall and Pmll sites of the initial vector.

Four oligonucleotide sequences for the miR-NRF2-based shRNAs
targeting have been designed using the Block-iT RNAi designer tools
(Invitrogen; Supplementary Data 5) and cloned into the pcDNA6.2-GW
/EmGFPmiR plasmid according to the manufacturer protocol (BLOCK-
iT polII miR RNAi, invitrogen #K4936-00). miR-NRF2 #4 and a miR-ctl⁹⁶
(Supplementary Data 5) have been further cloned by the Gibson
Assembly technique in the H-RasV12-shp53-(IRES)-GFP-(2a)-Firefly-
luciferase vector.

Stereotaxic injection

We stereotaxically performed lentiviral intracranial injection of mice to
induce de novo tumorigenesis. The mice were anesthetized with iso-
flurane (2-3%, 1 L/min oxygen), and subcutaneously injected in the
head with lidocaine (60 µL, 2.133 mg/mL). Analgesia was injected
intraperitoneally (i.p.) during presurgery and up to 24 h after surgery
(buprenorphine, 100 µL, 15 µg/mL). The HRasV12-shp53-(IRES)-GFP (lv)
or the HRasV12-shp53-(IRES)-GFP-(2a)-Firefly-luciferase lentivirus (lv-
luc, 1 µL, 6 × 10⁸ PFU/mL) was injected in the right subventricular zone
(SVZ) of the brain (x = 1 mm, y = 1 mm, z = -2.3 mm from the bregma).
We used a Hamilton 30 G needle with a silica fiber tip (MTI-FS) and an
automatic injector (Harvard Apparatus). After injection, the skin wound
was closed with surgical glue (SurgiBond®) and animals were placed
under an infrared lamp until they recover a vigil state. From the next
day, mice were injected i.p. with tamoxifen (TMX, 20 mg/mL in corn oil,
Sigma #T5648-IGi and Sigma #C8267) once per day for 5 consecutive
days to induce the recombination of the *Pten* locus and of the *loxP-RFP- loxP* cassette of the lentivirus allowing the expression of *H-RasV12*.

Bioluminescence imaging

We monitored tumor growth by in vivo bioluminescence twice a week
from 14 days post intracranial injection. The mice were i.p. injected

with Xenolight n-Luciferin (100 µL, 30 mg/mL, Perkin Elmer #122799),
anesthetized with isoflurane and their head and back were shaved.
Bioluminescence was recorded with an IVIS Spectrum In Vivo Imaging
System (Perkin Elmer) and the ratio was measured by normalizing the
head signal on the back signal. The onset of tumor growth corresponds
to a head/trunk bioluminescence ratio of 2 (see below) for the p16-3MR
+vhc and p16-3MR+GCV mice and to a head bioluminescence signal of
3e10⁶ for miR-ctl and miR-NRF2 GBM-bearing mice. The difference in
the evaluation of tumor growth was due to a point mutation in the P2A
sequence in the HRasV12-shp53-(IRES)-GFP-(2a)-Firefly-luciferase
vector.

Mouse treatments

Mice were treated with vehicle (PBS, DMSO 20%, Sigma #D8418-50ML) or GCV (25 mg/kg/day, Selleckchem, #S1878) prepared in PBS, 20% DMSO at 21 days post injection (DPI). During the course of the study, we implemented bioluminescence-monitored GBM growth for the two paradigms p16-3MR+vhc vs. p16-3MR+GCV and the WT+vhc vs. WT +ABT263 (see below). The mice were treated when head to back bioluminescence ratio was superior or equal to 2 (around 24 DPI; *n* = 43). GCV was administered via daily i.p. injections for 5 consecutive days per cycle, for two cycles with a 2-week interval between the two cycles. ABT263 (Selleckchem, #S1001) was prepared as previously described¹³ and was administered to mice by gavage at 50 mg/kg/day for 5 days.

Kaplan–Meier mice survival studies

Kaplan–Meier survival analysis was done using Prism (Graphpad soft-
ware v.8.2.1). In accordance with EU guidelines, mice were sacrificed
when reaching endpoints (20% body weight decrease, deterioration of
general condition). Mice were injected by batch. One batch always
included control and experimental mice injected the same day. When
control mice survival extended more than 57 DPI, the entire batch was
removed from the analysis to exclude technical bias linked to intra-
cranial injection.

Mice brain collection

When reaching endpoints, mice were sedated with CO₂ inhalation
followed by intracardiac perfusion with cold HBSS 1X. After harvesting
the brain, the GFP+ tumor was cut into two parts, under the MZFL II
stereomicroscope (Leica). The anterior part of the GFP+ tumor and the
GFP- parenchyma were chopped and stored in TRI-reagent (Molecular
Center Research, #TR 118) at -80 °C or directly snap-frozen in liquid
nitrogen for RNA isolation. The posterior part was snap-frozen in dry
ice-cooled-isopentane for histological studies. Brains were cryosec-
tioned at a 12-µm thickness (Leica cryostat).

SA-β-gal and immunohistochemical staining

For SA-β-gal staining, GBM sections were fixed in 2% PFA and 0.02%
glutaraldehyde (Sigma, #340855) for 10 min at RT. Note, that the
concentration of glutaraldehyde was dropped from 0.2% to 0.02%
from the original protocol⁹⁷ to allow combined immunohistochemical
staining. Sections were washed twice in PBS pH 7.0 and once in PBS pH 5.5 for 30 min. Slides were incubated in the X-gal solution as previously described⁹⁷ for 5 h 30 min at 37 °C for mouse sections and overnight (O/N) for patient GBM sections. Slides were then washed in PBS and post-fixed in 4% PFA for 10 min at RT.

For Immunohistochemical staining, GBM sections were fixed in 4% PFA for 10 min and washed in PBS. Endogenous peroxidases were inactivated in 1% H₂O₂ (in H₂O) solution for 5 min and sections were incubated in the blocking solution (PBS 1X, 10% NGS, 3% BSA, and 0.25–0.5% Triton) for 30 min. Sections were then incubated with the primary antibody (Supplementary Data 5) in the blocking solution for either 2 h at RT or O/N at 4 °C. Slides were rinsed and incubated with biotinylated secondary antibodies for 45 min at RT. An amplification step was performed using VECTASTAIN Elite ABC HRP Kit (Vector Laboratories, #PK-6100-NB) for 30 min at RT and staining was revealed by a DAB reaction. Images were acquired using an Axio Scan.ZI (Zeiss) and extracted using the ZEN 2.0 blue edition (Zeiss) software.

## Surface area quantification

Quantifications were performed using Fiji software (v.2.1.0/1.53c)⁹⁸. Region of interest (ROI) corresponding to the tumor, was selected using the ellipse tool. IHC images were then color deconvoluted according to the “Giemsa” or “Hematoxylin and DAB (H DAB)” vector to assess a threshold of the SA-β-gal or DAB signal, respectively. The signal threshold was adjusted in order to remove the unspecific background signal without clearing the specific one. Number of pixels was measured and the values were normalized on the GFP+ tumor surface area. For mice brain tissue, four slides with three sections on each (n = 12) for SA-β-gal quantification and three slides with three sections on each (n = 9) for IHC quantification were analyzed per sample. For the patient sample, four sections were analyzed for SA-β-gal quantification.

## Western blotting

Total proteins were extracted from tumor samples following the TRI-reagent protocol (Molecular Center Research, #TR 118). Protein pellets were solubilized in 1% SDS, 10 M urea, and stored at –80 °C. Protein concentration was assayed using the Pierce BCA protein assay kit (Thermo Fisher, #23225). Proteins were separated on 4–20% stain-free polyacrylamide gels (Mini-PROTEAN TGX Protein Gels, Biorad, #4568096) and transferred on a nitrocellulose membrane 0.45 µm (Thermo Fisher, #88018). Membranes were probed with primary antibodies (Supplementary Data 5) diluted in Super Block Blocking buffer in TBS (Thermo Fisher, #37535) and incubated O/N at 4 °C under gentle agitation. The secondary antibodies were incubated for 1 h at RT. Fluorescence was detected using the Odyssey CLx (Li-cor), and specific bands were quantified using Fiji software (v.2.1.0/1.53c)⁹⁸ and normalized against the corresponding β-TUBULIN band.

## Cell culture

Glioma 261 murine cell line (GL261) was cultured in DMEM (Thermo Fisher #31966021) 1% fetal bovine serum (Thermo Fisher #A3160801). The GL261 cell line was not authenticated. These cells were transfected with the pcDNA6.2-GW /EmGFPmiR plasmids containing miR sequences (miR-NRF2 #1, #2, #3, #4 and miR-ctl²⁶ using the FUGENE HD transfection reagent (Promega #E2311). Seven days later GFP positive cells were isolated by flow cytometry (Biorad S3e cell sorter), cultured for 2 more days and their total RNAs were extracted.

## RNA extraction and RT-qPCR

Total RNAs were extracted from tumor and parenchyma samples and GL261 cells following either the TRI-reagent (Molecular Center Research, #TR 118), the Maxwell RSC simplyRNA Tissue (Promega,

#AS1340), and the Macherey-Nagel Mini kit Nucleospin protocol (Macherey-Nagel, #740955.50).

cDNAs were generated using the Maxima Istr cDNA Synth Kit (LifeTechnologies, K1642). Quantitative PCR was performed using LightCycler 480 SYBR Green I Master Mix (Roche, #4707516001) on a LightCycler® 480 Instrument II (Roche). Samples were run in duplicate or triplicate, transcript levels were normalized to TBP and GAPDH, and analysis was performed using the 2⁻ΔΔCT method⁹⁹. Primers used in this study are listed in Supplementary Data 5.

## Bulk RNA-seq and analysis

The quantity and quality of the total RNAs extracted were assessed by the Tapestation 2200 (Agilent) and sequenced with the Illumina NextSeq 500 Sequencing system using NextSeq 500/550 High Output Kit v2 (150 cycles, # 20024907), 400 millions of reads, 50Gbases.

Quality of raw data was evaluated with FastQC (v.0.11.5). Poor quality sequences were trimmed or removed with Fastp software to retain only good-quality paired reads. Star (v2.5.3a) was used to align reads on mm10 reference genome using default parameters except for the maximum number of multiple alignments allowed for a read which was set to 1. Quantification of gene and isoform abundances was done with rsem (v1.2.28) on RefSeq catalog, prior to normalization with edgeR Bioconductor package (v3.28.0). Finally, differential analysis was conducted with the glm framework likelihood ratio test from edgeR. For malignant samples, a batch effect was detected in PCA representation. To correct it, we performed the analysis by using an additive model which includes this batch variable. Multiple hypothesis-adjusted *p*-values were calculated with the Benjamini–Hochberg procedure to control FDR.

Functional enrichment analysis was performed with clusterPro-filer (v3.14.3) Bioconductor package on the differentially deregulated genes with over-representation analysis (enricher function) and on all the genes with Gene Set Enrichment Analysis (GSEA function¹⁰⁰). Hallmark, Transcription factor targets (TFT), and Canonical pathways (CP) gene sets from MSigDB collections have been used, complete with some custom gene sets (Supplementary Data 1). CIBERSORT¹⁰¹ was used to accurately quantify the relative abundances of six distinct immune cell types according to the ImmGen immune cell genes sig-nature (reference GSE124829).

## Tumor dissociation for scRNAseq

After brain harvest, GFP+ tumors were dissected under a Leica MZFL II stereomicroscope. Tumor pieces were chopped and incubated for 5 min at 37 °C in an HBSS-papain-based lysis buffer (Worthington PAP) containing DNAse (0.01%, Worthington #LS002139) and L-Cystein (124 µg/mL, Sigma #C78805). Papain digestion was inhibited by ovo-mucoid (7 mg/mL, Worthington #LS003085). Tissue was further dissociated mechanically and centrifuged 300 × g, 10 min at 4 °C. Cells were resuspended in cold HBSS, a debris removal step was performed (Miltenyi #130-109-398) and blood cells were removed using a blood lysis buffer (Roche 11814 389001). After centrifugation, cells were resuspended in cold HBSS and incubated with the eBiosciences Fixable Viability Dye Fluor 450 or 660 (Invitrogen 65-0863), to label dead cells, and washed. Cells were then sorted using the MoFlo Astrios cell sorter (Beckman Coulter) or the S3e cell sorter (Biorad). Live cells were collected in HBSS 0.1% BSA precoated tubes, centrifuged, and resuspended in HBSS–0.1% BSA at a concentration of 1200 cells/µL. GFP+ and GFP– cells were collected separately for the scRNAseq ddSeq experiment.

## Single-cell RNA sequencing and analysis–ddSeq data

Cell suspension of one dissociated GBM was loaded in 4 wells (3 wells with GFP+ cells or malignant cells and 1 well with GFP– cells or non-malignant cells) on the ddSEQ Single-Cell Isolator (Biorad). A library was generated using SureCell™ Whole Transcriptome Analysis 3'
Library Prep Kit for the ddSEQ System (Illumina, #20014280) and was sequenced on a Nextseq 500 Illumina sequencing system, using a High Output Kit (150 cycles), with the following parameters: 400 million reads depth, 50 Gbases, and 70 million reads per sample.

Cutadapt 1.18 was used to trim nextera adapters in 3' on reads, then a quality control of sequences was done with FastQC. Cellular and UMIs barcodes were extracted with the ddSeeker (v 0.9.0) tool with default parameters. The following steps were done with Drop-seq tool (v 2.0.0). Trimming of 5' adapter sequences and of polyA tails was performed. Unaligned BAM was transformed to fastq with the Picard tool, prior to alignment with STAR on mm10 reference genome. Ddseeker bam outputs previously tagged with molecular/cell barcode were merged with aligned BAM files, according to the Drop-seq tool cookbook. Finally, TagReadWithGeneExonFunction was used to annotate each read with the gene it belongs to, and DigitalExpression was used to count gene transcripts in each cell. The output DGE matrix file is a matrix with a row for each gene, and a column for each cell, containing the number of transcripts observed. This output was loaded into the Seurat (v2) R package for further analysis keeping only cells where at least 200 features were detected, and genes detected in at least 5 cells. The final dataset contains 1740 cells and 15,448 genes. To normalize the data, we applied the global-scaling normalization method "LogNormalize" which normalizes the feature expression measurements for each cell by the total expression, multiplies this by a scale factor (10,000 by default), and log-transforms the result. Then, highly variable genes were detected prior to scaling transformation. To cluster cells, we computed a principal components analysis (PCA) on scaled variable genes, as determined above, using Seurat's *RunPCA* function, and visualized it by computing a Uniform Manifold Approximation and Projection (UMAP) using Seurat's *RunUMAP* function on the top 10 PCs. We used *FindClusters* function with a resolution of 0.6 resulting in 8 clusters. TME and malignant cluster cells were identified according to the expression of the GFP transgene. Then, the *FindAllMarkers* function was used to extract the top differentially expressed genes of each cluster, and to annotate them.

## Single-cell RNA-seq and analysis–10X data

Cells suspension of four dissociated GBMs (2 WT+GCV and 2 p16-3MR +GCV) were loaded with the Chromium Next GEM Chip G Single Cell Kit (10X Genomics, #PN-1000120) and a library was generated using Chromium Next GEM Single Cell 3' Reagent Kits v3.1 (10X Genomics, #20012850). The library was sequenced on an Illumina NovaSeq 6000 instrument using a 100-cycle S2 flow cell in XP mode, with the following parameters: 2050 million reads depth, 200 Gbases per run, and 50,000 reads per cell.

The Cell Ranger Single-cell Software suite (v3.0.2) was used to process the data. First, a custom reference genome was created with the `mkref` function to include 3'LTR and 3MR sequences into the mm10 reference genome. Count function was used on each GEM well that was demultiplexed by `mkfastq` to generate gene-cell matrices. Then, filtered feature bc matrix output was loaded into the Seurat Bioconductor package (v.3.2.3) to filter the datasets and identify cell types using R (v.3.6). Genes expressed in at least five cells and cells with at least 200 features were retained for further analysis. To remove likely dead or multiplet cells from downstream analyses, cells were discarded when they had <500 unique molecular identifiers (UMIs), >60,000 UMIs, or expressed over 8% mitochondrial genes.

All samples were merged together for downstream analysis. As no batch effects were observed among the four samples, no integration step was performed. Gene expression matrix was normalized using the negative binomial regression method implemented in the Seurat *SCTransform* function, via the selection of the top 3000 variable genes and regressed out the mitochondrion expression percentage. The final dataset was composed of 20,293 genes and 26,237 cells.

To cluster cells, we computed a principal components analysis (PCA) on scaled variable genes, as determined above, using Seurat's *RunPCA* function, and visualized it by computing a Uniform Manifold Approximation and Projection (UMAP) using Seurat's *RunUMAP* function on the top 30 PCs. We also computed the *k*-nearest neighbor graph on the top 30 PCs, using Seurat's *FindNeighbors* function with default parameters, and in turn, used Seurat's *FindClusters* function with varying resolution values. We chose a final value of 0.5 for the resolution parameter at this stage of clustering. Clusters were assigned preliminary identities based on the expression of combinations of known marker genes for major cell types. TME clusters were identified with the expression of the *Ptprc* (*Cd45*) gene marker. In order to better identify other cell types, TME cells were removed and a second clustering with a resolution of 0.6 was applied.

The *FindMarkers* function with the default parameters (min.-LogFC = 0.25, min.pct = 0.25, test.use = Wilcox) was used to identify differentially expressed genes in different conditions: (i) p16-3MR +GCV vs. WT+GCV in each cluster; (ii) cells from astrocyte and NP clusters with *Cdkn2a* expression ≥4 (307 cells) vs all the other cells (10,280 cells) in WT+GCV GBMs. Functional enrichment analysis was done with clusterProfiler (v3.14.3) Bioconductor package on the differentially deregulated genes with over-representation analysis (enricher function) and with Gene Set Enrichment Analysis (GSEA function¹⁰⁰) on all the genes. We searched for ligand/receptor interactions between cluster 0 et *Cd45* positive clusters at 0.5 resolution in our single-cell data, using CellPhoneDB (v.2.1.4). Copy number variations (CNVs) were inferred with inferCNV package (v.1.6.0) with the following parameters: "denoise" and a value of 0.1 for "cutoff".

## Signature expression analyses

We analyzed the senescence signature through our tumoral SCT normalized dataset and three datasets corresponding to patient GBMs²³,²⁶,⁵⁴. The Neftel et al., dataset was processed via Seurat (v3.2.3), 10X samples were normalized via the SCT method and Smartseq2 samples were retrieved in log2(TPM + 1). Fastq files for the Bhaduri et al., 10X dataset were processed by cellranger, followed by Seurat analysis with the SCT method. Finally, we retrieved the normalized expression matrix of Johnson and colleagues via synapse (https://www.synapse.org/#!Synapse:syn22257780/wiki/604645). For both murine and patient GBM datasets, we filtered out transcriptomes expressing the *CD45* (*PTPRC*) gene and pediatric GBMs²³ and we calculated a senescence score resulting from the single-sample GSEA¹⁰² using the R package GSVA (v.1.40.1). For these three datasets, we computed the z-scores of the resulting enrichment scores and sliced the signature score distribution into deciles to determine the HIGH senescence cells (last decile), the LOW senescence cells (1st decile), and the others with an average senescence potential (MEDIUM).

## Cox regression analysis

Normalized intensities from TCGA microarray data were obtained from cBioPortal (cbioportal.org), filtering for GBM TCGA, Firehose Legacy dataset. First, single-sample GSEA scores were calculated with R package GSVA (v1.32.0) for senescence genes signature and NRF2 targets signature. Secondly, we fitted a Cox proportional hazards regression model with the *coxph* function from the survival R package (v2.44.1-1), with additional covariates such as p16 copy number alteration (CNA) status, and age of patients. Plots were done with *ggforest* R function.

## Statistical analysis

Data are presented as mean with standard error to the mean (SD) unless otherwise specified. Statistical comparisons were performed using Wilcoxon-Mann-Whitney, *p*-values unless otherwise specified (*p* < 0.05; **p** < 0.01; ***p** < 0.001, **p** < 0.001). For Kaplan-Meier
survival curves, statistical significance was determined by the Mantel-Cox log-rank test (*p* < 0.05). Comparisons between conditions were performed using log-rank tests with a Benjamini-Hochberg correction for multiple comparisons.

## Reporting summary

Further information on research design is available in the Nature Portfolio Reporting Summary linked to this article.

## Data availability

The mm10 reference genome was retrieved from https://hgdownload.cse.ucsc.edu/goldenpath/mm10/chromosomes/. Gene Set Enrichment analysis gene sets came from MSigDB collections (https://software.broadinstitute.org/cancer/software/gsea/wiki/index.php/MSigDB_collections). The Neftel et al., the dataset was retrieved via the single cell portal (singlecell.broadinstitute.org), fastq files for the Bhaduri et al., dataset, were retrieved from SRA bioproject PRJNA579593, and the normalized expression matrix of Johnson and colleagues was retrieved via synapse (https://www.synapse.org/#!Synapse:syn22257780/wiki/604645). Normalized intensities from TCGA microarray data were obtained from cBioPortal (cbioportal.org), filtering for GBM TCGA, Firehose Legacy dataset. The raw data generated in this study are provided in a Source data file and a Source supplementary data file. The data generated in this study have been deposited in the Gene Expression Omnibus database under the accession code GSE168040. All data are available in the article, Supplementary Information and Source Data. Further information and material requests should be addressed to Isabelle Le Roux (isabelle.leroux@icm-institute.org). Source data are provided with this paper.

## Code availability

Computer codes used for the analysis of the senescence score are available following https://github.com/bellenger-l/glioblastoma_senescence/; https://doi.org/10.5281/zenodo.7525302.

## References

1. Ostrom, Q. T. et al. CBTRUS statistical report: primary brain and central nervous system tumors diagnosed in the United States in 2006–2010. *Neuro Oncol.* **15**, ii1–ii56 (2013).

2. Stupp, R. et al. Radiotherapy plus concomitant and adjuvant temozolomide for glioblastoma. *N. Engl. J. Med.* **352**, 987–996 (2005).

3. Gorgoulis, V. et al. Cellular senescence: defining a path forward. *Cell* **179**, 813–827 (2019).

4. Coppé, J.-P., Desprez, P.-Y., Krtolica, A. & Campisi, J. The senescence-associated secretory phenotype: the dark side of tumor suppression. *Annu. Rev. Pathol. Mech. Dis.* **5**, 99–118 (2010).

5. Faget, D. V., Ren, Q. & Stewart, S. A. Unmasking senescence: context-dependent effects of SASP in cancer. *Nat. Rev. Cancer* **19**, 439–453 (2019).

6. Sieben, C. J., Sturmlechner, I., van de Sluis, B. & van Deursen, J. M. Two-step senescence-focused cancer therapies. *Trends Cell Biol.* **28**, 723–737 (2018).

7. Xue, W. et al. Senescence and tumour clearance is triggered by p53 restoration in murine liver carcinomas. *Nature* **445**, 656–660 (2007).

8. Kang, T. W. et al. Senescence surveillance of pre-malignant hepatocytes limits liver cancer development. *Nature* **479**, 547–551 (2011).

9. Ruscetti, M. et al. NK cell-mediated cytotoxicity contributes to tumor control by a cytostatic drug combination. *Science* (80-.) **362**, 1416–1422 (2018).

10. Yoshimoto, S. et al. Obesity-induced gut microbial metabolite promotes liver cancer through senescence secretome. *Nature* **499**, 97–101 (2013).

11. Toso, A. et al. Enhancing chemotherapy efficacy in pten-deficient prostate tumors by activating the senescence-associated anti-tumor immunity. *Cell Rep.* **9**, 75–89 (2014).

12. Ruhland, M. K. et al. Stromal senescence establishes an immunosuppressive microenvironment that drives tumorigenesis. *Nat. Commun.* **7**, 1–18 (2016).

13. Chang, J. et al. Clearance of senescent cells by ABT263 rejuvenates aged hematopoietic stem cells in mice. *Nat. Med.* **22**, 1–9 (2016).

14. Zhu, Y. et al. The Achilles' heel of senescent cells: from transcriptome to senolytic drugs. *Aging Cell* **14**, 644–658 (2015).

15. Kohlhapp, F. J. et al. Venetoclax increases intratumoral effector T cells and antitumor efficacy in combination with immune checkpoint blockade. *Cancer Discov.* **11**, 68–79 (2021).

16. Baker, D. J. et al. Clearance of p16INK4a-positive senescent cells delays ageing-associated disorders. *Nature* **479**, 232–236 (2011).

17. Demaria, M. et al. An essential role for senescent cells in optimal wound healing through secretion of PDGF-AA. *Dev. Cell* **31**, 722–733 (2014).

18. Demaria, M. et al. Cellular senescence promotes adverse effects of chemotherapy and cancer relapse. *Cancer Discov.* **7**, 165–177 (2017).

19. Wang, H. et al. Targeting interleukin 6 signaling suppresses glioma stem cell survival and tumor growth. *Stem Cells* **27**, 2393–2404 (2009).

20. Banasavadi-Siddegowda, Y. K. et al. PRMT5-PTEN molecular pathway regulates senescence and self-renewal of primary glioblastoma neurosphere cells. *Oncogene* **36**, 263–274 (2017).

21. Marumoto, T. et al. Development of a novel mouse glioma model using lentiviral vectors. *Nat. Med.* **15**, 110–116 (2009).

22. Balakrishnan, I. et al. Senescence induced by BMI1 inhibition is a therapeutic vulnerability in H3K27M-mutant DIPG. *Cell Rep.* **33**, 108286 (2020).

23. Neftel, C. et al. An integrative model of cellular states, plasticity, and genetics for glioblastoma. *Cell* **178**, 835–849.e21 (2019).

24. Wang, L. et al. The phenotypes of proliferating glioblastoma cells reside on a single axis of variation. *Cancer Discov.* **9**, 1708–1719 (2019).

25. Couturier, C. P. et al. Single-cell RNA-seq reveals that glioblastoma recapitulates a normal neurodevelopmental hierarchy. *Nat. Commun.* **11**, 3406 (2020).

26. Bhaduri, A. et al. Outer radial glia-like cancer stem cells contribute to heterogeneity of glioblastoma. *Cell Stem Cell* **26**, 48–63.e6 (2020).

27. Richards, L. M. et al. Gradient of developmental and injury response transcriptional states defines functional vulnerabilities underpinning glioblastoma heterogeneity. *Nat. Cancer* **2**, 157–173 (2021).

28. Verhaak, R. G. W. et al. Integrated genomic analysis identifies clinically relevant subtypes of glioblastoma characterized by abnormalities in PDGFRA, IDH1, EGFR, and NF1. *Cancer Cell* **17**, 98–110 (2010).

29. Wang, Q. et al. Tumor evolution of glioma-intrinsic gene expression subtypes associates with immunological changes in the microenvironment. *Cancer Cell* **32**, 42–56.e6 (2017).

30. Castellan, M. et al. Single-cell analyses reveal YAP/TAZ as regulators of stemness and cell plasticity in glioblastoma. *Nat. Cancer* **2**, 174–188 (2020).

31. Bhat, K. P. L. et al. Mesenchymal differentiation mediated by NF-kB promotes radiation resistance in glioblastoma. *Cancer Cell* **24**, 331–346 (2013).

32. Darmanis, S. et al. Single-cell RNA-Seq analysis of infiltrating neoplastic cells at the migrating front of human glioblastoma. *Cell Rep.* **21**, 1399–1410 (2017).
33. Schmitt, M. J. et al. Phenotypic mapping of pathological crosstalk between glioblastoma and innate immune cells by synthetic genetic tracing. *Cancer Discov.* CD-20-0219 https://doi.org/10.1158/2159-8290.cd-20-0219 (2020).

34. Chen, Z. et al. Crucial role of p53-dependent cellular senescence in suppression of Pten-deficient tumorigenesis. *Nature* **436**, 725–730 (2005).

35. Acosta, J. C. et al. Chemokine signaling via the CXCR2 receptor reinforces senescence. *Cell* **133**, 1006–1018 (2008).

36. Chien, Y. et al. Control of the senescence-associated secretory phenotype by NF-κB promotes senescence and enhances chemosensitivity. *Genes Dev.* **25**, 2125–2136 (2011).

37. Friedmann-Morvinski, D. et al. Dedifferentiation of neurons and astrocytes by oncogenes can induce gliomas in mice. *Science* (**80-**) **338**, 1080–1084 (2012).

38. Louis, D. N. et al. The 2021 WHO classification of tumours of the central nervous system: a summary. *Neuro-Oncology* **23**, 1231–1251 (2021).

39. Coppé, J. P. et al. Senescence-associated secretory phenotypes reveal cell-nonautonomous functions of oncogenic RAS and the p53 tumor suppressor. *PLoS Biol.* **6**, 2853–68 (2008).

40. Acosta, J. C. et al. A complex secretory program orchestrated by the inflammasome controls paracrine senescence. *Nat. Cell Biol.* **15**, 978–990 (2013).

41. Basisty, N. et al. A Proteomic Atlas of senescence-associated secretomes for aging biomarker development. *SSRN Electron.* J. 1–26 https://doi.org/10.2139/ssrn.3380253 (2019).

42. Weng, Q. et al. Single-cell transcriptomics uncovers glial progenitor diversity and cell fate determinants during development and gliomagenesis article single-cell transcriptomics uncovers glial progenitor diversity and cell fate determinants during development and Gliom. https://doi.org/10.1016/j.stem.2019.03.006 (2019).

43. Tirosh, I. et al. Single-cell RNA-seq supports a developmental hierarchy in human oligodendroglioma. *Nature* **539**, 309–313 (2016).

44. Hara, T. et al. Interactions between cancer cells and immune cells drive transitions to mesenchymal-like states in glioblastoma. *Cancer Cell* **39**, 779–792.e11 (2021).

45. Bowman, R. L. et al. Macrophage ontogeny underlies differences in tumor-specific education in brain malignancies. *Cell Rep.* **17**, 2445–2459 (2016).

46. Kuleshov, M. V. et al. Enrichr: a comprehensive gene set enrichment analysis web server 2016 update. *Nucleic Acids Res.* **44**, W90–W97 (2016).

47. Hiebert, P. et al. Nrf2-Mediated fibroblast reprogramming drives cellular senescence by targeting the matrisome. *Dev. Cell* **46**, 145–161.e10 (2018).

48. Rojo de la Vega, M., Chapman, E. & Zhang, D. D. NRF2 and the hallmarks of cancer. *Cancer Cell* **34**, 21–43 (2018).

49. Gilder, A. S. et al. The urokinase receptor induces a mesenchymal gene expression signature in glioblastoma cells and promotes tumor cell survival in neurospheres. *Sci. Rep.* **8**, 1–16 (2018).

50. Brösicke, N. & Faissner, A. Role of tenascins in the ECM of gliomas. *Cell Adhes. Migr.* **9**, 131–140 (2015).

51. Osswald, M. et al. Brain tumour cells interconnect to a functional and resistant network. *Nature* **528**, 93–98 (2015).

52. Hayashi, M. et al. Microenvironmental activation of Nrf2 restricts the progression of Nrf2-activated malignant tumors. *Cancer Res.* **80**, 3331–3344 (2020).

53. Frohlich, D. A., McCabe, M. T., Arnold, R. S. & Day, M. L. The role of Nrf2 in increased reactive oxygen species and DNA damage in prostate tumorigenesis. *Oncogene* **27**, 4353–4362 (2008).

54. Johnson, K. C. et al. Single-cell multimodal glioma analyses identify epigenetic regulators of cellular plasticity and environmental stress response. *Nat. Genet.* **53**, 1456–1468 (2021).

55. Inoue, D. et al. Accumulation of p62/SQSTM1 is associated with poor prognosis in patients with lung adenocarcinoma. *Cancer Sci.* **103**, 760–766 (2012).

56. Menegon, S., Columbano, A. & Giordano, S. The dual roles of NRF2 in cancer. *Trends Mol. Med.* **22**, 578–593 (2016).

57. Harris, I. S. & DeNicola, G. M. The complex interplay between antioxidants and ROS in cancer. *Trends Cell Biol.* **30**, 440–451 (2020).

58. Pölönen, P. et al. Nrf2 and SQSTM1/p62 jointly contribute to mesenchymal transition and invasion in glioblastoma. *Oncogene* **38**, 7473–7490 (2019).

59. Mitsuishi, Y. et al. Nrf2 redirects glucose and glutamine into ana-bolic pathways in metabolic reprogramming. *Cancer Cell* **22**, 66–79 (2012).

60. Denicola, G. M. et al. Oncogene-induced Nrf2 transcription promotes ROS detoxification and tumorigenesis. *Nature* **475**, 106–110 (2011).

61. Mirzaei, R. et al. Brain tumor-initiating cells export tenascin-C associated with exosomes to suppress T cell activity. *Oncoimmunology* **7**, 1–14 (2018).

62. Angel, I., Pilo Kerman, O., Roussou-Noori, L. & Friedmann-Morvinski, D. Tenascin C promotes cancer cell plasticity in mesenchymal glioblastoma. *Oncogene* **39**, 6990–7004 (2020).

63. Zhou, P. et al. CD151-3β1 integrin complexes are prognostic markers of glioblastoma and cooperate with EGFR to drive tumor cell motility and invasion. *Oncotarget* **6**, 29675–29693 (2015).

64. Yue, H., Xu, Q. & Xie, S. High EMP3 expression might independently predict poor overall survival in glioblastoma and its expression is related to DNA methylation. *Medicine* **97**, 1–5 (2018).

65. Yuan, Z. S., Cao, Y. & Li, Z. Y. IGFBP2 induces SPRY1 expression via NF-κB signaling pathway in glioblastoma multiforme (GBM). *Eur. Rev. Med. Pharmacol. Sci.* **21**, 5072–5080 (2017).

66. Sharanek, A. et al. Transcriptional control of brain tumor stem cells by a carbohydrate binding protein. *Cell Rep.* **36**, 109647 (2021).

67. Morita, T. & Hayashi, K. Tumor progression is mediated by Thymosin-b4 through a TGFb/MRTF signaling axis. *Mol. Cancer Res.* **16**, 880–893 (2018).

68. Shi, Q. et al. Targeting SPARC expression decreases glioma cellular survival and invasion associated with reduced activities of FAK and ILK kinases. *Oncogene* **26**, 4084–4094 (2007).

69. Friedmann-Morvinski, D. et al. Targeting NF-κB in glioblastoma: a therapeutic approach. *Sci. Adv.* **2**, e1501292 (2016).

70. Tagscherer, K. E. et al. Apoptosis-based treatment of glioblastomas with ABT-737, a novel small molecule inhibitor of Bcl-2 family proteins. *Oncogene* **27**, 6646–6656 (2008).

71. Karpel-Massler, G. et al. TIC10/ONC201 synergizes with Bcl-2/Bcl-xL inhibition in glioblastoma by suppression of Mcl-1 and its binding partners in vitro and in vivo. *Oncotarget* **6**, 36456–36471 (2015).

72. Karpel-Massler, G. et al. Combined inhibition of Bcl-2/Bcl-xL and Usp9x/Bag3 overcomes apoptotic resistance in glioblastoma in vitro and in vivo. *Oncotarget* **6**, 14507–14521 (2015).

73. Karpel-Massler, G. et al. Induction of synthetic lethality in IDH1-mutated gliomas through inhibition of Bcl-xL. *Nat. Commun.* **8**, 1–13 (2017).

74. Guerrero, A. et al. Cardiac glycosides are broad-spectrum senolytics. *Nat. Metab.* **1**, 1074–1088 (2019).
75. Triana-Martínez, F. et al. Identification and characterization of Cardiac Glycosides as senolytic compounds. *Nat. Commun.* **10**, 1–12 (2019).

76. L'Hôte, V. et al. Ouabain and chloroquine trigger senolysis of BRAF-V600E-induced senescent cells by targeting autophagy. *Aging Cell* **20**, 1–14 (2021).

77. Amor, C. et al. Senolytic CAR T cells reverse senescence-associated pathologies. *Nature* **583**, 127–132 (2020).

78. Moiseeva, O. et al. Metformin inhibits the senescence-associated secretory phenotype by interfering with IKK/NF-κB activation. *Aging Cell* **12**, 489–498 (2013).

79. Wyld, L. et al. Senescence and cancer: a review of clinical implications of senescence and senotherapies. *Cancers* **12**, 1–20 (2020).

80. Kaefer, A. et al. Mechanism-based pharmacokinetic/pharmacodynamic meta-analysis of navitoclax (ABT-263) induced thrombocytopenia. *Cancer Chemother. Pharmacol.* **74**, 593–602 (2014).

81. Schwarzenbach, C. et al. Targeting c-iap1, c-iap2, and bcl-2 eliminates senescent glioblastoma cells following temozolomide treatment. *Cancers* **13**, 3585 (2021).

82. Rahman, M. et al. Selective vulnerability of senescent glioblastoma cells to BCL-XL inhibition. *Mol. Cancer Res.* **20**, 938–948 (2022).

83. Musi, N. et al. Tau protein aggregation is associated with cellular senescence in the brain. *Aging Cell* **17**, e12840 (2018).

84. Zhang, P. et al. Senolytic therapy alleviates Aβ-associated oligodendrocyte progenitor cell senescence and cognitive deficits in an Alzheimer's disease model. *Nat. Neurosci.* **22**, 719–728 (2019).

85. Liu, Y. et al. Combining an aurora kinase inhibitor and a death receptor ligand/agonist antibody triggers apoptosis in melanoma cells and prevents tumor growth in preclinical mouse models. *Clin. Cancer Res.* **21**, 5338–5348 (2015).

86. Fleury, H. et al. Exploiting interconnected synthetic lethal interactions between PARP inhibition and cancer cell reversible senescence. *Nat. Commun.* **10**, 2556 (2019).

87. Saleh, T. et al. Clearance of therapy-induced senescent tumor cells by the senolytic ABT-263 via interference with BCL–XL–BAX interaction. *Mol. Oncol.* **14**, 2504–2519 (2020).

88. Aasland, D. et al. Temozolomide induces senescence and repression of DNA repair pathways in glioblastoma cells via activation of ATR–Chk1, p21, and NF-κB. *Cancer Res.* **79**, 99–113 (2019).

89. Litak, J., Mazurek, M., Grochowski, C., Kamieniak, P. & Roliński, J. PD-L1/PD-1 axis in glioblastoma multiforme. *Int. J. Mol. Sci.* **20**, 5347 (2019).

90. Reardon, D. A. et al. Effect of nivolumab vs bevacizumab in patients with recurrent glioblastoma: the CheckMate 143 Phase 3 Randomized Clinical Trial. *JAMA Oncol.* **6**, 1003–1010 (2020).

91. Quail, D. F. & Joyce, J. A. The microenvironmental landscape of brain tumors. *Cancer Cell* **31**, 326–341 (2017).

92. Di Stefano, A. L. et al. Clinical, molecular, and radiomic profile of gliomas with FGFR3-TACC3 fusions. *Neuro Oncol.* **22**, 1614–1624 (2020).

93. Mori, T. et al. Inducible gene deletion in astroglia and radial glia—a valuable tool for functional and lineage analysis. *Glia* **54**, 21–34 (2006).

94. Marino, S. et al. PTEN is essential for cell migration but not for fate determination and tumourigenesis in the cerebellum. *Develop-ment* **129**, 3513–3522 (2002).

95. Gibson, D. G. et al. Enzymatic assembly of DNA molecules up to several hundred kilobases. *Nat. Methods* **6**, 343–345 (2009).

96. Morán, I. et al. Human β cell transcriptome analysis uncovers lncRNAs that are tissue-specific, dynamically regulated, and abnormally expressed in type 2 diabetes. *Cell Metab.* **16**, 435–448 (2012).

97. Le Roux, I., Konge, J., Le Cam, L., Flamant, P. & Tajbakhsh, S. Numb is required to prevent p53-dependent senescence following skeletal muscle injury. *Nat. Commun.* **6**, 8528 (2015).

98. Schindelin, J. et al. Fiji: an open-source platform for biological-image analysis. *Nat. Methods* **9**, 676–682 (2012).

99. Livak, K. J. & Schmittgen, T. D. Analysis of relative gene expression data using real-time quantitative PCR and the 2-ΔΔCT method. *Methods* **25**, 402–408 (2001).

100. Subramanian, A. et al. Gene set enrichment analysis: a knowledge-based approach for interpreting genome-wide expression profiles. *Proc. Natl Acad. Sci. USA* **102**, 15545–15550 (2005).

101. Newman, A. M. et al. Robust enumeration of cell subsets from tissue expression profiles. *Nat. Methods* **12**, 453–457 (2015).

102. Barbie, D. A. et al. Systematic RNA interference reveals that oncogenic KRAS-driven cancers require TBK1. *Nature* **462**, 108–112 (2009).

103. Weng, Q. et al. Single-cell transcriptomics uncovers glial progenitor diversity and cell fate determinants during development and gliomagenesis. *Cell Stem Cell* **24**, 707–723.e8 (2019).

## Acknowledgements

The authors acknowledge the service of the ICM platforms: iGenSeq, iVector especially Phillippe Ravassard, icm-QUAN-Imaging, Pheno-ICMice, Histomics-Histology, Celis-Cell culture, the Data Analysis Core platform for their expertise in statistics, in particular, François-Xavier Lejeune and the service of the animal facility UMS28 (Sorbonne University) and the service of the flow cytometry platform CYTO-ICAN especially Florence Deknuydt (Hôpital Pitié-Salpétrière). The authors are grateful to Judith Campisi, Inder Verma, Magdalena Götz, Anton Berns, and the Centro Nacional de Investigaciones Oncológicas (CNIO) for providing reagents and Carol Schuurmans for comments on the manuscript. We also thank members of the Huillard-Sanson laboratory for discussions, and help with the experiments, in particular Sophie Paris for teaching the intracranial injection. This work was supported by institutional fundings from Paris Brain Institute (ICM), the Institut National de la Santé Et de la Recherche Médicale and the Centre National de la Recherche Scientifique and grants from the ATIP-AVENIR (E.H.), the Cancéropole Ile de France (Emergence Program, I.L.R.), the Ligue contre le cancer, comité Ile de France (I.L.R.); the Fondation ARC pour la Recherche sur le Cancer (E.H., I.L.R.); the SIRIC-CURAMUS (joint Emergence program I.L.R. and C.A.). R. Salam was supported by fellowships from the French Ministry of Education and Research and the Ligue Nationale Contre le Cancer; A. Saliou was supported by a fellowship from the Ligue Nationale Contre le Cancer.

## Author contributions

R.S.: Conceptualization, methodology, experiments, formal analysis, methodology, writing—original draft. A.S.: Conceptualization, methodology, experiments, formal analysis, methodology, writing—original draft. F.B.: Methodology, formal analysis, resources. M.B.: Formal analysis, bioinformatics. C.A.: Formal analysis, bioinformatics, funding acquisition, writing—original draft. C.C.: Resources. A.A.: Bioinformatics. L.C.: Resources. M.S.: Funding acquisition, resources. E.H.: Funding acquisition, resources. L.B.: Formal analysis, bioinformatics, writing—original draft. J.G.: Formal analysis, bioinformatics, writing—original draft. I.L.R.: Conceptualization, methodology, experiments, formal analysis, methodology, funding acquisition, supervision,

## Competing interests

F.B. reports employment of next-of-kin from Bristol-Myers Squibb; research grants from Sanofi and Abbvie outside the submitted work; travel, and accommodations expenses from Bristol-Myers Squibb for travel expenses, outside the submitted work. The remaining authors declare no competing interests.
**Additional information**

**Supplementary information** The online version contains supplementary material available at https://doi.org/10.1038/s41467-023-36124-9.

**Correspondence** and requests for materials should be addressed to Isabelle Le Roux.

**Peer review information** *Nature Communications* thanks the anonymous reviewers for their contribution to the peer review of this work. Peer reviewer reports are available.

**Reprints and permissions information** is available at http://www.nature.com/reprints

**Publisher's note** Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Open Access** This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this license, visit http://creativecommons.org/licenses/by/4.0/.

© The Author(s) 2023