Prof Marie-Agnes Dillies

Prof Marc

Dr Jean-Jacques Diaz

Dr Thierry Virolle

Dr Florence Cavalli

Dr Marie-Agnès Dillies

Prof Monika Hegi

Mirca S Saurty-Seerunghen

Thomas Daubon

Léa Bellenger

Virgile Delaunay

Gloria Castro

Joris Guyon

Ahmed Rezk

Sylvie Fabrega

Ahmed Idbaih

Jean-Michel Peyrin

Christophe Antoniewski

Hervé Chneiweiss

Elias E El-Habr

Marie-Pierre Junier

EGF: Epidermal Growth Factor EGFR: Epidermal Growth Factor Receptor ELOVL2: ELOVL Fatty Acid Elongase

Keywords: glioma, scRNA-seq, functioning states, metabolism, computational analyses Actinin Alpha 4 ALDOA: Aldolase, Fructose-Bisphosphate Aryl Hydrocarbon Receptor Nuclear Translocator Glioma, scRNA-seq, MST, persulfidation, invasion, migration, computational analyses

The diversity of cancer cells populating glioblastoma is a major challenge in treating these brain cancers, which are the most devastating form of malignant brain tumors in adults. Cells in the same functioning state are present in tumors with distinct genomic anomalies but also among genetic subclones within each tumor. These functioning states are crucial for sustaining tumor development, and contribute to poor patient prognosis. Downstream of all signaling pathways regulating cell functioning states is metabolism. Hence, therapeutic targeting of relevant metabolism elements should help overcome the inter-and intra-tumor genetic heterogeneity in glioblastoma.

Understanding the metabolic weaknesses of glioblastoma cells in link with their functioning state in fully-grown patient tumors as encountered at the time of diagnosis, opens a path to overcome tumor heterogeneity, and develop more effective therapies.

To highlight the metabolic underliers of glioblastoma cells in distinct functional states, we analyzed public single-cell RNA-sequencing data from glioblastoma surgical resections, which offer the closest available view of tumor cell heterogeneity as encountered at the time of patients' diagnosis. We combined computational analyses of unicellular transcriptomes, modeling of the metabolic pathways at play in specific glioblastoma cell functioning states, and validation of the in silico findings using experimental manipulation of in vitro and in vivo human glioblastoma models.

Unsupervised without a priori analyses proved unsuitable for identifying cell functioning states from single-cell transcriptomic data. Identity of each tumor is encoded by information dispersed throughout the cell transcriptome, and masks other biological information. We thus developed a data reduction approach based on experimentally-defined molecular signatures which overcame this hurdle, unmasked information related to cell functioning states and allowed us to reconstruct the metabolic pathways at play in distinct cell functioning states. We validated our analytical strategy on two cell functioning states, using independent datasets of glioblastoma cell transcriptomes, patient-derived cell lines and orthotopic xenografts. The first one is the most aggressive state of cancer cells, their tumorigenic state which corresponds to their ability to initiate tumors. We found that glioblastoma cells with high tumorigenic potential are characterized by an enrichment in lipid and amino acid metabolism. Gene expression network modeling highlighted the very long chain polyunsaturated fatty acid synthesis pathway at the core of the network, notably its most downstream enzymatic component, ELOVL Fatty Acid Elongase 2 (ELOVL2). We demonstrated that ELOVL2 is required for glioblastoma cell tumorigenicity in vivo and that it may act through the formation/release of extracellular vesicles mediating cell-cell communication. The second cell functioning state which is a key contributor to glioblastoma malignancy and poor prognosis is the motile state. Computational analyses integrating cell trajectory modeling highlighted motile vi glioblastoma cells characterized by enhanced oxidative stress coupled with mobilization of the enzyme 3-Mercaptopyruvate sulfurtransferase (MPST). We found that enhanced ROS production was instrumental for the cells' motility in vitro. We also demonstrated experimentally that expression and activity of MPST was required for glioblastoma cell motility. Biochemical assays indicated that MPST acts by protecting protein cysteine residues from dismal hyperoxidation. In vivo, alteration of GBM cell properties by MPST knockdown translated in reduced tumor burden, and a robust increase in mice survival. Taken together, these results show that enhanced oxidative stress coupled with MPST mobilization plays a key role in GBM cell motility.

Altogether, our results demonstrate that single-cell transcriptomes from patients' glioblastoma can be harnessed to obtain an overview of metabolic pathways at play, within patient tumors, in link with cell functioning states. They also support the potential of our methodological strategy for the discovery of key druggable elements of metabolic modules whose therapeutic targeting is likely to overcome the inter-and intra-tumor heterogeneity of genomic anomalies in glioblastoma.

I would like to thank all the people who have supported me to complete this PhD work. It has been a long journey, and I certainly wouldn't be where I am now if it wasn't for them. It was an incredibly interesting and fruitful experience that completely renewed my motivation for my work.

First of all, I would like to express my deepest gratitude to my project supervisor, Dr. Marie-Pierre Junier (Research director). She was involved in the project design, participated to the experimental work and spent much time on the analysis of the bio-informatics results as well as on the writing and correction of my articles and manuscript. She was an invaluable help, providing continuous support and advice throughout these four years.

I would also like to thank Prof. Hervé Chneiweiss (Research director, director of the NPS-IBPS laboratory) for his availability, encouragement and critical insights into my research work, which pushed me to sharpen my thinking and brought my work to a higher level.

I am also extremely grateful to Dr. Elias El-Habr (Associate Professor), who co-supervised the metabolic aspect of my research work. His expertise in metabolism and his continuous support and advice were instrumental for this work. In addition, he provided great support for the experiments, including writing of experimental protocols, their setup, results analysis and even participated to some experiments, notably FACS analyses.

I also wish to thank the other members of my team. Virgile Delaunay (Technician) took in charge the xenografts of GBM cells in mice, together with Dr. Marie-Pierre Junier, and ensured tumor growth follow-up by bioluminescence. He also helped with cell culture. Gloria Castro Davalos (Master2 student) was in charge of investigating the effect of MPST knockdown on cell migration in microfluidics chips, which she set up in collaboration with Jean-Michel Peyrin (Research director). I would also like to acknowledge the great support from former members of our team. In addition to constant encouragement, Dr. Ashwin Narayanan and Dr. Alexandra Bogeas (former post-doctoral research fellows) contributed to the acquisition of experimental skills such as cell culture, transduction, western blot, RT-Q-PCR among others. Dr. Ghislaine Morvan-Dubois (Researcher, former team member) encouraged me into acquiring bio-informatics skills. Stecy Mienanzambi (former BSc student) helped with several experiments, including plasmid preparation and amplification and cell motility assays.

During the weekly lab meetings, all team members helped in understanding the results, in analyzing them and in proposing hypotheses/suggestions to progress in the project. Thus, the contribution of each member has been precious for the completion of this project. Thank you all for the countless times you have helped me with work as well as for the great discussions we had together. iv I had the great pleasure of working with Dr. Christophe Antoniewski (Head of the ArtBio platform at the IBPS) and Léa Bellenger (Bio-informatician). I thank them for their invaluable support for the acquisition of bio-informatic skills, for their patience, for the numerous insightful discussions we had, and for their assistance in this research work.

I would like to thank other collaborators as well. I thank Dr. Soazig Lelay (Researcher) and Dr.

Ramaroson Andriantsitohaina (Research director) for sharing the Nanosight equipment as well as Dr.

Delphine Garnier (Associate Professor) for helping out with experiments regarding the evaluation of extracellular vesicle size and quantity. I thank Dr. Thomas Daubon (Researcher) for his help in the GBM cell motility project, notably during definition of the molecular signature and the setup of migration and invasion assays. His PhD student Joris Guyon took in charge the quantification of the extent of migration on Matrigel and invasion in collagen. I thank Dr. François-Xavier Lejeune (Biostatistician) who though was not a collaborator in my research work, was always available for all my questions and provided helpful advice and practical suggestions.

I also thank Vincent Junier (Graphic designer) for defining a palette of colors that are very distinguishable from one another, which proved helpful for enhancing the readability of the results from computational analyses. I would also like to thank all researchers that participated to Open Science, by making their data available to everyone. I would like to express my deepest appreciation to the members of my mid-thesis committee, Dr. Table 1:

List of Tables

Thesis Overview

Gliomas are the most frequent tumors developing from the cancerous transformation of neural cells of the central nervous system. Glioblastoma (GBM) is the most aggressive form of gliomas and the most common and aggressive form of primary brain tumors in adults. The median survival of patients does not exceed 20 months despite aggressive treatments combining surgical removal of the tumor with chemotherapy and radiotherapy. GBM is a paradigm of intra-tumor heterogeneity. In addition to the intermingling of normal and cancerous cells, the tumor tissue is characterized by the co-existence of cancer cells bearing differing mutational loads and genomic rearrangements as well as functionally divergent cancer cell populations (e.g. with or without stem-like properties, migratory or static, proangiogenic or not, tumorigenic or not, sensitive or resistant to therapies).

The aim of my thesis work was to explore, at the single-cell level, cell functioning states as present in patients' GBM at the time of diagnosis, and to identify their eventual metabolic weaknesses. The methodological approach combines computational analyses of single-cell RNA sequencing (scRNA-seq) datasets, modeling of the metabolic pathways at play in specific GBM cell functioning states, and validation of the in silico findings using experimental manipulation of in vitro and in vivo human GBM models.

The general introduction comprises therefore a first part including an overview of glioma, a description of GBM pathophysiological features and available experimental models, and a summary of the available knowledge on three GBM cell characteristics that defy therapies, namely the tumor microenvironment, GBM cell motility and GBM cell heterogeneity and plasticity. This part is completed with an overview of the knowledge available on GBM cell metabolism in link with cell functional states, in accordance with the focus of my research work on metabolism. With respect to the methodologies implemented in this thesis work, the second part of the introduction presents scRNA-seq techniques and introduces the general workflow of scRNA-seq computational analysis.

The experimental results section is divided into three parts. The core PhD project involves work on two GBM cell functioning states, tumorigenicity and motility, which are presented as a published article (part 1), and as a manuscript (part 3), respectively. The second part of the results section includes the manuscript presenting the work related to the tumor identity card, which is a follow-up of the first published article on tumorigenicity. I participated in parallel to works on a third GBM cell functioning state, cell senescence, within the frame of another study from our team on the links between metabolism and GBM cell functioning states (data not presented).

The thesis concludes with a discussion on the contributions of this PhD work and their relevance for the field of study, and suggests future research directions.

GENERAL INTRODUCTION

INTRODUCTION

18

I. GLIOBLASTOMA

I.1 Glioma

Gliomas are part of primary brain and other central nervous system (CNS) tumors, the latter constituting a heterogeneous set of neoplasms associated with variable behaviors, symptoms, origins and malignancy. Primary brain and other CNS tumors are among the most frequent cancer sites and also among the most common causes of cancer death (Figure 1A-B). CNS tumors can be classified based on their behavior as non-malignant (70.3%) or malignant (29.7%), or based on their histology (Figure 1C) [START_REF] Ostrom | CBTRUS Statistical Report: Primary Brain and Other Central Nervous System Tumors Diagnosed in the United States in 2013-2017[END_REF]. They accounted for 1.6% of cancer diagnoses and 2.5% of cancerrelated deaths worldwide in 2020 [START_REF] Sung | Global Cancer Statistics 2020: GLOBOCAN Estimates of Incidence and Mortality Worldwide for 36 Cancers in 185 Countries[END_REF].

Gliomas are the most common malignant primary CNS tumors in adults. They are a heterogeneous group of tumors. World Health Organization (WHO) initially classified gliomas into 4 grades of malignancy (grade I to IV) based on their histopathological properties (Figure 2) [START_REF] Louis | The 2007 WHO classification of tumours of the central nervous system[END_REF]. Grade I gliomas (e.g. pilocytic astrocytoma) are low-proliferative lesions that can be cured by surgical resection alone. Grade II gliomas (e.g. low-grade diffuse astrocytomas) are more infiltrative, often recur after surgical resection and may progress to higher grades of malignancy. Grade III gliomas (e.g. anaplastic astrocytoma) correspond to intermediate-high-grade lesions with

Box 1: Cell of origin in glioblastoma

The identity of the cell of origin of glioblastoma remains ambiguous, though several genetic and epigenetic alterations are known to initiate or promote gliomagenesis. Studies using rodent models have shown that multiple cell types could give rise to glioblastoma if genetically manipulated, including neural stem cells (NSC) [START_REF] Koso | Transposon mutagenesis identifies genes that transform neural stem cells into glioma-initiating cells[END_REF] [START_REF] Llaguno | Malignant astrocytomas originate from neural stem/progenitor cells in a somatic tumor suppressor mouse model[END_REF] and oligodendrocyte progenitor cells (OPC) [START_REF] Liu | Mosaic analysis with double markers reveals tumor cell of origin in glioma[END_REF] [START_REF] Sugiarto | Asymmetry-defective oligodendrocyte progenitors are glioma precursors[END_REF]. In a coherent manner, comparison of normal brain tissues located away from the tumor with the tumor tissue itself resulted in the identification in the neurogenic sub-ventricular zone of low level drivermutations similar to those found in the matching tumor (Lee et al 2018a). More mature cells like neurons and astrocytes have also been shown to give rise to glioblastoma [START_REF] Friedmann-Morvinski | Dedifferentiation of neurons and astrocytes by oncogenes can induce gliomas in mice[END_REF]. Noteworthily, cancerous transformation of mature astrocytes has been shown to involve a transition of the mature cell towards an immature state akin to Neural Progenitor/Stem cells (NPC/NSC) [START_REF] Dufour | Astrocytes reverted to a neural progenitor-like state with transforming growth factor alpha are sensitized to cancerous transformation[END_REF] [START_REF] Friedmann-Morvinski | Dedifferentiation of neurons and astrocytes by oncogenes can induce gliomas in mice[END_REF]. Altogether, data gathered by the numerous studies on the subject tend to indicate that glioblastoma are more likely to derive from neural cells endowed with either innate or acquired immature features. Nevertheless, the cell of origin of brain neoplasms will be a matter of speculation as long as the identification of a pre-neoplastic state in the brain remains impossible. histological malignant characteristics such as nuclear atypia and abrupt mitotic activity. Patients with such tumors are treated with adjuvant radiation and/or chemotherapy. Finally, grade IV gliomas (e.g. glioblastoma) are the most aggressive tumors with cytologically malignant, mitotically active, necrosis-prone, angiogenic and highly infiltrative lesions. These tumors usually evolve rapidly and are associated to a fatal outcome. The identity of their cell of origin remains ambiguous, though several genetic and epigenetic alterations are known to initiate or promote gliomagenesis (Box 1). Tumor classification aims at grouping patient tumors into meaningful subgroups so as to guide treatment protocols.

Recently, CNS tumor classification has evolved from a classification based only on histological and structural similarities with normal glial cells to a more molecular and genetic one in combination with phenotypic information [START_REF] Louis | The 2016 World Health Organization Classification of Tumors of the Central Nervous System: a summary[END_REF]. All diffusely infiltrating gliomas are now termed as diffuse gliomas and are grouped on the basis of their growth pattern, behavior, and mutation in Isocitrate dehydrogenase 1 or 2 (IDH1/2) genes (Figure 3).

I.2 Glioblastoma physio-pathological features and molecular classification

Glioblastoma (GBM, also called Glioblastoma multiforme) are the most frequent and most aggressive forms of primary CNS tumors in adults, accounting for 14.5% of primary CNS tumors and 48.6% of primary malignant CNS tumors) [START_REF] Ostrom | CBTRUS Statistical Report: Primary Brain and Other Central Nervous System Tumors Diagnosed in the United States in 2013-2017[END_REF]. The median survival of GBM patients is only 10-15 months after diagnosis [START_REF] Louis | The 2016 World Health Organization Classification of Tumors of the Central Nervous System: a summary[END_REF]. These tumors are highly heterogeneous, with their genetic, epigenetic and functional characteristics differing between patient tumors (inter-tumor heterogeneity) but also within the same tumor (intra-tumor heterogeneity). The resulting extensive heterogeneity is one of the main challenges to effective treatment in GBM.

WHO classifies GBM into IDH-wildtype or IDH wt (~90%), IDH-mutant or IDH mut (~10%) and Not

Otherwise Specified or NOS (Table 1) [START_REF] Louis | The 2016 World Health Organization Classification of Tumors of the Central Nervous System: a summary[END_REF]. IDH wt GBM corresponds to primary GBM appearing de novo without clinical and histological evidence of a precursor lesion. The median overall survival of patients is about 10-15 months. In contrast, IDH mut GBM develop from low-grade gliomas, in younger patients (median age of 44 years compared to 62 years for IDH wt GBM) and have a better prognosis (median survival of 24-31 months). GBM, NOS refers to tumors for which full IDH evaluation cannot be performed.

An additional GBM classification has been proposed based on bulk expression profiles obtained from tumor tissue fragments (Figure 4). Phillips and colleagues were the first to distinguish high-grade gliomas into proneural, mesenchymal and proliferative transcriptional subtypes based on expression of defined gene signatures [START_REF] Phillips | Molecular subclasses of highgrade glioma predict prognosis, delineate a pattern of disease progression, and resemble stages in neurogenesis[END_REF]. Later, Verhaak and colleagues proposed Treatment of several cancers like thyroid cancer, breast cancer or melanoma has significantly improved during the past two decades, resulting in a better outcome for patients (Dal Maso et al 2020). On the opposite, GBM remains a lethal cancer with high rates of treatment resistance and a survival of less than two years despite aggressive treatment protocol combining maximal safe debulking surgery, radiotherapy and concomitant or adjuvant chemotherapy, usually with Temozolomide (TMZ) [START_REF] Stupp | Radiotherapy plus concomitant and adjuvant temozolomide for glioblastoma[END_REF] (Stupp et al 2017). Several factors can explain GBM therapy failure besides the limitations of the GBM experimental models (see Box 2 for a summary of the currently available models). Among them, tumor microenvironment (TME), GBM cell invasiveness and cancer cell heterogeneity and plasticity appear as especially challenging to overcome. Each GBM presents with a unique TME, composed of malignant cells interacting with the brain microenvironment (immune cells, neural cells, endothelial cells, …), thus forming a complex niche.

The non-malignant cells can influence progression of GBM tumors, as well as their ability to resist current therapies. Moreover, GBM cells are highly infiltrative in nature, making complete surgical resection impossible. Another major obstacle to effective treatment is the extensive heterogeneity displayed by GBM cells at different levels: genetic, epigenetic, transcriptional or functional between tumors but also within the same tumor. GBM cells are also highly plastic -they are able to adapt, and survive, to different microenvironmental cues and to therapies. These challenges are discussed in the following paragraphs. To investigate the mechanisms related to GBM cell functioning, there is a need for GBM models, which faithfully recapitulate the myriad complexities of GBM, including the diversity of cell functioning states.

I.3 GBM microenvironment

GBM cell line cultures in serum (classic established lines) have, for long, been extensively used. However, over the last decades, their use has become increasingly controversial because of important differences between the established lines and the original tumor [START_REF] Lee | Tumor stem cells derived from glioblastomas cultured in bFGF and EGF more closely mirror the phenotype and genotype of primary tumors than do serum-cultured cell lines[END_REF] [START_REF] Gillet | Redefining the relevance of established cancer cell lines to the study of mechanisms of clinical anti-cancer drug resistance[END_REF] [START_REF] Allen | Origin of the U87MG glioma cell line: Good news and bad news[END_REF]. Lee and colleagues extensively compared traditional serum-containing culture versus serum-free culture using cells derived from the same primary GBM tissues [START_REF] Lee | Tumor stem cells derived from glioblastomas cultured in bFGF and EGF more closely mirror the phenotype and genotype of primary tumors than do serum-cultured cell lines[END_REF]. They reported that genetic and epigenetic profiles of cells maintained in serum-containing medium dramatically changed over time compared to the original patient tumors. More importantly, they showed that cells maintained in serum-free medium, when xenografted into immunodeficient mice, yielded tumors with similar histology and expression profiles as their parental tumors, which is not the case with cells maintained in serum-containing medium. Cell culture models have thus evolved towards culture in defined serum-free media, which better preserves the genetic anomalies and histological characteristics of the original patient tumors, thus providing a more physiologically relevant model to study GBM [START_REF] Lee | Tumor stem cells derived from glioblastomas cultured in bFGF and EGF more closely mirror the phenotype and genotype of primary tumors than do serum-cultured cell lines[END_REF] [START_REF] Xie | The Human Glioblastoma Cell Culture Resource: Validated Cell Models Representing All Molecular Subtypes[END_REF].

While in vitro models have allowed important discoveries, they have their limits. Cell cultures are relatively homogeneous and do not recapitulate the unique and complex tumor architecture and tumor microenvironment (TME) (immune cells, vasculature, …) found in patients. Several efforts have been undertaken to develop 3D in vitro models mimicking part of these aspects. These 3D models are often designated as organoids, whether they correspond to simple cellular spheres obtained from cultures of patient-derived cells (PDC) (see for example [START_REF] Hubert | A Three-Dimensional Organoid Culture System Derived from Human Glioblastomas Recapitulates the Hypoxic Gradients and Cancer Stem Cell Heterogeneity of Tumors Found In Vivo[END_REF] or to complex tissues such as fragments of the tumor obtained after surgical resection (see for example (Jacob et al 2020)). Giant spheres obtained from culturing PDC in serum-free media under the form of cellular spheres embedded or not in an extracellular matrix recapitulate nutrient and oxygen gradients and some original histological features. More complex cell heterogeneity has been achieved with dissociated patient tumor cells bioprinted with brain ECM and supporting cells (Bioprinted GBM organoids), with cerebral organoid genetically engineered to develop GBM-like tumors (Neoplastic Cerebral Organoid or neoCOR) and with co-culturing cerebral organoids with dissociated patient tumor cells (Glioblastoma Cerebral Organoid or GLICO). The last organoid model developed correspond to cultures of minced pieces of resected patient tumor (GBM organoids or GBO). Organoid-GBM models have been successfully used to investigate TME importance in maintaining the cellular diversity observed in GBM [START_REF] Pine | Tumor Microenvironment Is Critical for the Maintenance of Cellular States Found in Primary Glioblastomas[END_REF], to study GBM invasion [START_REF] Krieger | Modeling glioblastoma invasion using human brain organoids and single-cell transcriptomics[END_REF] and to evaluate drug response [START_REF] Loong | Patient-derived tumor organoid predicts drugs response in glioblastoma: A step forward in personalized cancer therapy?[END_REF]. While each model has their strengths and weaknesses, GBOs have several strengths, compared to the other models, for studying functional heterogeneity in GBM. They remarkably resemble the patient samples from which they are derived -several features are recapitulated, including histological features, nutrient and oxygen gradients, cellular heterogeneity, transcriptional profiles, invasiveness in vivo and responses to treatment [START_REF] Hubert | A Three-Dimensional Organoid Culture System Derived from Human Glioblastomas Recapitulates the Hypoxic Gradients and Cancer Stem Cell Heterogeneity of Tumors Found In Vivo[END_REF] (Jacob et al 2020). They may thus be good models to test drug responses as well. They can be rapidly established. They come directly from patient tumors, are not limited to defined mutations or oncogenes and thus, provide the closest view of the original tumors [START_REF] Rybin | Organoid Models of Glioblastoma and Their Role in Drug Discovery[END_REF]. Nonetheless, GBOs lose over time the neurons they contain and lack a functional vasculature.

One way to overcome these limits is to use genetically engineered mouse models (GEMMs) based on main driver mutations found in primary GBM. GEMMs are useful for studying mechanisms linked to specific driver mutations and role of TME [START_REF] Hicks | Contemporary Mouse Models in Glioma Research[END_REF]. However, because GEMMs are usually composed of cells with homogeneous genetic background, they may not reflect the full intra-tumoral genomic and phenotypic heterogeneity observed in GBM. Patient-derived xenografts (PDX) have the advantage of preserving both the genetic and histological features of the primary tumor from which it was derived like invasiveness [START_REF] Wang | A reproducible brain tumour model established from human glioblastoma biopsies[END_REF] [START_REF] Kijima | Mouse models of glioblastoma[END_REF]. However, PDX models also have some disadvantages. They cannot be established from all patients, it usually takes long for tumors to develop (2-11 months) and they require immunodeficient mice (to prevent rejection of xenografted human cells), hence not recapitulating the tumor immune landscape as in patients. A recent analysis of whole-genome and RNA sequencing of over one hundred PDC with tumor initiating properties, and matched xenografts and parental tumors was recently published [START_REF] Shen | Comprehensive genomic profiling of glioblastoma tumors, BTICs, and xenografts reveals stability and adaptation to growth environments[END_REF]. The results show that PDC and xenografts are akin to their parental tumor at the genomic level but can differ at the mRNA expression and epigenomic levels, most likely because of the different growth environment for each sample type. One example of these bi-directional interactions between malignant and non-malignant cells is the dynamic cross-talk between GBM cells and endothelial cells (EC [START_REF] Lathia | Direct in vivo evidence for tumor propagation by glioblastoma cancer stem cells[END_REF] and that anti-angiogenic drugs eliminate most stem-like GBM cells from the tumor mass [START_REF] Calabrese | A perivascular niche for brain tumor stem cells[END_REF]. GBM cells also co-opt local vasculature for migration (Scherer 1938) (Watkins et al 2014) (Cuddapah et al 2014). In turn, GBM cells can regulate the vasculature and extent of angiogenesis in the tumor through several mechanisms. They can secrete cytokines and chemokines that activate EC [START_REF] Bao | Stem cell-like glioma cells promote tumor angiogenesis through vascular endothelial growth factor[END_REF]) [START_REF] Folkins | Glioma tumor stem-like cells promote tumor angiogenesis and vasculogenesis via vascular endothelial growth factor and stromal-derived factor 1[END_REF]) [START_REF] Thirant | Differential proteomic analysis of human glioblastoma and neural stem cells reveals HDGF as a novel angiogenic secreted factor[END_REF] or communicate through gap junctions [START_REF] Thuringer | Transfer of functional microRNAs between glioblastoma and microvascular endothelial cells through gap junctions[END_REF]. They can also transdifferentiate into EC or pericytes (Wang et al 2010) (Soda et al 2011) [START_REF] Cheng | Glioblastoma stem cells generate vascular pericytes to support vessel function and tumor growth[END_REF]. Moreover, GBM cells have vasculogenic mimicry capacities, i.e. they can form fluid-conducing channels [START_REF] Yao | Vascular endothelial growth factor receptor 2 (VEGFR-2) plays a key role in vasculogenic mimicry formation, neovascularization and tumor initiation by Glioma stem-like cells[END_REF].

Organoid-GBM models

GBM TME also contains numerous immune cells (brain-resident microglia, infiltrating monocytes that differentiate into macrophages, mast cells, neutrophils and lymphocytes), with monocyte-derived macrophages and lymphocytes being the most and the least abundant immune cell types in IDH wt -GBM respectively [START_REF] Klemm | Interrogation of the Microenvironmental Landscape in Brain Tumors Reveals Disease-Specific Alterations of Immune Cells[END_REF]. GBM cells secrete cytokines and chemokines that induce recruitment and/or stimulation of immune cells in the tumor to help sustain GBM progression. For instance, soluble factors released from GBM cells trigger MMP9 expression in microglia or macrophages [START_REF] Hu | Glioma-associated microglial MMP9 expression is upregulated by TLR2 signaling and sensitive to minocycline[END_REF], which can subsequently influence GBM growth and invasiveness. Other examples of GBM-secreted factors acting on immune cells include GDNF [START_REF] Ku | GDNF mediates glioblastomainduced microglia attraction but not astrogliosis[END_REF], CCL2 [START_REF] Zhang | A dialog between glioma and microglia that promotes tumor invasiveness through the CCL2/CCR2/interleukin-6 axis[END_REF], Versican (Hu et al 2015a), periostin (POSTN) [START_REF] Zhou | Periostin secreted by glioblastoma stem cells recruits M2 tumour-associated macrophages and promotes malignant growth[END_REF] and CSF1 [START_REF] Nijaguna | Glioblastoma-derived Macrophage Colony-stimulating Factor (MCSF) Induces Microglial Release of Insulin-like Growth Factorbinding Protein 1 (IGFBP1) to Promote Angiogenesis[END_REF]. Furthermore, GBM-secreted factors (e.g. TGF-β, PGE2, FasL) can act on surrounding immune cells to establish an immunosuppressive microenvironment [START_REF] Li | The molecular profile of microglia under the influence of glioma[END_REF]. Immune evasion can also be enhanced through inhibition of T-cell activity in GBM via PD-L1

found on glioblastoma extracellular vesicles [START_REF] Ricklefs | Immune evasion mediated by PD-L1 on glioblastoma-derived extracellular vesicles[END_REF]. Besides, recruitment of immune cells into the tumor helps per se to establish an immunosuppressive microenvironment. For instance, microglia in GBM release IL-10, a cytokine with a broad range of immunosuppressive roles [START_REF] Wagner | Microglial/macrophage expression of interleukin 10 in human glioblastomas[END_REF]. In addition, the capacity of these cells to present antigens to helper and cytotoxic T cells and their tumor-phagocytic function are impaired -they seem to rather adopt a tumor-supportive role [START_REF] Li | The molecular profile of microglia under the influence of glioma[END_REF]. Moreover, recruited immune cells can secrete factors that modulate GBM cell functioning. Microglial stimulation by GBM-released CCL2 induces IL-6 release from microglia, which then enhances GBM cell motility [START_REF] Zhang | A dialog between glioma and microglia that promotes tumor invasiveness through the CCL2/CCR2/interleukin-6 axis[END_REF]. GBM cell motility is also stimulated by microglia/macrophage-released TGF-β1 [START_REF] Wesolowska | Microglia-derived TGF-beta as an important regulator of glioblastoma invasion--an inhibition of TGF-betadependent effects by shRNA against human TGF-beta type II receptor[END_REF]) [START_REF] Ye | Tumor-associated microglia/macrophages enhance the invasion of glioma stem-like cells via TGF-beta1 signaling pathway[END_REF]. Interactions between GBM cells and immune cells are also instrumental in the remodeling of the extracellular matrix. GBM-released versican triggers MT-MMP expression in microglia via activation of the TLR2 receptor (Hu et al 2015a). This transmembrane metalloprotease (also called MMP14) cleaves the precursor form of MMP2 released by GBM cells into its active form, thereby promoting GBM cell invasion and growth [START_REF] Broekman | Multidimensional communication in the microenvirons of glioblastoma[END_REF]. In addition, CSF1 secreted by GBM cells induces IGFBP1 release from microglia, which can induce angiogenesis [START_REF] Nijaguna | Glioblastoma-derived Macrophage Colony-stimulating Factor (MCSF) Induces Microglial Release of Insulin-like Growth Factorbinding Protein 1 (IGFBP1) to Promote Angiogenesis[END_REF]. Cross-talk between GBM cells and immune cells is further supported by recent single-cell studies. Darmanis and colleagues reported that tumor-residing myeloid cells may have anti-inflammatory and proangiogenic roles in the tumor core, thus being endowed with the ability to enhance GBM growth, survival, and dissemination (Darmanis et al 2017). Caruso and colleagues inferred interactions between GBM cells and immune cells from single-cell RNA sequencing (scRNA-seq) data [START_REF] Caruso | A map of tumor-host interactions in glioma at single-cell resolution[END_REF]. They obtained evidences that these interactions regulate ECM remodeling, cell chemotaxis There is also an important interplay between GBM cells and surrounding astrocytes. It has been shown that the conditioned media of GBM cells inhibits expression of the tumor suppressor gene p53 in normal astrocytes, which subsequently changes ECM composition and promotes GBM cell malignancy [START_REF] Biasoli | Glioblastoma cells inhibit astrocytic p53-expression favoring cancer malignancy[END_REF]. Extracellular vesicles released from GBM cells alter cytokine and chemokine production and intracellular signaling in astrocytes, generating a growth-stimulating environment for GBM cells [START_REF] Oushy | Glioblastoma multiformederived extracellular vesicles drive normal astrocytes towards a tumour-enhancing phenotype[END_REF]. Moreover, surrounding astrocytes communicating with GBM cells overexpress numerous factors and activate signaling pathways known to regulate proliferation, motility and angiogenesis in GBM [START_REF] Zhang | Novel insights into astrocyte-mediated signaling of proliferation, invasion and tumor immune microenvironment in glioblastoma[END_REF]. For example, communication between astrocytes and GBM cells through gap junctions has been shown to be crucial for GBM cell functioning. Astrocytic connexin43 (Cx43), a gap junction protein, was required for the formation of an invasive niche in a mouse GBM model (Sin et al 2016). Modulation of GBM cell invasiveness appears also as a result from microRNA transfers between GBM cells and astrocytes via gap junctions [START_REF] Hong | Gap junctions modulate glioma invasion by direct transfer of microRNA[END_REF]. The GBM-astrocyte gap junctions also protect against cytotoxic effects of chemotherapy and enhance chemoresistance (Gielen et al 2013) (Chen et al 2015) [START_REF] Lin | Astrocytes protect glioma cells from chemotherapy and upregulate survival genes via gap junctional communication[END_REF].

Besides, astrocytes can enhance radioresistance in GBM by producing cytokines and chemokines that change gene expression profiles of GBM cells with stem-like properties, notably through STAT3 mobilization [START_REF] Rath | Coculture with astrocytes reduces the radiosensitivity of glioblastoma stem-like cells and identifies additional targets for radiosensitization[END_REF]. Moreover, astrocytes may provide nutrients to GBM cells so as to meet their high metabolic needs. For instance, surrounding astrocytes supply glutamine to malignant cells to support GBM growth [START_REF] Tardito | Glutamine synthetase activity fuels nucleotide biosynthesis and supports growth of glutamine-restricted glioblastoma[END_REF]. Astrocytes can also communicate with other cells in the TME so as to sustain GBM cell growth. For example, there exists a cross-talk between GBM cells, astrocytes and immune cells, whereby astrocyte-microglia interaction within GBM TME mediates specific transcriptional reprogramming in microglia and myeloid cells towards an anti-inflammatory, immunosuppressive state [START_REF] Heiland | Tumor-associated reactive astrocytes aid the evolution of immunosuppressive environment in glioblastoma[END_REF].

Oligodendrocyte progenitor cells (OPCs) present in the TME have been shown to secrete factors that enhance GBM cell stem-like properties and also trigger acquisition of chemo-resistant properties [START_REF] Hide | Oligodendrocyte Progenitor Cells and Macrophages/Microglia Produce Glioma Stem Cell Niches at the Tumor Border[END_REF]. Interestingly, increased OPC number has been observed at the tumor border, together with microglia/macrophages, and in vitro evidences support a pro-survival effect of GBM cell and macrophage secretomes on OPCs [START_REF] Hide | Oligodendrocyte Progenitor Cells and Macrophages/Microglia Produce Glioma Stem Cell Niches at the Tumor Border[END_REF]. In addition, oligodendrocytes may promote GBM cell invasiveness via angiopoietin-2 signaling pathway [START_REF] Kawashima | Oligodendrocytes Up-regulate the Invasive Activity of Glioblastoma Cells via the Angiopoietin-2 Signaling Pathway[END_REF].

Communication between GBM cells and oligodendrocytes is further supported by a study inferring interactions between GBM cells and oligodendrocytes from scRNA-seq data [START_REF] Caruso | A map of tumor-host interactions in glioma at single-cell resolution[END_REF]. The authors found that interactions between GBM cells and oligodendrocytes may regulate cell growth, cell-cell adhesion, angiogenesis and tumorigenesis.

Cross-talk between neurons and GBM cells has also been observed. Neurons favor GBM proliferation and growth through a neuronal activity-dependent secretion of the synaptic protein neuroligin 3 (NLGN3) (Venkatesh et al 2015) [START_REF] Venkatesh | Targeting neuronal activityregulated neuroligin-3 dependency in high-grade glioma[END_REF]. NLGN3 is a postsynaptic adhesion molecule secreted by neurons and OPCs [START_REF] Broekman | Multidimensional communication in the microenvirons of glioblastoma[END_REF]. Neuronal activity also leads to the secretion of neurotrophins like BDNF [START_REF] Venkatesh | Neuronal Activity Promotes Glioma Growth through Neuroligin-3 Secretion[END_REF], which play a role in GBM proliferation and growth (Venkatesh et al 2015) [START_REF] Lawn | Neurotrophin signaling via TrkB and TrkC receptors promotes the growth of brain tumor-initiating cells[END_REF]. Besides, neurons may communicate with GBM cells through neurotransmitters. For example, modulation of dopaminergic, serotonergic and cholinergic pathways has been shown to affect GBM cell functioning, with dopamine signaling being important for their proliferation and survival [START_REF] Dolma | Inhibition of Dopamine Receptor D4 Impedes Autophagic Flux, Proliferation, and Survival of Glioblastoma Stem Cells[END_REF]. Dopamine is a catecholamine primarily synthesized by neurons [START_REF] Dolma | Inhibition of Dopamine Receptor D4 Impedes Autophagic Flux, Proliferation, and Survival of Glioblastoma Stem Cells[END_REF]. Glutamate is another neurotransmitter that may be involved in neural regulation of GBM cell functioning since glutamate signaling has been shown to promote GBM cell proliferation and migration [START_REF] Ishiuchi | Blockage of Ca(2+)-permeable AMPA receptors suppresses migration and induces apoptosis in human glioblastoma cells[END_REF]) [START_REF] Ishiuchi | Ca2+-permeable AMPA receptors regulate growth of human glioblastoma via Akt activation[END_REF]. However, whether neuronal glutamate release contributes to GBM progression or whether it would signal in the same manner as GBM-derived glutamate has not yet been explored.

Moreover, ligands present on neurons may also support GBM cell functioning. For instance, a NOTCH ligand expressed along axons, Jagged, has been shown to support GSC invasion along unmyelinated white matter tracts (Wang et al 2019a). GBM cells can also regulate the functioning of neurons.

Glioma-secreted glutamate leads to hyper-excitability of neural circuits [START_REF] Campbell | Human glioma cells induce hyperexcitability in cortical networks[END_REF], and thus may increase release of activity-dependent mitogens by other cells in the TME (e.g. NLGN3 from neurons and OPC and neurotrophins from neurons). Finally, functional glutamatergic synapses formed by neurons on microtubes extending from GBM cells have been identified in patients' GBM [START_REF] Venkatesh | Electrical and synaptic integration of glioma into neural circuits[END_REF]) [START_REF] Venkataramani | Glutamatergic synaptic input to glioma cells drives brain tumour progression[END_REF].

ECM, the non-cellular component of the TME, functions as a source of biochemical signals [START_REF] Zhang | Novel insights into astrocyte-mediated signaling of proliferation, invasion and tumor immune microenvironment in glioblastoma[END_REF]. In GBM, ECM dynamics are deregulated so as to remodel ECM and favor invasiveness of GBM (see § I.4). ECM may also support maintenance of cell functioning states in GBM (e.g. TNC for stemness [START_REF] Sarkar | Activation of NOTCH Signaling by Tenascin-C Promotes Growth of Human Brain Tumor-Initiating Cells[END_REF], integrin signaling for growth and invasion [START_REF] Martin | Caveolin-1 regulates glioblastoma aggressiveness through the control of alpha(5)beta(1) integrin expression and modulates glioblastoma responsiveness to SJ749, an alpha(5)beta(1) integrin antagonist[END_REF]) [START_REF] Haas | Integrin alpha7 Is a Functional Marker and Potential Therapeutic Target in Glioblastoma[END_REF]).

In summary, GBM cells seem to hijack almost all cells in the TME: they can co-opt existing blood vessels and stimulate angiogenesis, can disarm immune cells and favor an immunosuppressive environment, can manipulate all neural-lineage cells into supporting tumor progression, and can remodel the ECM. Given the multiple modes of communication between malignant and nonmalignant cells, the bidirectionality of the cross-talks and the cell plasticity in GBM, a complex and unique TME is created that allows GBM cells to adapt to microenvironmental changes, whether these changes result from the natural tumor growth or from therapeutic assaults. Cell motility encompasses two processes: migration and invasion. Cell migration is a cell's ability "to move around in a space that is freely available" (de Gooijer et al 2018). In contrast, cell invasion is a more active process corresponding to a cell's ability "to remodel its microenvironment, to penetrate and move into its surroundings" (de Gooijer et al 2018). Motile cells have to develop protrusions in the direction of migration, disrupt adhesion sites at the cell rear, degrade and remodel the extracellular matrix (ECM) and ultimately make new connections at the cell front by remodeling their cytoskeleton so that they are pushed forward (Armento et al 2017) (de Gooijer et al 2018). Cell motility is thus a highly integrated process which involves changes in the motile cell itself as well as in its microenvironment (Figure 6). Regardless of their genetic heterogeneity, motile GBM cells need to undergo these changes and adapt to the unique brain environment as well as to changing environments during the migratory/invasive process. Three major migration/invasion routes have been described. GBM cells invade the narrow and tortuous extracellular space of the brain parenchyma by actively changing the ECM composition (Scherer 1938). They also use the perivascular spaces surrounding blood vessels as migratory routes (Scherer 1938) (Montana & Sontheimer 2011) (Watkins et al 2014) (Cuddapah et al 2014). Finally, GBM cells can use white matter tracts, notably the corpus callosum, as a highway to invade the opposite brain hemisphere (Scherer 1938) (Pedersen et al 1995). GBM cells can invade either as individual cells (single cell migration, when cell-cell junctions are absent) or multicellular groups (collective migration, when cell-cell adhesions are retained) (Vollmann-Zwerenz et al 2020) (Volovetz et al 2020). Collective GBM cell migration is facilitated by the formation of an interconnected network of tumor microtubes (Osswald et al 2015).

I.4 Infiltrative nature of GBM cells

Cell motility is initiated by chemo-attractants or pro-migratory factors binding to cell surface receptors, which activate or inactivate small GTPases and result in cytoskeleton reorganization and formation of protrusions (filopodia, lamellipodia and invadopodia/podosomes) (Armento et al 2017).

Several small GTPases associated with cytoskeletal remodeling have been demonstrated to be required for GBM cell motility. Increased expression or activity of Rac1 [START_REF] Chuang | Role of synaptojanin 2 in glioma cell migration and invasion[END_REF] uPAR [START_REF] Mohanam | Biological significance of the expression of urokinasetype plasminogen activator receptors (uPARs) in brain tumors[END_REF]) [START_REF] Gilder | The Urokinase Receptor Induces a Mesenchymal Gene Expression Signature in Glioblastoma Cells and Promotes Tumor Cell Survival in Neurospheres[END_REF]. Besides, the interaction of hyaluronan with its receptor CD44 is required for GBM cell invasion [START_REF] Park | Role of hyaluronan in glioma invasion[END_REF] as well as secretion of invasionpromoting ECM glycoproteins (e.g. SPARC [START_REF] Rempel | SPARC modulates cell growth, attachment and migration of U87 glioma cells on brain extracellular matrix proteins[END_REF]) [START_REF] Kunigal | SPARC-induced migration of glioblastoma cell lines via uPA-uPAR signaling and activation of small GTPase RhoA[END_REF] (Lee et al 2013) [START_REF] Gerarduzzi | The Matrix Revolution: Matricellular Proteins and Restructuring of the Cancer Microenvironment[END_REF], TNC (Xia et al 2016), THBS1 (Daubon et al 2019)). ECM stiffness is a major regulator of GBM cell motility. Studies have shown that a rigid ECM, as in the perivascular space, promotes GBM cell migration and that invasive GBM produces stiffness-promoting factors like collagen, fibronectin (FN), and laminin [START_REF] Ulrich | The mechanical rigidity of the extracellular matrix regulates the structure, motility, and proliferation of glioma cells[END_REF]. Changes in oxygen and nutrient levels also affect cell migration. GBM cells activate or produce pro-migratory and pro-invasive factors in response to hypoxic and nutrient-starvation conditions in order to reach a more favorable environment (See (Armento et al 2017) for all references).

I.5 GBM cell heterogeneity and plasticity

In addition to the complex and unique TME in GBM and to their highly infiltrative nature, the extensive heterogeneity observed at different levels in GBM is another key GBM hallmark contributing to treatment failure. GBM tissues present a wide genomic alteration landscape (genetic In contrast to the cell's identity based on the cell lineage (e.g. neuronal, oligodendroglial, astroglial) and its level of developmental differentiation (e.g. stem cell, progenitor, mature cells), a cell functioning/functional state corresponds to the condition of a given cell which translates its physiological function in the whole tissue such as stem-like or not, proliferative or quiescent, motile or static, pro-angiogenic or not, tumorigenic (i.e. able to initiate tumors) or nontumorigenic, sensitive or resistant to therapies. In this manuscript, I chose to use the term functioning state to avoid confusions, the term phenotype being often used in the literature to designate cell lineages as well as functional states.

Heterogeneity is observed between GBM cells from distinct tumors (inter-tumor heterogeneity) but also within the same tumor (intra-tumor heterogeneity). Within the same tumor, cells bearing different genetic alterations co-exist. For example, distinct Receptor Tyrosine Kinase (RTK) alterations have been found in different cells from the same tumor. EGFR-amplified glioblastomas also contain PDGFRA-amplified tumor cells and different cells from the same tumor can express different EGFR or PDGFRA variants [START_REF] Snuderl | Mosaic amplification of multiple receptor tyrosine kinase genes in glioblastoma[END_REF]. Joint analysis of DNA-and RNAsequencing data from GBM tumors revealed that multiple EGFR transcript variants were present in a tumor and the same variant could be generated by different DNA alterations [START_REF] Francis | EGFR variant heterogeneity in glioblastoma resolved through single-nucleus sequencing[END_REF].

Besides, using single-nucleus sequencing of GBM tumors, Francis and colleagues showed that EGFR copy number varied greatly between individual cells and that different EGFR variants characterized different subclonal GBM cell populations [START_REF] Francis | EGFR variant heterogeneity in glioblastoma resolved through single-nucleus sequencing[END_REF]. These results suggest that in addition to multiple RTK alterations within one GBM [START_REF] Snuderl | Mosaic amplification of multiple receptor tyrosine kinase genes in glioblastoma[END_REF], clonal diversification may arise due to distinct variants in a single RTK [START_REF] Francis | EGFR variant heterogeneity in glioblastoma resolved through single-nucleus sequencing[END_REF] Three main models have been proposed to explain how heterogeneity arises in GBM (Figure 8) [START_REF] Inda | Glioblastoma multiforme: a look inside its heterogeneous nature[END_REF]. The first model to have been proposed is the clonal evolution model based on the Darwinian theory of evolution. According to this model, genetic or epigenetic alterations appear randomly, generating new phenotypes. These are maintained by natural selection pressures, that is, only the best fit clone(s) will survive and expand. This would result in the co-existence of different most fit clones within the same tumor as well as of weaker clones generated during tumor expansion [START_REF] Inda | Glioblastoma multiforme: a look inside its heterogeneous nature[END_REF]. This intra-tumoral heterogeneity may be enhanced by environmental changes. For instance, chemo-or radiotherapy induces environmental changes that may lead to the acquisition of a resistant phenotype [START_REF] Inda | Glioblastoma multiforme: a look inside its heterogeneous nature[END_REF]. A more recent perspective of this model views clonal evolution of tumors as a growing tree, whereby the trunk represents the original clone carrying driver mutations shared by all cells in the tumor, and the branches represent the other clones that emerged due to acquisition and accumulation of new mutations [START_REF] Inda | Glioblastoma multiforme: a look inside its heterogeneous nature[END_REF]. These branches reflect the genetic intra-tumoral heterogeneity. The second model is that of cancer stem cell (CSC).

This model is based on a hierarchical organization during tumorigenesis akin to that observed during Upon selective pressure, resistant CSCs will survive. In these two models, the surviving cells will then expand and cause tumor relapse. In the last model, heterogeneity occurs as a result of cellular plasticity in response to microenvironment cues. This model implies that cells adapting to selective pressure are the source of tumor relapse. [START_REF] Larjavaara | Incidence of gliomas by anatomic location[END_REF].

Data supporting each of these three models have been obtained from GBM studies.

Reconstructions of clonal evolution trees from genomic profiling using bulk and single-cell RNA-seq support a clonal evolution model (Sottoriva et [START_REF] Suva | The Glioma Stem Cell Model in the Era of Single-Cell Genomics[END_REF]. This model usually co-exists with other sources of heterogeneity even in cancers containing a hierarchy of tumorigenic and non-tumorigenic cells (See [START_REF] Meacham | Tumour heterogeneity and cancer cell plasticity[END_REF] for all references). In concordance with the cell plasticity model, several studies have shown that GBM cells display great genetic, epigenetic and functional plasticity in response to microenvironmental cues. Comparing the genetic profiles of GSCs from the tumor core to those from the tumor periphery, Piccirillo and colleagues found that even though they shared some common anomalies suggesting that they derived from common ancestor cells, they exhibited distinct genetic profiles indicating that the microenvironment altered their genetic profiles [START_REF] Piccirillo | Distinct pools of cancer stem-like cells coexist within human glioblastomas and display different tumorigenicity and independent genomic evolution[END_REF]. Combining FACS sorting with functional assays to evaluate self-renewal, proliferation and multipotency, Dirkse and colleagues demonstrated that GBM cell subpopulations underwent reversible functional state transitions over time and also in response to microenvironmental changes like hypoxia [START_REF] Dirkse | Stem cell-associated heterogeneity in Glioblastoma results from intrinsic tumor plasticity shaped by the microenvironment[END_REF]. Moreover, exposing GSCs to serum causes their conversion to a differentiated, non-tumorigenic state while withdrawing serum causes them to revert back to their stem-like tumorigenic state (Natsume et al 2013) [START_REF] Turchi | Tumorigenic potential of miR-18A* in glioma initiating cells requires NOTCH-1 signaling[END_REF]. Conversely, differentiated GBM cells can acquire stem-like and tumorigenic properties upon glucose deprivation (Flavahan et al 2013), low pH conditions [START_REF] Hjelmeland | Acidic stress promotes a glioma stem cell phenotype[END_REF], or low oxygen levels [START_REF] Heddleston | The hypoxic microenvironment maintains glioblastoma stem cells and promotes reprogramming towards a cancer stem cell phenotype[END_REF]. Chemotherapeutic agents (Auffinger et In the absence of oxygen (anaerobic condition), pyruvate is metabolized to lactate. In contrast, cancer cells seemed to preferentially metabolize glucose to lactate, despite oxygen availability (Figure 9). He termed this phenomenon as aerobic glycolysis also known as the "Warburg" effect [START_REF] Warburg | On the origin of cancer cells[END_REF]). Later, studies on proliferative tissues showed that this phenomenon was not restricted to cancer cells but was common to all proliferative cells [START_REF] Vander Heiden | Understanding the Warburg effect: the metabolic requirements of cell proliferation[END_REF].

Glycolysis is essential for GBM growth -it is significantly increased in GBM as in other cancers (Marin-Valencia et al 2012), and knockdown of several glycolytic genes (HK2, PFKP, ALDOA, PGAM1, ENO1, ENO2, and PDK1) impedes GBM growth [START_REF] Sanzey | Comprehensive analysis of glycolytic enzymes as therapeutic targets in the treatment of glioblastoma[END_REF].

The glycolytic shift occurs at the expense of mitochondrial OXPHOS, which is ∼9-fold more efficient for ATP production than aerobic glycolysis [START_REF] Vander Heiden | Understanding the Warburg effect: the metabolic requirements of cell proliferation[END_REF]. GBM cells compensate for the lower efficiency of ATP production of glycolysis in several ways. They consume more glucose than normal cells, as shown by the upregulation of glucose transporters [START_REF] Bache | mRNA expression levels of hypoxiainduced and stem cell-associated genes in human glioblastoma[END_REF] and by the increased glucose uptake (Marin-Valencia et al 2012) in GBM tissues compared with adjacent normal tissues. Moreover, they can mobilize additional energetic pathways. In GBM cells, glutamate produced during glutaminolysis may serve as a major substrate to refuel the TCA cycle (Oizel et al 2017). Glutamine metabolism in GBM cells is critical for their proliferation, survival and drug response [START_REF] Tanaka | Compensatory glutamine metabolism promotes glioblastoma resistance to mTOR inhibitor treatment[END_REF]. Fatty acids are also major bioenergetic substrates critical for GBM growth, with overexpression of several lipid carriers and enzymes critical for fatty acid oxidation (FAO), for example the CTP1 regulators of long-chain fatty acid transport into the mitochondria or the de novo fatty acid synthesis enzyme FASN (Lin et al 2017a). Carbohydrate, lipid and amino acid metabolisms all generate intermediates that converge to the TCA cycle to fuel energy production.

These alternative energy sources can also drive GBM growth and survival since in vivo, glucose provides only a minority of the acetyl-CoA pool in TCA cycle in GBM [START_REF] Maher | Metabolism of [U-13 C]glucose in human brain tumors in vivo[END_REF]. The energetic metabolic pathways mobilized in GBM cells may depend on their functioning state. For instance, supply cells with ribose 5-phosphate (R5P), which is a building block for nucleic acid synthesis (e.g. purine, pyrimidine). Glucose metabolism also appears instrumental in GBM cell therapeutic resistance. Escape from the anti-angiogenic Bevacizumab treatment has been linked to a further increase in glycolysis and its uncoupling from OXPHOS in favor of lactate production in in vivo GBM models as well as in patients [START_REF] Fack | Bevacizumab treatment induces metabolic adaptation toward anaerobic metabolism in glioblastomas[END_REF].

Hoang

Lipid metabolism, not only plays an important role in energy metabolism, but also provides raw materials for maintaining structural integrity of membranes and for cellular signaling (Garcia et al 2021). Citrate-derived acetyl CoA is used for lipid production. Isocitrate dehydrogenase 1 (IDH1), which catalyzes citrate formation from α-ketoglutarate in the cytoplasm, is upregulated in primary IDH wt -GBM, favoring tumor growth and therapy resistance (Wahl et al 2017) [START_REF] Calvert | Cancer-Associated IDH1 Promotes Growth and Resistance to Targeted Therapies in the Absence of Mutation[END_REF].

Amino acids such as glutamine also fuel bioenergetic reactions and macromolecule formation that support GBM cell functioning (See [START_REF] Natarajan | Glutamine Metabolism in Brain Tumors[END_REF] and (Garcia et al 2021) for reviews).

Metabolic variations in cancer cells have long been viewed as passive adaptations to cell

functioning. There is however evidence that changes in metabolic activities not only accompany changes in cell functioning but can also drive these changes. Glucose restriction was reported to promote acquisition of tumor-initiating properties by non-tumorigenic GBM cells through enhanced expression of the high affinity glucose transporter type 3 GLUT3 (SLC2A3) (Flavahan et al 2013).

Existence of variations in SLC2A3 expression within patients' tumors was however not sought for. Altogether the results acquired to this day support the idea that changes in metabolic activities are not only essential to sustain GBM cell functioning states, but can also regulate them.

They also highlight the heterogeneity of the metabolic activities at play within GBM cells, a heterogeneity likely to be associated with differing cell behaviors and/or genetic alterations.

Metabolism is downstream of all signaling pathways regulating cell functioning states irrespective of genomic specificities of tumors. Therefore, metabolic enzymes required for adopting and/or maintaining cell functional states crucial for tumor development constitute therapeutic targets able to overcome tumor heterogeneity. There are still only few studies focusing on GBM cell functioning states and their metabolic characteristics. A systemic view of the metabolic pathways accompanying and/or driving GBM cell functioning states is lacking. it is also possible to perform multi-modal sequencing (Table 2). Several assays for simultaneous measurement of RNA and proteins in single cells have been published over the last six years. In 2016,

II: SINGLE

Frei and colleagues developed a proximity ligation assay for RNA (PLAYR) and combined it to CyTOF for measuring transcripts and their corresponding proteins in the same single cells [START_REF] Frei | Highly multiplexed simultaneous detection of RNAs and proteins in single cells[END_REF].

In 2017, two higher-throughput assays combining scRNA-seq and measurement of proteins using published. Lareau and colleagues published the mtscATAC-seq method which allows simultaneous profiling of single-cell mitochondrial DNA (mtDNA) mutations and chromatin accessibility (Lareau et al 2021). An adaptation of this method gave rise to ATAC with Select Antigen Profiling by sequencing or ASAP-seq [START_REF] Mimitou | Scalable, multimodal profiling of chromatin accessibility, gene expression and protein levels in single cells[END_REF]. ASAP-seq allows to simultaneously profile three distinct modalities in single cells: accessible chromatin using scATAC-seq, cell surface and intracellular protein levels using oligo-labeled antibodies and optionally mitochondrial DNA using mtDNA genotyping. DOGMAseq is another single-cell triple omics technique, assessing chromatin accessibility, gene expression and protein [START_REF] Mimitou | Scalable, multimodal profiling of chromatin accessibility, gene expression and protein levels in single cells[END_REF]. It couples the protein barcode detection method used in CITE-seq to the 'Multiome' product from 10X Genomics which captures the chromatin accessibility and transcriptomic profiles from the same cells. Another technique capturing three distinct modalities from the same individual cells, a lower-throughput one in terms of number of cells analyzed, is the single-cell triple omics sequencing or scTrio-seq which simultaneously yields genomic (CNVs), DNA methylomic, and transcriptomic information from the same cell [START_REF] Hou | Single-cell triple omics sequencing reveals genetic, epigenetic, and transcriptomic heterogeneity in hepatocellular carcinomas[END_REF]. Despite all the technological advances, scRNA-seq remains the most accessible and widely used high-throughput approach. This section aims at presenting scRNA-seq techniques, the general analytical workflow and at discussing some related concepts.

II.1 Single-cell RNA sequencing techniques

Ensemble-based approaches (microarrays and next-generation sequencing (NGS) of highthroughput RNA sequencing (RNA-seq)) have long been used to profile gene expression from bulk cell populations. The resulting expression value for a given gene corresponds to an average of its expression levels across millions of cells (Figure 10). When analyzing tumor fragments, it corresponds to an average of its expression levels across millions of malignant as well as non-malignant cells. To better understand biological phenomena like cellular heterogeneity in cancer, it is crucial to access and analyze gene expression profiles of individual cells.

Measuring gene expression at the single-cell level started with low-throughput approaches such as reporter constructs or microscopy-coupled immunohistochemistry (IHC) at the protein level, and single-cell quantitative PCR (qPCR) or single-molecule RNA fluorescence in situ hybridization (RNA FISH) at the RNA level [START_REF] Stegle | Computational and analytical challenges in single-cell transcriptomics[END_REF]. These approaches provided important insights into transcriptional and translational kinetics and cell differentiation potential. However, given that only a single gene or a small number of genes could be studied, these approaches reflected only a very small part of the global active transcriptional landscape within a cell. Recent technological advances in the microfluidics and nanotechnological fields have allowed to have access to the wider picture, with high-throughput profiling in each cell, with large numbers of cells simultaneously profiled and with dramatically reduced costs [START_REF] Tirosh | Dissecting human gliomas by single-cell RNA sequencing[END_REF]) [START_REF] Li | Dropout imputation and batch effect correction for single-cell RNA sequencing data[END_REF]. Single-cell sequencing techniques have been now extended to high-throughput DNA and epigenome profiling, and even multi-omics profilings [START_REF] Mossner | The role of single-cell sequencing in studying tumour evolution[END_REF]. It is important to remember that the obtained data correspond only to a temporal and spatial snapshot of the sample, and thus may not reflect the whole intra-sample heterogeneity and the dynamics of gene expression and cell functioning states.

Over the last decades, scRNA-seq has gained popularity because of the relative ease of profiling RNA, and also because it is less affected by the low amount of starting material within a single cell and the partial capture rate of experimental protocols compared to other single-cell genomics profiling methods [START_REF] Tirosh | Dissecting human gliomas by single-cell RNA sequencing[END_REF]. Generating scRNA-seq data from a biological sample requires multiple steps (Figure 11). The first step is isolation of viable single cells from the tissue of interest, which has to be dissociated into individual cells prior single-cell isolation, except if using laser capture microdissection for cell isolation. Single-cell isolation can be achieved using with a charged couple device by exciting the fluorophores with appropriate lasers [START_REF] Heather | The sequence of sequencers: The history of sequencing DNA[END_REF].

Several scRNA-seq technologies have been developed (Table 2), which differ in two key aspects: single-cell isolation approach (plate-based or microdroplet-based) and transcript coverage.

Plate-based methods usually use FACS or Fluidigm C1 to isolate and place single cells into individual wells on a single plate that contain all chemicals necessary for cell lysis, mRNA capture, RT and amplification. Commonly used plate-based methods include technologies such as SMART-seq [START_REF] Ramskold | Full-length mRNA-Seq from single-cell levels of RNA and individual circulating tumor cells[END_REF], SMART-seq2 [START_REF] Picelli | Full-length RNA-seq from single cells using Smart-seq2[END_REF] and MARS-Seq [START_REF] Jaitin | Massively parallel single-cell RNA-seq for marker-free decomposition of tissues into cell types[END_REF]. A main disadvantage of plate-based methods is the limited number of single cells that can be profiled at the same time (e.g. maximum of 800 cells for Fluidigm C1 system) [START_REF] Li | Dropout imputation and batch effect correction for single-cell RNA sequencing data[END_REF]. Microdroplet-based methods overcome this issue, by allowing simultaneous capture of tens of thousands of cells, with a reduced reagent cost. Cells and gel beads are loaded at controlled flow rates into channels of a microfluidic chip so as to generate droplets. Each droplet contains one single cell, reagents for cell lysis and RT, and a microbead with DNA sequences coupled to a cell barcode and a RNA barcode (UMI). Each transcript having a unique cell barcode, the remaining steps of the library preparation can thus be performed in bulk, hence reducing cost. Another advantage is that the use of UMI allows to distinguish between original molecules and amplification duplicates, thus improving quantification of mRNA molecules [START_REF] Ziegenhain | Comparative Analysis of Single-Cell RNA Sequencing Methods[END_REF]. Examples include Drop-seq [START_REF] Macosko | Highly Parallel Genome-wide Expression Profiling of Individual Cells Using Nanoliter Droplets[END_REF], inDrop [START_REF] Klein | Droplet barcoding for single-cell transcriptomics applied to embryonic stem cells[END_REF], and 10X Genomics Chromium (10X) [START_REF] Zheng | Massively parallel digital transcriptional profiling of single cells[END_REF]. Microdroplet-based methods gained popularity because of their high throughput, their fast and simple workflow, their consistent data quality and their relatively low cost. Recently, new scRNA-seq technologies have been developed. For example, Rosenberg and colleagues proposed a method called SPLiT-seq that uses combinatorial barcoding to profile single-cell transcriptomes without physically isolating each cell [START_REF] Rosenberg | Single-cell profiling of the developing mouse brain and spinal cord with split-pool barcoding[END_REF]. They can simultaneously sequence tens of thousands of cells from multiple biological samples in one single experiment, hence reducing batch effects that are introduced during library preparation and sequencing. The second main difference between sequencing technologies is in the coverage of transcript. Technologies like SMART-seq offer fulllength transcript coverage, which is typically accompanied by a higher transcript capture rate.

However, for accurate quantification of transcripts in full length, higher sequencing depth is required. In contrast, the more commonly used technologies such as MARS-seq, Drop-seq, 10X

Genomics and SPLiT-seq rather capture only the 3'-end of transcripts, thus requiring relatively shallow sequencing depth, and reducing sequencing cost [START_REF] Li | Dropout imputation and batch effect correction for single-cell RNA sequencing data[END_REF].

Ziegenhain and colleagues compared the performance of six scRNA-seq technologies using homogeneous cell populations as well as poly-adenylated synthetic RNA transcripts of known concentration [START_REF] Ziegenhain | Comparative Analysis of Single-Cell RNA Sequencing Methods[END_REF]. Their results showed substantial technical variations between methods in terms of sensitivity (i.e. the probability to capture and convert a mRNA molecule present in a single cell into a cDNA molecule), of accuracy (i.e. how well does the quantification of detected read counts reflect the actual mRNA levels), and of precision (i.e. the degree of variability in gene detection and quantification of detected read counts). SMART-seq2 was the most sensitive and accurate technology, and also showed the most even transcript coverage. Of note, all methods performed well in terms of accuracy. Precision was assessed based on two factors: how often a transcript is measured (dropout probability) and with how much variation it is measured. While SMART-seq2 performed better in terms of dropout rates, UMI-based technologies performed better in terms of variability of quantification. The latter quantified mRNA levels with less amplification noise. Altogether, these results show that each technology has its advantages and weakpoints.

These factors, together with number of cells analyzed, determine the capacity of a scRNA-seq method to detect relative differences in expression levels.

II.2 Data analysis

Development of analytical tools for scRNA-seq data remains the object of an intense area of research. While some user-friendly 'plug-and-play' tools or pipelines exist (e.g. Loupe from 10X Genomics), they may behave as a 'black box' [START_REF] Haque | A practical guide to single-cell RNA-sequencing for biomedical research and clinical applications[END_REF] since the exact algorithmic details and parameters employed are not always accessible and cannot always be modified to suit the study.

Analytical tools have to be adapted for the types of samples being studied (e.g. homogeneous or heterogeneous samples, malignant or non-malignant cells) and to take into account technical drawbacks ranging from poor quality samples to stochastic failures in transcript detection. Despite the current belief that big data analysis will by itself provide "spontaneously" unbiased information, it remains necessary to tailor the organization of the analytical workflow (Figure 12) using carefully chosen analytical tools according to the type of samples analyzed (e.g. cancer or normal tissues) and the information sought for (e.g. cell lineage identity, cell developmental stage, molecular substrates of cell-cell interactions, cell functional states, ...).

Once sequencing data are obtained, the quality of entire sequencing reads needs to be assessed by evaluating base-call errors and GC content for example, and identifying low-quality bases [START_REF] Hwang | Single-cell RNA sequencing technologies and bioinformatics pipelines[END_REF]. This quality control can be achieved using FastQC tool (BabrahamBioinformatics)

for instance. Adapter sequences should also be removed. Then, sequencing reads are aligned to a reference genome using tools such as HISAT2 (Kim et abundance can be quantified as described by Islam and colleagues [START_REF] Islam | Quantitative single-cell RNA-seq with unique molecular identifiers[END_REF]. The gene expression matrix -observations (i.e. cells) × features (i.e. genes) -is generated using only reads mapping to exonic loci with high quality [START_REF] Hwang | Single-cell RNA sequencing technologies and bioinformatics pipelines[END_REF]. Of note, most of the tools used for the abovementioned steps were originally designed for bulk RNA-seq data.

Raw read counts need to be pre-processed in order to make samples comparable prior downstream analyses. For scRNA-seq derived from cancer tissues, inference of chromosome copy number variations (CNV) from RNA expression profiles allows to distinguish cancer cells from normal cells. Then, the high-dimensional data are usually reduced to a low dimensional space, which can then be used for downstream analyses. Clustering is the most commonly used downstream analysis.

It allows to group cells with similar expression profiles, thereby providing a view of the cell heterogeneity. The use of trajectory modeling in downstream analyses is constantly increasing.

Initially adopted for identifying the differentiation paths taken by the cells in normal developing tissues, trajectory modeling is now extended to cancer tissues. Next, genes characterizing the cell clusters or the path embranchments of the cell trajectories can be identified using differential expression analysis. To understand the biological relevance of these genes, one can proceed with combinations of various approaches such as literature review, gene ontology and pathway analyses, gene expression network modeling or patient survival analysis in disease context. These gene-level analyses allow to identify the molecular characteristics associated with a given cell population.

II. 3 Data pre-processing

Most common pre-processing steps in scRNA-seq data analysis are presented in this paragraph. They can be used in a sequential pipeline, and adapted depending on downstream analyses.

a. Quality control (QC) is performed before analyzing scRNA-seq data. Cell transcriptome quality is evaluated based on three main criteria, though there is no consensus on the values of filtering parameters: (1) number of read counts detected per cell, (2) number of genes detected per cell, and (3) fraction of counts from mitochondrial genes detected per cell. Cells with low number of read counts and genes detected, and a high fraction of mitochondrial counts are likely to be cells whose cytoplasmic mRNA has leaked out through a broken membrane, thus retaining mostly mitochondrial mRNA. These cells are usually excluded from further analyses. Likewise, cells with very high number of read counts and genes detected may correspond to cell doublets, and are also excluded. These QC filters should be adapted to each dataset and should be set as permissive as possible to avoid filtering out viable, relevant cell populations. Data from heterogeneous tissues (e.g. cancers like GBM) may exhibit multi-modal distributions of total number of genes and read counts detected, reflecting the heterogeneous nature of the sample even at the level of QC metrics. Genes detected in too few cells should also be excluded to avoid potential analytical bias due to scarcely detected genes. While there is no consensus on the ideal number of cells in which a gene should be detected, the threshold should not be too high since it would prevent the identification of cell groups with very few cells. In our analyses, we used 3 as the threshold as we consider the minimum possible number of cells in a group to be 3.

b. Data normalization is a key step. Given that raw read counts (measured data) are directly influenced by technical variations like differences in sequencing depth, cell lysis and RT efficiency, normalization is required to account for these technical biases, prior to comparing gene expression between cells. A commonly used normalization is that based on feature variance, which standardizes data to z values such that their mean expression will equal to 0 and their standard deviation to 1

[z=(x-μ)/σ, where z is the standardized value, x the raw value, μ the mean expression value across all cells and σ its standard deviation]. For example, this step is part of the widely used analytical Seurat workflow [START_REF] Satija | Spatial reconstruction of single-cell gene expression data[END_REF]. While such normalization is mandatory when the units of variables are In that case, the number of gene-mapped reads is divided by the transcript's length so as to take into account the length of transcripts in the estimation of read counts. Gene-length normalization is required for within-sample comparisons (i.e. when comparing expression of different genes within a cell) [START_REF] Evans | Selecting between-sample RNA-Seq normalization methods from the perspective of their assumptions[END_REF]. Gene-length normalization may be beneficial for full-length sequencing methods which exhibit gene-length bias, with shorter genes having lower counts and higher dropout rate [START_REF] Phipson | Gene length and detection bias in single cell RNA sequencing protocols[END_REF]. In contrast, UMI-based methods do not exhibit gene length bias and have a mostly uniform dropout distribution across genes of varying length [START_REF] Phipson | Gene length and detection bias in single cell RNA sequencing protocols[END_REF]. Therefore, gene-length normalization of UMI counts is not appropriate since it will artificially inflate the expression values of shorter genes relative to longer genes [START_REF] Phipson | Gene length and detection bias in single cell RNA sequencing protocols[END_REF]. Once normalized, data are usually log(x+1)-transformed so as to reduce the skewness of their distribution.

c. Regressing out biological confounders can also be performed. Transcriptional bursts, cell cycle changes or cell size, for example, lead to variability in transcript levels, and thus may be considered as biological noise or confounders (Jia et al 2017) [START_REF] Eling | Challenges in measuring and understanding biological noise[END_REF]. In that case, they are regressed out, with the assumption that they might mask relevant biological signals. The most common biological confounder that is regressed out is cell cycle effects on the transcriptome. This can be done using simple linear regressions against a cell cycle score as implemented in the Seurat algorithm [START_REF] Satija | Spatial reconstruction of single-cell gene expression data[END_REF] or using more complex mixture models as implemented in the scLVM algorithm [START_REF] Buettner | Computational analysis of cell-to-cell heterogeneity in single-cell RNA-sequencing data reveals hidden subpopulations of cells[END_REF]. These methods can also be used to regress out other known biological effects such as mitochondrial gene expression, which may be indicative of cell stress [START_REF] Chen | Controlling for Confounding Effects in Single Cell RNA Sequencing Studies Using both Control and Target Genes[END_REF]). However, these biological factors may also be informative of the cells' biology (e.g. cell cycle status of proliferative cells in cancer (Tirosh et al 2016b)). Within an organism, several biological processes co-occur. Hence, there may be dependencies between different processes, and correcting for one process may mask or modify the signal of another. For these reasons, we didn't correct for cell cycle effects or other biological factors in our analyses.

d. Taking into account dropouts is linked to a single-cell data characteristic, the presence of "dropouts". In scRNA-seq data, low amounts of RNA molecules in individual cells, inefficient mRNA capture, stochastic mRNA expression and amplification failures result in transcripts being missed [START_REF] Qiu | Embracing the dropouts in single-cell RNA-seq analysis[END_REF]. Consequently, a gene that is not detected can either be a gene that is actually not expressed (true zero) or a gene whose transcripts have been missed (dropouts). These dropout events cause scRNA-seq data to be highly sparse, with only a small fraction of the whole cell transcriptome being captured [START_REF] Qiu | Embracing the dropouts in single-cell RNA-seq analysis[END_REF]. One way to address dropouts is to adopt a two-step strategy, i.e. selecting genes that are assumed to be less affected by dropouts such as the most variable ones, followed by dimension reduction, as implemented in Seurat [START_REF] Satija | Spatial reconstruction of single-cell gene expression data[END_REF]. However, subtle gene expression variations can also have biological roles, and focusing on highly variable genes would ignore these aspects. Another way of dealing with dropouts is dropout imputation. This involves calculating the probability that a zero value for a given gene is due to its transcripts having been missed (i.e.

determining if a zero value is a dropout) and then replacing the dropouts by appropriate values as done by CIDR (Lin et al 2017b). This approach is based on the assumption that cells exhibiting similar overall transcriptome profiles should express homogeneous levels of each gene, which is not always the case in real data. During the course of this thesis work, we tested the influence of dropouts on subsequent analysis results. We showed that artificially imputing values to dropouts doesn't change results of downstream analysis, and that non-imputed data yield biologically coherent results ((Saurty-Seerunghen et al 2019), part 1 of the Results section). Also, biological factors such as differing tissues, spatial locations or species can be regarded as a batch effect because of eventual differences in their cell contents or technical processing. Thus, multiple datasets cannot be simply pooled together and analyzed -appropriate data integration methods that take into account the abovementioned factors are required. Two main strategies are adopted for joint analysis of multiple scRNA-seq datasets. They are mainly based on the assumption that the same cell populations (e.g. cell lineage types) should be found in different datasets, even though the proportion of each cell population may vary. Merging strategies use common cell populations as anchors. In other words, they identify the most similar cell groups between the differing datasets, and pool them in a single group (Figure 13A). As a result, we obtain a novel dataset where all cells from the starting datasets are merged. These approaches are implemented by several algorithms, including mutual nearest neighbor or MNN [START_REF] Haghverdi | Batch effects in single-cell RNA-sequencing data are corrected by matching mutual nearest neighbors[END_REF], canonical correlation analysis or CCA [START_REF] Butler | Integrating single-cell transcriptomic data across different conditions, technologies, and species[END_REF]) [START_REF] Stuart | Comprehensive Integration of Single-Cell Data[END_REF] and Sub-Type

Anchor Correction for Alignment in Seurat or STACAS (Andreatta & Carmona 2021a). In the second strategy called projection strategy, cells from a query dataset are projected onto existing references

II.4 Dimensionality reduction

Technological advances over the last decades have allowed for thousands of simultaneous measurements collected for a single sample, with a spectacular increase in the number of samples collected. Data are gathered into vectors whose dimensions correspond to the number of measurements (features), thus yielding high-dimensional biological datasets. The problem is that as dimensionality increases, the volume of the space increases exponentially, causing the available data to become sparse. All objects in that space appear dissimilar to one another, the space becomes practically impossible to explore exhaustively, and data organization strategies (e.g. clustering)

become inefficient. This phenomenon is referred to as the "Curse of Dimensionality" (Figure 14).

In addition, in single-cell datasets, the low capture efficiency and low sequencing depth per cell make the collected gene measurements very noisy, reflected by low expression counts and an excessive number of zeros (dropouts). For efficient and effective downstream analysis of scRNA-seq data, it is important to simplify high-dimensional noisy data in such a way that only the most meaningful information is preserved. This is where dimensionality reduction (DR) comes to the rescue. DR compresses data into a low dimensional space that captures the dominant sources of Reduced data are usually used by exploratory methods that require data summaries such as clustering and for computationally complex downstream analysis tools such as trajectory inference.

In contrast, measured data (raw data or processed data that retain the structure of zeros) are usually used for statistical testing as in differential expression analysis.

II.5. Clustering

Broadly speaking, clustering aims at organizing data into groups whose members share some similarities (Figure 15). When applied on single-cell transcriptomic data, clustering organizes cells into groups with similar transcriptional profiles based on distance metrics such Euclidean, cosine and correlation-based distances. Cells are assigned to clusters by minimizing intra-cluster distances or finding dense regions in the reduced expression space [START_REF] Luecken | Current best practices in single-cell RNA-seq analysis: a tutorial[END_REF]. Two approaches can be applied to cluster and annotate cells: supervised and unsupervised, the latter being most similarity (compactness) to inter-cluster dissimilarity (separation) like the Silhouette index [START_REF] Rousseeuw | Silhouettes: A graphical aid to the interpretation and validation of cluster analysis[END_REF]). Cluster number is thus determined based on the structure of the data, rather than by inference. K-means clustering is then used to group cells into N clusters.

Community detection methods are graph-partitioning algorithms which rely on a graph representation of scRNA-seq data obtained using a K-Nearest Neighbour approach (KNN graph) [START_REF] Luecken | Current best practices in single-cell RNA-seq analysis: a tutorial[END_REF]. Cells are represented as nodes and each cell is connected to its K most similar cells. Densely-connected regions of the graph correspond to densely-sampled regions of the expression space. Communities/Clusters correspond to groups of cells with a higher number of links between them in the graph than expected from the total number of links. Given that in community detection methods, only neighboring cell pairs have to be considered as belonging to the same cluster, the search space for possible clusters is greatly reduced, hence making this approach faster than other clustering methods [START_REF] Luecken | Current best practices in single-cell RNA-seq analysis: a tutorial[END_REF]. The community detection method applied to scRNA-seq data is the Louvain algorithm implemented in PhenoGraph [START_REF] Levine | Data-Driven Phenotypic Dissection of AML Reveals Progenitor-like Cells that Correlate with Prognosis[END_REF], in Seurat [START_REF] Stuart | Comprehensive Integration of Single-Cell Data[END_REF] or in SCANPY [START_REF] Wolf | SCANPY: large-scale single-cell gene expression data analysis[END_REF].

All the above-described unsupervised methods rely on intrinsic properties of the samples, usually without a priori assumptions. Unsupervised without a priori strategies were among the first considered due to their simplicity and the fact that clustering is performed solely based on the observed samples without the need for labeled training data. They have been repetitively used with success to identify differing cell-type lineages from both cancer and normal scRNA-seq datasets (Chung et [START_REF] Rosenberg | Single-cell profiling of the developing mouse brain and spinal cord with split-pool barcoding[END_REF]. However, unsupervised without a priori clustering analysis appears as a poorly-suited method to identify biological characteristics involving low-, rather than high-, amplitude changes in gene expression. For example, cells characterized as being in a proliferative state are dispersed between differing clusters following unsupervised without a priori clustering, and their identification requires the use of a molecular classifier (Tirosh et al 2016a). We were confronted to this issue when seeking for cell functioning states from GBM scRNA-seq data. As described in the first part of the Results section, we overcame this issue by developing an unsupervised with a priori clustering approach based on the use of a classifier corresponding to a molecular signature with experimentally-validated elements (signature-based approach).

Once data are clustered, the next step is to assess the validity of the identified clusters and to annotate them with a meaningful label. For omics data, this corresponds to successful annotation of the represented biology. To this aim, genes characterizing each cluster are identified and analyzed to determine if they are involved in specific biological cell processes, pathways or identities, if they are co-expressed or co-regulated and if they may have some clinical relevance. These aspects will be discussed later.

II.6 Trajectory modeling

In contrast to clustering which regards omics data as being discrete, trajectory inference (TI, also called pseudo-time analysis) regards them as a snapshot of a dynamic cellular process, which can be reconstructed by finding paths that minimize transcriptional changes between neighboring cells [START_REF] Luecken | Current best practices in single-cell RNA-seq analysis: a tutorial[END_REF]. Cell ordering along these paths is described by a pseudo-time variable. By revealing gene expression dynamics, TI thus aims at capturing transitions between cell identities or states (Figure 16). It includes two steps: dimensionality reduction and trajectory modeling.

Over the last decades, the number of available TI methods has exploded -they usually differ in the complexity of the modeled paths or topologies. A recent study assessing the performance, scalability, robustness and usability of 45 TI methods reported that no method performs optimally for all types of trajectories, that they are rather quite complementary and that their optimal performances depend on the characteristics of the data [START_REF] Saelens | A comparison of single-cell trajectory inference methods[END_REF]. PAGA, implemented in SCANPY python package [START_REF] Wolf | PAGA: graph abstraction reconciles clustering with trajectory inference through a topology preserving map of single cells[END_REF], outperformed other methods when topology is unknown.

For simpler topologies (linear or bifurcation), Slingshot [START_REF] Street | Slingshot: cell lineage and pseudotime inference for single-cell transcriptomics[END_REF] outperformed other methods. Altogether, Saelens and colleagues showed that seven topologies can be inferred: cycle, linear, bifurcation, multifurcation, tree, connected and disconnected. They concluded that TI methods should be chosen based on the expected topology [START_REF] Saelens | A comparison of single-cell trajectory inference methods[END_REF].

Both the static and dynamic nature of the data can be reconciled in coarse-grained graph representations (e.g. as implemented in PAGA), with clusters represented as nodes, and trajectories between clusters as edges [START_REF] Luecken | Current best practices in single-cell RNA-seq analysis: a tutorial[END_REF].

Once a trajectory is established, gene expression dynamics along the trajectory are analyzed, since trajectory-associated genes may correspond to genes regulating the modeled process. To add directionality to the trajectory, RNA velocities can be overlaid on inferred trajectories. The RNA velocity for a given gene in a given cell at a given time point is the instantaneous rate of change of spliced mRNA abundance, which can be determined by modeling the spliced and unspliced mRNA abundances in a cell [START_REF] Manno | RNA velocity of single cells[END_REF]. It can be used to infer a 'direction of change' and thus, the future state of individual cells in the gene expression space. Overlaying it onto the trajectory thus adds directionality to the trajectory and can be used to determine whether a biological process is The significance of a difference is evaluated using p-values (P). A result is said to be significant if P ≤ α, where α is the type I error rate chosen to control false positives. Thus, by setting α at 0.05, we accept a 5% chance for a result to be a false positive. As the number of tests performed increases, the number of false positives also increases. To mitigate this problem, multiple-testing correction methods are applied such as Bonferroni or Benjamini-Hochberg correction. Bonferroni is a simple method that controls the probability of incorrectly rejecting the null hypothesis (H0) in a single test among all tests done, i.e. it controls the Family-Wise Error Rate (FWER) [START_REF] Krzywinski | Comparing samples-part II[END_REF]. If H0 for all tests are true, then the probability of false positives due to chance in all tests is expected to be 0.05. This method is very stringent and is applied when it is imperative not to have any false positives. A less stringent alternative is the Benjamini-Hochberg adjustment, which controls for False Discovery Rates (FDR). FDR is the expected proportion of false discoveries among all discoveries. It is important to remember that p-values just give the probability that an observed difference could be due to random chance (false positives) and do not infer on its biological significance. The lowest p-values do not necessarily correspond to the most meaningful information [START_REF] Wasserstein | The ASA Statement on p-Values: Context, Process, and Purpose[END_REF].

One could filter genes based on fold change (FC) of expression to look at the most discriminant changes in gene expression between groups. Such filtering appears well-suited when dealing with well-distinguished groups such as cells from different lineages characterized by ON-OFF systems of gene expression. However, this filter should be applied carefully when dealing with more continuous variations in gene expression as expected for modifications in cell functioning states. In addition, more subtle changes in gene expression (small FC or coefficient of variation (CV)) might also convey important biological information.

II.8 Biological interpretation and annotation

Once clusters or trajectories are identified, their biological relevance is assessed and a biological identity is assigned by analyzing genes characterizing each cluster or path embranchment of trajectories. This analysis requires external sources of information describing the expected expression profiles of cells with a given identity. For instance, comparing genes characterizing cell groups to literature-derived marker genes or reference databases, one can annotate normal cell types or assign other cell characteristics like their differentiation stage [START_REF] Luecken | Current best practices in single-cell RNA-seq analysis: a tutorial[END_REF]. This strategy relies on well-defined prior knowledge, and can be quite time-consuming.

Another way of analyzing deregulated genes is to group genes based on their biological roles so as to have an indication on the processes or pathways that are enriched in specific cell groups.

Different methods are available to do so: Over-representation Analysis (ORA), Functional Class Scoring (FCS) and Single-Sample (SS) methods. In an ORA, the overlap between deregulated genes and genes in specific gene sets is tested based on assumption that the probability of this overlap will follow a hypergeometric law. For example, gene sets can be biological processes (BP) or molecular functions (MF) from the Gene Ontology (GO) database [START_REF] Ashburner | Gene ontology: tool for the unification of biology. The Gene Ontology Consortium[END_REF] (The Gene Ontology 2017) or pathways from KEGG database (Kanehisa et al 2017). A link between a cell group and a given list of genes can then be established. Such tests can be implemented using publicly available web-based tools such as DAVID [START_REF] Huang | DAVID Bioinformatics Resources: expanded annotation database and novel algorithms to better extract biology from large gene lists[END_REF] and enrichR (Chen et al 2013) (Kuleshov et al 2016) (Xie et al 2021b) or R packages like clusterProfiler [START_REF] Yu | clusterProfiler: an R package for comparing biological themes among gene clusters[END_REF]. ORA methods are relatively simple, reliable and fast to implement [START_REF] Geistlinger | Toward a gold standard for benchmarking gene set enrichment analysis[END_REF]. However, given that genes in a gene set may be co-expressed and that the hypergeometric test is based on the assumption that differentially expressed genes are independent, ORA methods tend to result in higher numbers of false positives than other methods [START_REF] Geistlinger | Toward a gold standard for benchmarking gene set enrichment analysis[END_REF]. The second group of methods, FCS, requires all genes to be first ranked whether or not they are differentially expressed. In general, genes are ordered by direction and magnitude of expression change. Based on the ranked gene vector, FCS methods then assess whether genes in a gene set accumulate at the top or bottom of the gene vector (or more generally, accumulate unevenly along the gene vector). The most widely used FCS method is Gene Set Enrichment Analysis or GSEA [START_REF] Subramanian | Gene set enrichment analysis: a knowledge-based approach for interpreting genome-wide expression profiles[END_REF]. A key shortcoming of FCS methods is that they require pre-ranked gene lists, which may bias data interpretation towards genes with large CV across the datasets considered while genes with small CV can also be functionally important. The last group of methods, the Single-Sample (SS) methods, uses the observed gene expression levels to compute a gene set score in each individual sample and test the association between a phenotype and the sample-level gene set scores with classical statistical models. One commonly used single-sample method is the Single Sample GSEA or SSGSEA [START_REF] Barbie | Systematic RNA interference reveals that oncogenic KRAS-driven cancers require TBK1[END_REF]. Studies have been carried out to compare the performance of enrichment tests in terms of prioritization (i.e. their ability to rank relevant gene sets at the top), sensitivity (i.e. produce small pvalues for the relevant gene sets), specificity (i.e. their false positive rate), computational runtime and functional relevance (i.e. link between top-ranked gene sets and investigated functional state) (Tarca et al 2013) [START_REF] Geistlinger | Toward a gold standard for benchmarking gene set enrichment analysis[END_REF]. They reported that ORA ranks better than FCS and SS methods and is among the most performant tools. While all these methods prove to be useful to gain general insight on main deregulations, it should be remembered that the yielded results are only suggestive and have to be independently validated for the biological object under scrutiny. Indeed, the list of genes associated to a given BP or pathway derive from the aggregation of experimental evidences obtained in different mammalian species, organs, cell types and physio-pathological conditions, and also in different ways (e.g. experimental demonstration or only correlation-based).

Of note, to ensure recognition of analyzed genes in the different databases, one should ensure that gene symbols in the analyzed list are the same as in the databases.

Cancer data poses additional challenges. Relevant reference databases are not available.

Besides, using external sources of information are not adapted for annotating cell functioning states in cancer for two main reasons. Firstly, there is no clear well-established marker gene for identifying There are different ways of calculating gene signature scores. Tirosh and colleagues proposed a scoring method that calculates the arithmetic mean expression of a specific gene set and then, normalizes it by the arithmetic mean expression of reference genes (Tirosh et al 2016a). One drawback of this approach is that the reference genes are identified based on their average expression across the whole dataset, making the method dependent on dataset composition. In addition, this method implies disregarding the ambiguity of the zero value significance, i.e. real lack of expression or detection failure (see § II.3.d). Another way of calculating signature score is to use the ratio of the sum of read counts for all genes in a signature and the sum of all detected read counts in a cell [START_REF] Pont | Single-Cell Signature Explorer for comprehensive visualization of single cell signatures across scRNA-seq datasets[END_REF]. However, comparing signature scores calculated with this method is made difficult by the variability of number of genes in each signature. We recently proposed to use

II.11 Correlation between gene expression and patient clinical outcomes

In cancer research, it is common practice to show clinical relevance of findings, for example by showing if expression of a gene is associated to patient survival. However, this type of analysis suffers important technical bias. Expression levels in bulk RNA-seq data correspond to an average expression of genes across millions of cells in a tumor fragment, therefore crushing crucial information notably when dealing with tumors as heterogeneous as GBM. Deconvolution strategies have been developed to overcome to some degree this problem (e.g. CIBERSORT [START_REF] Newman | Robust enumeration of cell subsets from tissue expression profiles[END_REF], MuSiC (Wang et al 2019d)). For example, CIBERSORT allowed identification of immune cell components from bulk RNA-seq data (Chen et al 2018). Another bias is that RNA-seq data corresponds to the transcriptomes of the majority of cells in a fragment of the whole tumor at the time of diagnosis. Thus, the dynamics of gene expression in time and space are ignored. Since the data only provide a temporal and spatial transcriptional snapshot, it may not be clinically relevant to do correlation analysis between gene expression from such data and patient survival. Moreover, the biological significance of a few months difference in overall survival or progression-free survival is questionable, when considering that patients suffering of GBM have at best a 20-months median overall survival when benefiting from all therapeutic options (Stupp et al 2017). Such analyses can also be highly biased by the type of samples included in the datasets. Monitoring the updates of the TCGA RNA-seq data, we recently found that 8 of the 153 IDH wt -GBM samples in the bulk RNA-seq dataset from the Human Protein Atlas had been relabeled IDH mut -GBM (updated data accessible on (Gliovis)). Repeating previous analyses published by us (Saurty-Seerunghen et al 2019) and others (Gimple et al 2019b) with the updated dataset, we found that the negative association previously observed between our gene of interest and patients' survival was no longer true. This turn of events driven by such a low number of mislabeled samples is likely to stem from the huge difference in life expectancy between patients suffering from GBM bearing a wild-type or mutant form of IDHs (Brennan et al 2013) [START_REF] Ceccarelli | Molecular Profiling Reveals Biologically Discrete Subsets and Pathways of Progression in Diffuse Glioma[END_REF]. The limitation of the relevance of such correlative analyses is further highlighted when considering that two laboratories demonstrated in an independent manner and with distinct approaches the essential role played by this gene in sustaining GBM cell aggressiveness in vitro and in vivo, and provided evidences for its involvement in the patients' tumor setting. The prognostic value of gene expression should therefore be considered with caution.

II.12 Genome-wide GBM analyses at the single-cell level

The first single-cell GBM study was published by Patel and colleagues in 2014 [START_REF] Patel | Single-cell RNA-seq highlights intratumoral heterogeneity in primary glioblastoma[END_REF]. Using scRNA-seq data from five GBM patients, they showed that malignant cells from the same tumor are characterized by different transcriptional profiles, and that the different GBM subtypes identified by Verhaak and colleagues co-exist within the same tumor. They also observed heterogeneity in stemness, oncogenic signaling and immune response in malignant cells. More recently, scRNA-seq data were obtained from multiple regions of several GBM tumors [START_REF] Yu | Surveying brain tumor heterogeneity by single-cell RNAsequencing of multi-sector biopsies[END_REF]. Analysis of data from three GBM showed that each tumor region contains cells with differing CNV profiles, and that the combination of CNV profiles varied between regions of the same tumor.

Moreover, comparison of transcriptomic subtypes observed in eight GBM showed that each tumor region contains cells from the different transcriptomic subtypes, and that the proportion of the different transcriptomic subtypes varies between regions of the same tumor. No clear relation between CNV profile and transcriptomic profile was established in this study. Two other studies also confirmed intra-tumoral genomic heterogeneity in GBM using single-cell technologies. Francis and colleagues published a single-cell DNA sequencing method for identifying unique non-overlapping subclonal alterations and validated their method using in vitro and in vivo functional studies [START_REF] Francis | EGFR variant heterogeneity in glioblastoma resolved through single-nucleus sequencing[END_REF]. They showed that within one GBM, multiple variants of a single RTK gene were expressed and different variants characterized different subclonal GBM cell populations, which contributes to the diversification of deregulated signaling pathways and intra-tumoral heterogeneity in GBM.

Combining scRNA-seq and bulk exome-seq data from three GBM, Muller and colleagues reconstructed phylogenetic trees of clonal lineages for EGF-and PDGF-driven GBMs and identified a cell-type hierarchy during tumor evolution, alike to that observed in early brain development [START_REF] Muller | Single-cell sequencing maps gene expression to mutational phylogenies in PDGF-and EGF-driven gliomas[END_REF].

The majority of other single-cell studies designed exploited scRNA-seq to address GBM cell heterogeneity by assessing the analogy of cancer cells with neurodevelopmental cell types. Yuan and colleagues found that GBM cells resemble both neural and non-neural lineages and that the neural lineage identity is related to GBM cell proliferation [START_REF] Yuan | Single-cell transcriptome analysis of lineage diversity in high-grade glioma[END_REF]. Neftel and colleagues showed that GBM cells exist in four developmental states namely neural progenitor-like (NPC-like), oligodendrocyte progenitor-like (OPC-like), astrocyte-like (AC-like), and mesenchymal-like (MES-like)

states (Neftel et al 2019). Combining scRNA-seq to uniquely-barcoded single cells in vivo, they also demonstrated that GBM cells in different neurodevelopmental states could switch between the different states, thus highlighting cellular plasticity in GBM. In contrast to the plasticity model in GBM proposed by Neftel et al, findings from Couturier and colleagues rather support a hierarchy model (Couturier et al 2020). Using scRNA-seq data from 16 GBM, the authors showed that GBM develops along a conserved neuro-developmental tri-lineage hierarchy (astrocytic, neuronal, oligodendrocytic)

with glial progenitor cells (GPC) being at the apex of the hierarchy. They also demonstrated that different GBM cell populations had distinct proliferative, tumorigenic and chemo-resistant properties. A hierarchy model was supported by another GBM scRNA-seq study, in which trajectory modeling revealed a differentiation hierarchy, with GSCs at the top [START_REF] Castellan | Single-cell analyses reveal YAP/TAZ as regulators of stemness and cell plasticity in Glioblastoma[END_REF]. Gene regulatory network modeling allowed to identify master transcriptional regulators of the GSC state, notably transcriptional coactivators YAP1 and WWTR1 (also called TAZ). The authors experimentally showed that YAP/TAZ is required for establishment and maintenance of the stemness state in GBM.

Bhaduri and colleagues used scRNA-seq data to explore the composition of cell types present in GBM tumors and created a GBM cell-type atlas (Bhaduri et al 2020). They also identified an outer Radial Glial (oRG)-like GBM cell population with enhanced invasive potential and increased PTPRZ1 expression. Another paper also reported the presence of radial glial-like cells in adult GBM and confirmed this finding using scRNA-seq (Wang et al 2020). Analyzing scRNA-seq data from seven primary GBM patients, Richards and colleagues revealed that two major axes of variation can explain heterogeneity in GBM: a differentiation axis between stem-like and differentiated astrocyte-like tumor cells and a developmental-injury response gradient (Richards et al 2021). The authors reported that the first axis is most probable since they found most GBM cells to be organized along an orthogonal astrocyte maturation gradient while maintaining GSC transcriptional profiles.

scRNA-seq studies of gliomas other than GBM also mainly led to the identification of neurodevelopmental-like lineage hierarchies. Venteicher and colleagues showed that IDH mut gliomas (oligodendrogliomas and diffuse astrocytomas) follow a cellular hierarchy model with neural progenitor-like (NPC-like) cells at the top of the hierarchy, which differentiate into two cell subpopulations: OPC-like and AC-like [START_REF] Venteicher | Decoupling genetics, lineages, and microenvironment in IDH-mutant gliomas by single-cell RNA-seq[END_REF]. Filbin and colleagues showed that H3K27M gliomas affecting infants also seemed to follow the cellular hierarchy model, with two main differences compared to IDH mut gliomas [START_REF] Filbin | Developmental and oncogenic programs in H3K27M gliomas dissected by single-cell RNA-seq[END_REF]. First, there may be a differentiation block in H3K27M gliomas, since the majority of cancer cells exhibited a transcriptional profile akin to undifferentiated progenitors. Second, they identified OPC-like cells but not NPC-like cells at the top of the hierarchy in H3K27M gliomas. Another scRNA-seq study also reported differentiation block in pediatric brain tumors [START_REF] Jessa | Stalled developmental programs at the root of pediatric brain tumors[END_REF]. The authors found that normal cell lineage hierarchies are recapitulated in pediatric brain tumors, with differentiation of neural progenitors being impaired.

Two very recent studies explored GBM cell heterogeneity in a non-neurodevelopmental perspective. Analyzing scRNA-seq GBM data based on pathway activity rather than on gene expression profiles as done up to now, Garofano and colleagues uncovered not only a neurodevelopmental axis in GBM but also a metabolic axis [START_REF] Garofano | Pathway-based classification of glioblastoma uncovers a mitochondrial subtype with therapeutic vulnerabilities[END_REF]. They identified two neurodevelopmental states (proliferative/progenitor and neuronal) as well as two metabolic states (mitochondrial and glycolytic/plurimetabolic). They found that the metabolic cellular states were characterized by specific transcriptomic and genetic profiles, which translated into distinct metabolic dependencies. While the glycolytic/plurimetabolic state relied on aerobic glycolysis and amino acid and lipid metabolism for energy production, the mitochondrial state was sustained only by oxidative phosphorylation. Johnson and colleagues investigated the sources of intra-tumoral heterogeneity in gliomas, including GBM, by integrating single-cell DNA methylomes, scRNA-seq and and single-cell copy number profiles with bulk multi-omics data from patient tumors [START_REF] Johnson | TGF-beta is an inducer of ZEB1dependent mesenchymal transdifferentiation in glioblastoma that is associated with tumor invasion[END_REF]. They found that in gliomas, genomic alterations precede epigenetic diversification and promote epigenetic instability. Intra-tumoral epigenetic variability linked to DNA-methylation contributes to gene expression dysregulation and transcriptional heterogeneity, and may facilitate the adoption of distinct functional states (stem-like and differentiated states) in response to stress stimuli (e.g. therapy-induced stress). Thus, epigenetic heterogeneity in GBM may provide a plastic intermediate between genetic subclones and functional states.

Other single-cell GBM studies focused on studying TME. Darmanis and colleagues were the first to explore this aspect in GBM using single-cell data. They analyzed scRNA-seq data from tumor core and peritumoral samples from four GBM patients to assess the effect of the TME on immune cell populations in GBM as well as to characterize the infiltrating malignant cells (Darmanis et al 2017). Another study also used scRNA-seq data to depict the unique immune landscape in GBM [START_REF] Cui | Single-Cell Transcriptomics of Glioblastoma Reveals a Unique Tumor Microenvironment and Potential Immunotherapeutic Target Against Tumor-Associated Macrophage[END_REF]. Combining scRNA-seq and CITE-seq, Pombo Antunes and colleagues explored the heterogeneity of myeloid cells in GBM (Pombo Antunes et al 2021). They found that the myeloid compartment in GBM is large and diverse, with dendritic cells and macrophages being conserved across species (human versus mouse) and dynamic across disease stages (primary versus recurrent).

They also highlighted the heterogeneity of tumor-associated macrophages (TAM) in terms of spatial distribution, transcriptional profiles and cytokine secretion. Muller and colleagues had also explored TAM heterogeneity in gliomas, including GBM, using scRNA-seq, highlighting differences between TAM corresponding to blood-derived macrophages or to microglia [START_REF] Muller | Single-cell profiling of human gliomas reveals macrophage ontogeny as a basis for regional differences in macrophage activation in the tumor microenvironment[END_REF]. Chen and colleagues identified a pro-tumor subpopulation of macrophages in GBM characterized by MARCO expression using scRNA-seq data [START_REF] Chen | Single-cell characterization of macrophages in glioblastoma reveals MARCO as a mesenchymal pro-tumor marker[END_REF]. Mathewson and colleagues used scRNA-seq data to characterize the transcriptional and clonal landscape of T-cells in GBM [START_REF] Mathewson | Inhibitory CD161 receptor identified in glioma-infiltrating T cells by single-cell analysis[END_REF].

This characterization resulted in the discovery that antibody-mediated CD161 blockade enhances T cell-mediated killing of glioma cells. Caruso and colleagues inferred interactions between GBM cells and immune cells (myeloid and T-cells) or oligodendrocytes from scRNA-seq data and reported potential processes that are regulated by these interactions [START_REF] Caruso | A map of tumor-host interactions in glioma at single-cell resolution[END_REF]. Xie and colleagues were the first to characterize the heterogeneity of endothelial cells (EC) in GBM using scRNA-seq data (Xie et al 2021a). Analyzing data from four human glioblastoma and paired tumor peripheral tissues, they found that there exist different EC populations in GBM, with distinct transcriptional and phenotypic profiles and associated with different anatomical locations within and around the tumor.

Moreover, analysis of spatial scRNA-seq and single-cell proteomics data from GBM patient tumors provided evidences of the role of the TME in shaping transcriptional heterogeneity in GBM as well as of the dynamic adaptation of cellular states (lineage states (OPC-like, Radial-Glial-like, NPC-like) and reactive states (reactive hypoxia and reactive immune)) and spatial relationships within the TME [START_REF] Ravi | Functional Subclone Profiling for Prediction of Treatment-Induced Intratumor Population Shifts and Discovery of Rational Drug Combinations in Human Glioblastoma[END_REF]. Furthermore, using scRNA-seq data from five GBM patients across four experimental models (GBM spheres, GBOs, GLICO (GLIoblastoma Cerebral Organoids) and PDX), Pine and colleagues demonstrated that a neuroanatomically accurate human microenvironment is crucial and sufficient for recapitulating the cellular states found in human primary GBMs [START_REF] Pine | Tumor Microenvironment Is Critical for the Maintenance of Cellular States Found in Primary Glioblastomas[END_REF].

There are also other single-cell studies from GBM models including GBM cell lines and organoids (e.g. [START_REF] Krieger | Modeling glioblastoma invasion using human brain organoids and single-cell transcriptomics[END_REF], [START_REF] Hara | Interactions between cancer cells and immune cells drive transitions to mesenchymal-like states in glioblastoma[END_REF]), and mouse models (e.g. [START_REF] Ochocka | Single-cell RNA sequencing reveals functional heterogeneity of glioma-associated brain macrophages[END_REF], [START_REF] Alexander | Multimodal single-cell analysis reveals distinct radioresistant stem-like and progenitor cell populations in murine glioma[END_REF]). I focused on the abovementioned studies because they all used data obtained from patient tumors, which provide the closest view of the original GBM tumors.

Finally, scATAC-seq data from GBM patient tumors have become available recently.

Combining scRNA-seq and scATAC-seq data analysis, Wang and colleagues showed that there exists a gradient of diversity of proliferative GBM cells ranging from proneural to mesenchymal with stemlike cells at the extremes (Wang et al 2019b). Using scATAC-seq data from four GBM patient tumors, Guilhamon and colleagues reported heterogeneity of putative cancer stem cells in GBM in terms of chromatin accessibility profiles as well as functional diversity (invasive, constructive and reactive states) [START_REF] Guilhamon | Single-cell chromatin accessibility profiling of glioblastoma identifies an invasive cancer stem cell population associated with lower survival[END_REF]. Joint analysis of scRNA-seq and scATAC-seq data obtained at two time points of a GBM patient tumor (primary at the time of initial biopsy and recurrent at the time of autopsy) allowed to identify and compare transcriptional network states characterizing cell subpopulations in primary and recurrent GBM [START_REF] Park | A single-cell based precision medicine approach using glioblastoma patient-specific models[END_REF].

Altogether, most single-cell transcriptomic GBM studies done up to now led either to the understanding of GBM physio-pathology through analogies with neuro-developmental lineage hierarchies or to the exploration of the contribution of the tumor environment on GBM development.

I. Capture at the single-cell level of metabolic modules distinguishing aggressive and indolent glioblastoma cells

I.1. Introduction

As stated at the beginning of this document, the aim of my thesis work was to identify eventual GBM metabolic weaknesses at the single-cell level in articulation with the cell functioning states present in patients' tumors at the time of diagnosis. The underlying idea is to identify therapeutic targets of relevance for a broad range of patients despite the inter-and intra-tumor genetic heterogeneity. Growth of all tumors depends on the ability of cancer cells to proliferate, spread across the host tissue(s), and adapt to environment changes occurring along tumor development, regardless of their founding mutations and other genomic alterations. This provided the rationale for focusing on cell functioning states that are expected to be commonly found in all GBM. Likewise, targeting metabolism elements is expected to be effective across differing genomic landscapes, mobilization of metabolic activities being a mandatory step between any signaling pathway and cell behavior.

We chose to start from scRNA-seq because they give access to tumor cell diversity, and because they offer the closest view currently available of the original patient tumors at the time of diagnosis. We used single-cell transcriptomic datasets available when starting the studies. Of note, using publicly-available scRNA-seq datasets favors the reproducibility of the results and enhances their robustness since they are provided by independent laboratories, and contain data of an already high richness that keeps growing (as for now, tens of patients, thousands of cells, millions of gene expression values). By profiling gene expression in individual cells rather than across millions of cells as done in bulk sequencing data, single-cell transcriptomic data are endowed with the potential to deliver information on a cell functioning state and its underlying molecular networks.

To the best of our knowledge, our studies are the first ones to take advantage of single-cell transcriptomic data from patient tumors for identifying GBM cell functioning states in articulation with their associated metabolic pathways. They are based on the two following postulates: (1) a cell functioning state can be inferred from its transcriptomic profile; (2) metabolic pathways at play in a given cell or cells' group can be reconstructed and their key elements identified from the cell transcriptomic profiles.

Up to now, scRNA-seq analyses of GBMs and other cerebral tumors focused on the most dispersed (Darmanis et al 2017) or most expressed genes (Patel et al 2014) (Tirosh et al 2016b) (Venteicher et al 2017) [START_REF] Filbin | Developmental and oncogenic programs in H3K27M gliomas dissected by single-cell RNA-seq[END_REF] or on gene meta-modules highlighted by cell hierarchical clustering done tumor per tumor (Neftel et al 2019). The authors focused on the identification of cell lineage likenesses and genomic anomalies rather than cell functioning states. We therefore thought that applying unsupervised clustering methods on the full set of genes detected in the dataset would highlight at least in part cell functioning states. This turned out to be wrong. GBM cells group predominantly by tumor of origin while normal cells from different tumors are grouped together by cell type, an issue overlooked in the literature. Of note, the fact that normal cells from distinct GBM group by cell type, independently from tumor of origin, shows that scRNA-seq from different tissue samples can be analyzed together without prior batch effect correction.

We therefore studied first the source of the specific influence of the GBM tumor of origin on cancer cell grouping. The results are included in the following article and the extension of this study to other cancers is presented in part 2 of the Results section. We then developed a data reduction approach based on a signature of a functionally coherent set of genes for identifying cell functioning states. This approach was inspired by previously reported uses of sets of genes with correlated expression in specific cell types from normal tissues to identify proliferating cells and stem cells (Patel et al 2014) (Tirosh et al 2016b). We tested the approach using a molecular signature composed of five transcription factors (ARNT2, OLIG2, POU3F2, SALL2, SOX9), each demonstrated experimentally to be required for GBM cell tumorigenicity, that is, GBM cells' ability to initiate tumors (Ligon et al 2007) (Suva et al 2014) (Bogeas et al 2018). The coherence of this signature is further reinforced by the presence of ARNT2 binding sites in the regulatory regions of OLIG2, POU3F2 and SOX9, and the decreased expression in OLIG2, POU3F2, SALL2, SOX9 upon ARNT2 knockdown ((Bogeas et al 2018), and our unpublished results). This strategy allowed cell grouping into two contrasting groups, one with low, the other with high tumorigenic potential, and the modeling of metabolic pathways prevailing in tumorigenic GBM cells. This led to delineate a novel metabolic pathway and pinpoint its key element of potential importance for GBM cell tumorigenicity.

Importantly, we further provided in vitro and in vivo experimental evidences demonstrating that this key element, ELOVL2, a downstream enzyme in the synthesis pathway of polyunsaturated fatty acids (PUFA), is required for GBM cell tumorigenicity.

Altogether, the results obtained show that the two starting postulates onto which the study was based stand true; cell functioning states and their associated metabolic deregulations can be identified from cells' transcriptomic profiles. The detailed results are presented within the published article presented in the next pages. They provided the basis for the development of a more refined clustering strategy that resulted in the discovery of a key metabolic element, this time found to be essential for GBM cell motility (part 3 of the Results section).

I.2. Conclusions

This study led to several findings relevant for computational analysis of cancer cell transcriptomes as well as for the understanding of GBM pathophysiology.

First, it revealed that the tumor identity is encoded by information dispersed throughout the cell transcript repertoires. We now obtained evidences that such an encoding extends to a variety of brain and non-brain neoplasms (see part 2 of the Results section). Second, it showed that information regarding a cell functional state, as well as the associated metabolic pathways, could be extracted from the cell transcriptome profiles. A data-driven approach allows to unmask information related to cell functioning states from single-cell transcriptomes.

As for GBM pathophysiology, it resulted in the discovery of the essential role of the metabolic pathway ensuring synthesis of very-long PUFA, and more specifically of ELOVL2 an enzyme located downstream of this pathway, in controlling GBM cell tumorigenicity. Further computational modeling and experimental assays suggest that ELOVL2 is involved in the formation/release of extracellular vesicles mediating cell-cell communications, known to affect GBM growth (Fareh et al 2017) [START_REF] Broekman | Multidimensional communication in the microenvirons of glioblastoma[END_REF]. This finding is in accordance with the altered phospholipid composition in cell membrane, as found by metabolite profiling of GBM cells upon ELOVL2 knockdown, and reported in an article published during the writing of this article (Gimple et al 2019b). With an approach fully distinct from ours, this study identified ELOVL2 as a gene associated to GBM specific superenhancers and controlling glioblastoma cell stemness. The publication of this article reinforces the importance of ELOVL2 for GBM. It also highlights the strength of our strategy for identifying metabolic enzymes crucial for cell functioning by starting from the tumor setting as present in patients. Finally, it shows that a simple computational analysis of single-cell gene expression data is as efficient as a complex approach such as the one implemented by Gimple and colleagues, which is based on bulk chromatin profiling (chromatin immunoprecipitation followed by next-generation sequencing or ChIP-seq) of fresh tissue tumor samples and GBM cell cultures.

Our methodological strategy is therefore adapted to the discovery of key druggable elements of metabolic modules whose therapeutic targeting is likely to overcome the variability of genetic mutation patterns within and between tumors. So, we set out to develop it and verify its validity for identifying metabolic players crucial for another cell functioning state essential for GBM growth, cell motility (see part 3 of the Results section).

II. The identity of the tumor of origin of cancer cells is encoded by information spread throughout each single cell transcriptome

Introduction

We aimed at identifying metabolic vulnerabilities underlying cell functioning states that are expected to be commonly found in all GBM despite the inter-and intra-tumor genetic heterogeneity. As shown in our previous work, using unsupervised without a priori clustering analyses, GBM cells group predominantly by tumor of origin while normal cells coming from the tumors predominantly group by cell type ((Saurty-Seerunghen et al 2019), part 1 of Results section). Contribution of the different tumors to each cluster was scored by computing Normalized Mutual Information (NMI). This metric is used as an external clustering validation that compares clustering results to a known truth (Manning et al 2008). We used it to compare clustering results to tumor labels. A NMI value of 1 implies that clusters gather objects (cells) corresponding to a single label, whereas a value of 0 denotes that all labels are split across all clusters (i.e. each tumor contributes to each cluster). We tried to identify the source of the tumordriven cancer cell grouping so as to overcome it and unmask information related to cell functioning states. Since normal cells cluster independently from their tumor of origin, scRNA-seq technical biases are unlikely to underlie the dependency of cancer cell clustering on the tumor of origin. In spite of this, we tested the effect of some biases and confirmed that tumor-driven GBM cell grouping is not due to scRNA-seq technical biases. We next found that biological differences between tumors, linked to molecular specificities, are not responsible for this phenomenon either. Our results indicate noteworthily that inter-tumor biological differences cannot be alleviated when taking into account genomic alterations through CNVs and over-dispersed genes, although mutations and CNV are expected to alter the global transcriptome. Moreover, a cancer cell grouping driven by genomic alterations would be expected to highlight genomic subclones rather than the tumor of origin, in coherence with the coexistence of GBM cells with different genomic alterations within one tumor (Nobusawa et decreases with decreasing number of randomly-selected genes demonstrated that this postulate stands true, at least for the GBM cells coming from four tumors composing the dataset we analyzed ((Saurty-Seerunghen et al 2019) and Figure 19A). Here, we further questioned the existence of dispersed information encoding specific identities in additional GBM, other brain and non-brain cancers as well as in normal cells.

Methods

Unsupervised analyses using down-sampled numbers of randomly-selected genes were performed as previously described using three additional independent scRNA-seq datasets from adult For normal cells, we took advantage of the Darmanis dataset we previously used that contains more than 2000 normal cells and more than 1000 cancer cells derived from four GBM.

All datasets downloaded included single-cell transcriptomes of malignant cells from multiple patient tumors (Appendix 3). We performed iterative HCPC analyses on decreasing numbers of genes randomly selected among all genes detected in at least 3 malignant cells as previously described (Saurty-Seerunghen et al 2019). Ten analyses of distinct sets of randomly-selected genes were performed for each size of gene sets (n = 5000, 4000, 2000, 1000, 500, 250, 100, 50, 25, 10 and 5). Of note, gene lists considered in these analyses are dissimilar from one another, with <30% overlap for gene set sizes of 250-5000 and almost no overlap for gene set sizes ≤100 (Appendix 4). Contribution of each lineage subtype -for normal cells -or tumor -for cancer cells -to each cluster was evaluated using the NMI score, relative to NMI score of the grouping analysis performed with all detected genes.

Results and Discussion

As previously shown by several laboratories and us (see Fig1b in Saurty-Seerunghen et al 2019), normal cells from different tumors group predominantly by lineage subtype (e.g. oligodendrocytes, astrocytes, myeloid cells). Repeating unsupervised analyses with a decreasing number of randomlyselected genes on normal cells resulted in progressive reduction of the NMI score (Figure 19B). These results show that normal cells contain, like GBM cells, information dispersed throughout their transcriptome that, in this case, encodes the identity of their lineage subtype.

With respect to cancer cells, we first extended our exploration of the existence of a tumor identity card to additional GBM datasets. Coherently with our previous results, relative NMI scores of the clusterings of GBM cells decrease progressively with the reduction in the number of randomly-selected genes considered in the analyses of all datasets (Figure 20A-D). The influence of the tumor was suppressed only when reducing the number of analyzed genes to less than 50, as shown by the stabilization of NMI scores. Identical results were obtained when analyzing other types of gliomas (Figure 20E-H), i.e. oligodendroglioma, astrocytoma (diffuse astrocytoma, IDH-mutant in the WHO classification, Louis et al 2016) and pediatric gliomas (GBM and H3K27M-glioma). These results show that information encoding the identity of the tumor of origin can be found in all brain tumor cells, regardless of the brain cancer subtype and of the sequencing technology used.

To test the robustness of this information, we determined whether or not it is maintained in individual cells from different regions of a given GBM tumor. Using the results of our previous unsupervised analysis of the Darmanis dataset, we highlighted GBM cells coming from either the tumor core or its periphery on a UMAP representation. As shown in Figure 21, cancer cells from the core and periphery grouped together, indicating that the information related to the tumor identity card is maintained despites differences in the cells' environment.

Next, we determined whether a tumor identity card is also encoded in the transcriptome of cancer cells coming from non-brain tumors. We performed unsupervised analyses using down-sampled number of randomly-selected genes on scRNA-seq data from melanoma, breast cancers, head and neck (Figure 22). These results show that all cancer cells retain an identity card of their tumor of origin.

Altogether, our results show that all malignant cells from a given tumor share a common tumor identity card which is encoded by information dispersed throughout the cells' transcriptome, regardless of the tumor type considered. Our analysis of a GBM cell dataset containing cancer cells from two distinct tumor regions suggests that this information is robust, being maintained across totally distinct (Osswald et al 2015) (Broekman et al 2018) [START_REF] Dominiak | Communication in the Cancer Microenvironment as a Target for Therapeutic Interventions[END_REF]. Shared information could ease the establishment and the efficiencies of these interactions.

Whatever be its source and physical support, the information revealed by clustering analyses taking into account the whole cell transcriptome appears related first to the uniqueness of the evolution of cells, whether it is the evolution of normal cells along their lineage specification or the evolution of cancer cells along the tumor development. Consequently, extracting, from the cells' transcriptomes, more specific biological information -for example the functioning state in which cells are -demands to perform computational analyses on a priori selected gene subsets.

III. Human glioblastoma cell motility depends on activity of the cysteine catabolism enzyme 3-Mercaptopyruvate sulfurtransferase

We explored a second GBM cell functioning state contributing to malignancy and poor prognosis of GBM, cell motility. As discussed in the general introduction, GBM cells are highly infiltrative -they can invade the healthy brain parenchyma and can migrate far away from the original tumor, thus making complete surgical resection of GBM impossible. Moreover, currently used irradiation and chemical therapies have been shown to promote GBM cell invasiveness (Huber et al 2013) (Lu et al 2012). In this study, we aimed at capturing metabolic modules sustaining cell motility, and also at identifying metabolic elements underlying the cell transition in or out of this functional state.

We applied our signature-driven data reduction approach to highlight subpopulations of cells in a motile state present in all tumors in independent single-cell transcriptomic datasets. The genes composing the motility signature were selected based on their requirement for GBM cell motility being experimentally demonstrated in vitro and in vivo with recent patient-derived glioblastoma cells (PDC). We upgraded our previous grouping method (Saurty-Seerunghen et al 2019) in order to obtain more circumscribed groups of cells, the analyses now being based on the expression levels of each of the motility signature genes, and integrating UMAP and HCPC approaches. A motility score integrating expression of each signature's elements was used to identify two homogeneous clusters having the highest and lowest motile potential in independent datasets. We also integrated cell trajectory modeling in our analytical workflow so as to identify genes crucial for the switch from low to high motility.

We found that GBM cells with high motile potential are characterized by higher energetic needs and higher oxidative stress than cells with low motile potential, as well as an over-representation of the methionine and cysteine metabolic pathway. Cell trajectory modeling highlighted an element of the cysteine catabolism pathway -the enzyme 3-Mercaptopyruvate sulfurtransferase (MPST) -at the crossroad of the path from low to high motility. MPST requirement for GBM cell motility was experimentally demonstrated. Further biochemical explorations indicated that MPST acts by protecting protein cysteine residues from hyperoxidation. The detailed results are presented under the form of an article manuscript in the next pages.

These results further show the potential soundness and the power of our analytical strategy, which allows discoveries of targetable metabolic enzymes crucial for GBM cell functioning.

Introduction

Glioblastoma (GBM) cell heterogeneity is now recognized as a major issue to overcome in order to develop more efficient therapies for this primary brain tumor that remains with a very poor prognosis. The median survival of patients does not exceed twenty months despite aggressive treatments combining surgical removal of the tumor with chemotherapy and radiotherapy (Stupp et al 2017). GBM is a paradigm of intra-tumor heterogeneity. Normal and cancerous cells intermingle, and tumor territories differing in tumor cell densities, vascular perfusion and necrosis entangle [START_REF] Hambardzumyan | Glioblastoma: Defining Tumor Niches[END_REF]. Genome-wide DNA and RNA sequencing at the tissue and single-cell levels have further illustrated GBM heterogeneity at the genomic and ontogenic levels by showing the co-existence of cancer cells bearing differing mutational loads and genomic rearrangements and varying ontogenic similarities with neuro-developmental lineages (Snuderl et (Richards et al 2021). Notwithstanding these variable mutation loads, chromosome rearrangements, and developmental-like profiles, we reasoned that GBM heterogeneity would appear less overwhelming when considering the behavior or functional state of tumor cells in link with their underlying metabolic pathways. Tumor cells in identical functional states can be detected within each GBM and across patients, regardless of their genomic backgrounds. Downstream of all signaling pathways regulating cell functioning states is metabolism. Studies have shown that changes in metabolic activity not only accompany changes in cell functioning states but can also drive the changes (Flavahan et al 2013) (El-Habr et al 2017) (Oizel et al 2017). Thus, metabolic enzymes playing key roles in the adoption and/or maintenance of cell functional states crucial for tumor development, as observed at the time of patients' diagnosis, constitute therapeutic targets able to overcome tumor heterogeneity.

Cell motility is a key contributor to GBM malignancy and poor prognosis. GBM cells are highly infiltrative, making complete surgical resection of the tumor impossible. Moreover, currently used irradiation and chemical therapies have been shown to promote GBM cell invasiveness (Huber et al 2013) (Lu et al 2012). Cell motility encompasses two processes: migration, the cell's ability "to move around in a space that is freely available", and invasion, a process including microenvironment remodeling by the cells (de Gooijer et al 2018). Motile cells have to develop protrusions in the direction of migration, disrupt adhesion sites at the cell rear, degrade and remodel the extracellular matrix (ECM) and ultimately make new connections at the cell front by remodeling their cytoskeleton so that they are pushed forward, making cell motility a highly integrated process (Armento et al 2017). Three major migration/invasion routes have been described: the narrow and tortuous extracellular space of the brain parenchyma (Scherer 1938), the perivascular spaces surrounding blood vessels (Scherer 1938) (Montana & Sontheimer 2011) (Watkins et al 2014) (Cuddapah et al 2014) and white matter tracts, notably the corpus callosum, used as a highway for bilateral brain hemisphere invasion (Scherer 1938) (Pedersen et al 1995). Both single cell and collective migration have been reported (Vollmann-Zwerenz et al 2020) (Volovetz et al 2020), the later appearing to be facilitated by the formation of an interconnected network of tumor microtubes (Osswald et al 2015). So far, an integrated view of the metabolic pathways at play in a motile GBM cell is lacking.

Here, we applied a signature-driven data reduction approach to highlight subpopulations of cells in a motile state present in all tumors from four independent single-cell transcriptomic datasets. We previously demonstrated that cell functional state and their associated key metabolic players can thereby be uncovered from transcriptomes of single cells derived from surgical resections of the patients' tumors (Saurty-Seerunghen et al 2019). Computational analyses integrating cell trajectory modeling characterized cells with high motile potential as endowed with higher energetic needs and oxidative stress than their more static counterparts, and highlighted the 3-Mercaptopyruvate sulfurtransferase (MPST) enzyme at the crossroad of the path from low to high motility. Requirement of this element of the cysteine catabolism for GBM cell motility was experimentally demonstrated. Mechanistically, MPST acts by protecting protein cysteine residues from hyperoxidation. Our findings uncover a previously unsuspected involvement of MPST in the switch from low to high cell motility and maintenance of GBM cell motility.

Material and Methods

All figures were prepared using Adobe Illustrator. All bioinformatics analyses were performed using the R software version ≥ 3.6.1 (https://cran.r-project.org/). All resources and materials, R packages, corresponding websites and references are listed in Table S1.

scRNA-seq data pre-processing

Four publicly available single-cell transcriptomic datasets from GBM patient tumors were downloaded, including the only two publicly-available datasets obtained with the SMART-seq2 technique (Neftel et S1). The last two datasets were selected over others because of the availability of the raw UMI counts. They are thereafter designated as N-S, D-S, N-10X and PA-10X respectively. Only good-quality cancer cells were included in the analyses. Cells from Tumor105B-D in the N-10X dataset were excluded because of their lower numbers of detected genes and UMI counts per cell than the other tumor cells. Good-quality cancer cells were identified based on number of genes and transcripts detected per cell and also on percentage of mitochondrial genes for 10X datasets. For D-S dataset, cells with <1700 genes and <90000 transcripts were considered as being lowcomplexity transcriptomes. For 10X datasets, cells with <150 genes, <250 UMI counts and >20% mitochondrial genes were considered as being low-complexity transcriptomes. Therefore, we retained for further analyses 4916 malignant cells from twenty adult GBM patients in the N-S dataset, 1033 malignant cells from four adult GBM patients in the D-S dataset, 5797 malignant cells from six adult GBM patients in the N-10X dataset and 8666 malignant cells from six adult GBM patients in the PA-10X dataset. Log2-transformed Counts Per Million (log2(CPM+1)) were used for D-S and the 10X datasets, unless otherwise specified. CPM corresponds to the counts of gene-mapped reads normalized by the total number of mapped reads per cell library divided by one million, thus allowing comparison of read abundance across libraries of different sizes. For the N-S dataset, gene expression data currently available are in log2((TPM/10) + 1). TPM (Transcripts Per Million) normalization takes into account an eventual biased estimation of long transcript numbers, by dividing the number of mapped reads by the transcript's length. Since expression values in CPM could not be calculated from these data, log2((TPM/10) + 1) values were just transformed into log2(TPM + 1) to optimize the comparison between the different datasets. To avoid potential analytical bias due to scarcely detected genes, genes detected in less than 1% cells were filtered out.

Malignant and normal cells were distinguished either according to the cell annotations provided for both SMART-seq2 datasets, or when absent as for 10X Genomics datasets based on detection of copy-number variations (CNVs), a hallmark of malignant cells. Expression data in log2(CPM/100 + 1) values were used, which shrinks expression ranges, allowing to keep only the major gene expression variations, based on which malignant and normal cells can be discriminated. Data were processed using a three-step approach: CNV inference, marker gene expression and unsupervised cell clustering. CONICSmat R package used to infer CNVs has the advantage of not requiring an a priori reference cell dataset (Muller et al 2018). The default filtering and normalization procedures were followed, as outlined in https://github.com/diazlab/CONICS/wiki/Tutorial---CONICSmat;---Dataset:-SmartSeq2-scRNA-seq-of-Oligodendroglioma. CONICSmat fits a two-component Gaussian Mixture Model (GMM) to the average gene expression across all cells within each chromosomal region. As a result, genes of a given region present a lower expression in cells with a deletion of this region than cells without the deletion. The opposite is observed for genes of amplified regions. The posterior probabilities for each cell belonging to one of the two components of the model are then calculated. The copy number status across cells is predicted by the posterior probabilities for each cell belonging to the component with the higher mean. For each region, a CONICSmat likelihood ratio test adjusted p-value <0.001 and a difference in Bayesian Criterion >300 were retained as CNVs. The GMM-based CNV predictions were then used to group cells into potential malignant and non-malignant groups and were visualized using heatmaps (ComplexHeatmap R package (Gu et al 2016)) and UMAP (umap R package) plots based on 500 most variable genes (Fig. S1A-B, Fig. S2A-B). Second, expression of marker genes for pan-immune cells (PTPRC/CD45), macrophages (ITGAM, FCGR3A/CD16A, CD14), microglia (CSF1R, TMEM119), T-cells (CD2, CD3D) and oligodendrocytes (MOG, MAG) was highlighted on UMAP plots (Fig. S1C, Fig. S2C). Finally, a hierarchical clustering followed by a K-means clustering on the UMAP components (FactoMineR R package (Husson et al 2010)) was applied to identify cell groups that were most similar or different to one another (Fig. S1D, Fig. S2D). Combining CNV status prediction, marker gene expression and clustering results, cells harboring CNVs, clustering together on UMAP and devoid of normal cells markers were considered as malignant (Fig. S1E, Fig. S2E).

Cell grouping analyses

Clustering analysis was based on a molecular signature of ten elements (see results) and achieved using the classical Hierarchical Clustering on Principal Components (HCPC) approach (FactoMineR package (Husson et al 2010) modified to implement in a stepwise manner four methods for multivariate data analyses: Principal Component Analysis (PCA), Uniform Manifold Approximation and Projection (UMAP), hierarchical clustering, and partitioning clustering by the k-means method (Fig. S3A). The ten principal components (PCs) identified by PCA were reduced non-linearly to two dimensions using UMAP. Euclidean method was used to construct a cell-to-cell distance matrix based on the two UMAP components. Hierarchical clustering was then performed on this distance matrix using the Ward's criterion (ward.D2 algorithm) in order to determine the number of clusters. The resulting partitioning of the cells was improved by a K-means clustering with ten iterations. The cell grouping was visualized using UMAP (umap R package) or chord plots (circlize R package (Gu et al 2014)). A Normalized Mutual Information (NMI) score (ClusterR package (Manning et al 2008)) was calculated to determine the contribution of cells issued from distinct tumors to each cluster, as described (Saurty-Seerunghen et al 2019). Of note, THBS1 was excluded from the motility signature when applied to PA-10X dataset because it was detected in <1% GBM cells in that dataset.

Motility and other scores were obtained by computing the geometric mean of the expression values per cell of each element of the molecular signature corresponding to a given score. For each signature, only genes detected in ≥ 1% cells in ≥3 datasets were considered. When null, expression values were imputed a value of 1.

Analyses of genes differentially expressed between cell groups

Genes differentially expressed between cell groups with differing scores were identified following Mann-Whitney (Wilcoxon Rank Sum) test with p-values adjusted for multiple testing (Benjamini-Hochberg (BH), p-value < 0.01) (Soneson & Robinson 2018). Fold change (FC) was calculated as follows: FCi=xi-yi, where xi and yi are the log2 expression levels of gene i in conditions x and y, respectively. Only genes detected in at least 3% of GBM cells were considered for this analysis. Genes coding for metabolism enzymes were identified using the list from KEGG (Kanehisa et al 2017).

Comparing list of differentially expressed genes between distinct datasets was performed with venn R package on genes detected in ≥3% malignant cells in all datasets, and after updating gene symbols based on gene metadata files downloaded from HGNC and NCBI websites (Table S1). Genes absent from these metadata files, or with ambiguous symbols (e.g. genes whose current approved symbol is the previous symbol of another gene), or whose symbol was previously associated with more than one gene were excluded from the lists.

Gene ontology analyses were carried out on enrichR website (Chen et al 2013) (Kuleshov et al 2016) (Xie et al 2021b). We considered BH-adjusted p-value < 0.05 (Fisher's Exact test) as the cut-off criterion for significance. Redundant terms were excluded. Graphs were generated using Prism 7.0 (Graphpad) software.

Trajectory inference analysis

To model the path taken by cells with low motile potential to reach a high motile potential, STREAM python package was used (Chen et al 2019). Briefly, the expression values of the ten elements of the motility signature were first extracted. Next, dimensions were reduced to four components using the spectral embedding algorithm (dimension_reduction STREAM function, other parameters are set to default). The components were then used for simultaneous tree structure learning and fitting using ElPiGraph (seed_elastic_principal_graph STREAM function, default parameters). The resulting trajectory structure was represented in a 2D subway map plot, where straight lines represent branches and each dot represents a single cell (plot_stream_sc STREAM function, dist_scale = 0.5 and root = 'S5', other parameters to default). To identify genes differentially expressed between cell populations from two adjacent branches, the detect_de_markers STREAM function was used. The mean of scaled gene expressions in each branch was calculated. Then, the fold change (FC) of mean expression between pairs of branches was computed, and Mann-Whitney U test performed. The statistic U is then standardized to Z-score to assess in which branch the gene is overexpressed. Genes with Z-score greater than 1 and log2FC greater than 0.15 are considered as genes over-expressed between two adjacent branches. Multiple-testing correction was performed using the Benjamini-Hochberg method. Significance level was set at q value < 0.05.

For this analysis, only the N-S dataset was used because of the deeper sequencing depth of SMART-seq2 data compared to 10X Genomics data, and because of the higher number of malignant cells and distinct patient tumors compared to D-S dataset.

Biological material, lentiviral transduction and pharmacological treatment

Patient-derived cells (PDC) 6240**, R633, 5706** and P3 obtained in our and other laboratories from neurosurgical biopsy samples of distinct primary GBM were cultured in defined medium containing bFGF and EGF or bFGF only (P3), as described (El-Habr et al 2017) (Rosenberg et al 2017) (Eskilsson et al 2016). The PDC stably express luciferase, and 5706**-PDC express GFP as well. Cells were transduced with lentiviral vectors encoding a shRNA construct either neutral (shControl) or targeting MPST transcripts (shMPST) (Table S1). Non-transduced cells were eliminated following a 1-2 week treatment with puromycin (1-2 µg/mL). The lentiviral particles were produced by the Plateforme vecteurs viraux et transfert de gènes (Necker Federative structure of research, University Paris Descartes, France). GBOs were generated and cultured from a primary GBM patient tumor (Table S1) following the method developed by Jacob and collaborators (Jacob et al 2020). In relevant experiments, cells were treated with 200µM I3-MT-3 or 1mM NAC (MedChemExpress and Sigma respectively) or their vehicles (DMSO and culture media respectively).

Cell migration and invasion assays

Cell migration was assessed using spheroid-on-Matrigel assay, as described (Guyon et al 2020). Briefly, cells were counted as described (El-Habr et al 2017), and cell spheroids generated by seeding 2500 cells/100µL 0.4% methylcellulose (Sigma) in a U-bottom well (96 wells plate, Falcon). After a 24-48 hours culture, the spheres were seeded on Matrigel (0.2mg/mL culture media, Corning), and cell migration assessed after 7 hours (5706**), 24 hours (R633, P3) or 7 days (6240), unless specified otherwise. Images were acquired with a Zeiss microscope. The total area covered by GBM cells migrating from a spheroid or GBO, as well as the central spheroid core area, were delineated and measured with FIJI software. The migratory index calculated corresponds to the ratio between the total area and the spheroid core area. Values were normalized to shControl or Vehicle. Microfluidic chips were also used to evaluate cell migration. They were designed and fabricated, as described (Courte et al 2018). PDMS (Polydimethylsiloxane) blocks were modeled to contain 3mm-long and 50 µm-high reservoirs opening on two successive rows of 4.5µm-high straight microchannels linked by 40μm long arches (arch separation = 3µm). The PDMS blocks were bonded to glass coverslips using an oxygen or air plasma treatment to form a microfluidic chip and sterilized using UV. 2*10 5 cells PDC expressing shControl or shMPST were seeded at a 1:1 ratio into the cell reservoirs. The shMPST construct encoding GFP, in addition to the shRNA, allowed distinction between shControl and shMPST cells. Seven days post-seeding, cells were fixed with 4% paraformaldehyde, and nuclei stained using DAPI (Sigma). Cells were visualized with a fluorescent microscope (Zeiss). The distance travelled by the 10 cells having migrated into the microchannels the furthest away from the wall of the channel entrance was measured using FIJI software in each microfluidic chip.

Cell invasion was assessed using spheroid-in-collagen assay (Guyon et al 2020) and Matrigelcoated transwells (Renault-Mihara et al 2006). Spheroids generated as described above, were seeded in collagen (1mg/mL culture media, Corning), and cell invasion assessed 24 hours later. Images were acquired with a Zeiss microscope. The total area covered by GBM cells leaving the spheroid and invading the collagen, as well as the central spheroid core area, were delineated and measured with FIJI software. The invasive index calculated corresponds to the ratio of the total area and the spheroid core area. Values were normalized to shControl or Vehicle. Transwells (8µm pores, Corning) were coated for 30min with Matrigel (0.2mg/mL culture media, Corning). They were then placed in wells containing culture media with growth factors. 40000 (6240**, R633 and P3) and 20000 (5706**) cells were seeded in 100µl media without growth factors into the transwells, and further cultured for 24 hours. Cells having invaded the bottom part of the membrane were counted after paraformaldehyde fixation and nuclei staining with DAPI, as described (Renault-Mihara et al 2006).

Detection of intracellular ROS levels

Cells or spheroids were incubated at 37°C for 30min with CellROX Green (5µM/R633, 10µM/6240** and P3) or Deep Red reagent (1.25µM/P3, 2.5µM/R633 and 5706**) (ThermoFisher, France). Fluorescence signal was then assessed by microscopy (CellROX Green, Zeiss microscope) or by FACS (CellROX Deep Red, ARIA II, BD Biosciences). When indicated, cells were treated with a ROS scavenger, N-Acetyl Cysteine (NAC, 1h, 37°C) and/or a ROS generator, Menadione (Mena, 30min, 37°C), and ROS levels were measured by FACS analysis. The latter was performed on 10000-gated events using Vdye-labeled cells for detecting any non-viable cell (1µL/4-5*10 5 cells/mL, 30min, 4°C, eBiosciences). Violet laser (405nm) and Pacific Blue filter were used for Vdye detection, Blue laser (488nm) and FITC filter for CellROX Green reagent detection, and Red laser (640nm) and APC filter for CellROX Deep Red reagent detection. Data analysis and figure generation were performed using the FACS Diva software (BD Biosciences).

Measure of mitochondrial mass

Cells, spheroids or GBO were incubated with 1µM MitoTracker Green reagent (Invitrogen) at 37°C for 30min and the fluorescence signal was assessed by microscopy or FACS. For FACS analysis, it was performed on 10000-gated events. Violet laser (405nm) and Pacific Blue filter were used for Vdye detection, and Blue laser (488nm) and FITC filter for MitoTracker Green reagent detection. Data analysis and figure generation were performed using the FACS Diva software (BD Biosciences).

Gene expression analysis

Gene expression knockdown in response to shRNA expression was verified by RT-QPCR as previously described (Saurty-Seerunghen et al 2019) using the LightCycler480 (Roche, France) and the SYBR Green PCR Core Reagents kit (Bimake.com). The thermal cycling conditions comprised an initial denaturation step at 94 °C for 5 min, and 40 cycles at 94°C for 30sec, 60°C for 30sec and 72°C for 30sec. Transcripts of the TBP gene encoding the TATA-box binding protein (a component of the DNA-binding protein complex TFIID) were quantified as an endogenous RNA control. Quantitative values were obtained from the cycle number (Cq value), according to the manufacturer's manuals. Sequences of primers used for Q-PCR are given in Table S1.

Protein analyses

Cells were lysed using HEN Buffer (0.5M EDTA, 20% SDS, 1% NP-40, 10mM Neocuproine and 0.1M HEPES, pH 7.4) supplemented with 5mM 4-chloro-7-nitrobenzofurazan (NBF-Cl, Sigma-Aldrich), as described (Zivanovic et al 2019). Proteins were then precipitated twice in the presence of methanol and chloroform at a centrifugation of 14000rpm for 15min at 4°C, the pellets being resuspended in 50mM HEPES (pH 7.4) supplemented with 2% SDS. Once completely dissolved, protein concentration was determined by BCA assay and adjusted to 3mg/mL. For detection of MPST protein levels, 30μg protein extracts were separated by SDS-PAGE using 4-12% (w/v) NuPAGE gels (Invitrogen, France) and transferred to Transblot nitrocellulose membranes (Biorad, USA). MPST detection was achieved using immunoblotting with anti-MPST (Sigma, 1:4000). The secondary antibody was anti-rabbit IgG (GE Healthcare, 1:30000). Signal detection was performed with SuperSignal™ West Femto Maximum Sensitivity Substrate chemiluminescence detection system (ThermoFisher, France). Densitometric analysis was achieved using ImageJ software. MPST levels were assessed by normalizing the chemiluminescent signal by the amount of proteins loaded assessed using Red Ponceau. Persulfidated protein levels were assessed following the Dimedone-Switch protocol described by Zivanovic and colleagues (Zivanovic et al 2019). Briefly, this method relies on the indiscriminate labeling of persulfides, thiols, sulfenic acids, and amino groups with 4-chloro-7-nitrobenzofurazan (NBF-Cl, the NBF-adducts yielding a fluorescent signal detected at 488nm). Upon incubation with dimedone-based probes, only the persulfides exchange their labeling, thus enabling their identification (red signal detected at 800nm). The protein extracts previously prepared in presence of NBF-Cl were incubated for 30min at 37°C with 25μM click mix provided by Dr Milos Filipovic (1mM Daz-2, 1mM Cy7.5 alkyne, 2mM TBTA Cu Complex, 4mM Lascorbic acid, 30% acetonitrile and 20mM EDTA in PBS) and precipitated as described above. Proteins were subsequently resuspended in 50mM HEPES supplemented with 2% SDS and separated using SDS-PAGE on 4-12% (w/v) polyacrylamide gels. Gels were fixed in 12.5% methanol and 4% acetic acid for 30min. The Cy7.5 signal corresponding to persulfidated proteins was recorded at 800nm and the NBFadducts signal was recorded at 488nm on a Chemidoc Imager (BioRad, France). Persulfidation levels were assessed by normalizing the Cy7.5 intensity to Alexa488 (NBF-adducts) intensity.

Intracerebral xenografts

The animal maintenance, handling, surveillance and experimentation were performed in accordance with and approval from the Comité d'éthique en expérimentation animale Charles Darwin N°5 (Protocol #5379). R633, P3 and 5706** PDC transduced with lentiviruses encoding either a shControl or a shMPST, were used. 1.5*10 5 cells were injected stereotaxically into the striatum of anesthetized 8-10 weeks-old Nude female mice (Janvier Laboratories, France), using the following coordinates: 0mm posterior and 2.5mm lateral to the bregma, and 3mm deep with respect to the surface of the skull. Tumor formation and development were monitored by bioluminescence imaging performed on a Photon Imager Biospace (Biospace Lab, France), after intra-peritoneal injection of 150µL luciferin (20mM, ThermoFisher, France). Bioluminescent signals were visualized with M3 Vision software (Biospace Lab, France).

Statistical analyses

R version 3.6.1 or Prism 7.0 software (GraphPad) were used to generate plots and for statistical analyses. Significance level was set at p < 0.05, unless otherwise indicated (*p < 0.05). The type of statistical test used is provided in the figure legends. All experiments were performed using independent biological samples. Mean ± SD are shown.

Results

Grouping cells according to their potential motile state

We used four independent publicly available single-cell transcriptome datasets generated from surgical resections of patients bearing IDH-wildtype GBM. Two were obtained with the SMART-seq2 technology (Neftel et al 2019) (Darmanis et al 2017) and two with the 10X Genomics technology (Neftel et al 2019) (Pombo Antunes et al 2021). They are hereafter designated as N-S, D-S, N-10X and PA-10X, following the initial of the first author of the paper reporting its first description and the technology used. The computational analyses were first performed using the N-S dataset, and the robustness of the results evaluated by repeating the analyses with the other three datasets. This choice stems from the deeper sequencing depth of SMART-seq2 data compared to 10X Genomics data, and the higher number of malignant cells and distinct patient tumors in the N-S dataset compared to the D-S dataset.

For clustering cells in a motile state, we implemented a data reduction approach driven by a molecular signature (Fig. 1a), an approach we previously reported to group cells from distinct tumors according to their functional state (Saurty-Seerunghen et al 2019). Here, cells were grouped based on the combined expression levels of each gene in a motility signature, using the HCPC algorithm (Hierarchical Clustering on Principal Components) (Husson et al 2010) modified to integrate also UMAP components (Uniform Manifold Approximation and Projection) (Fig. S3A). The genes composing the motility signature were selected based on the in vitro and in vivo experimental demonstration of their requirement for GBM cell motility using recent patient-derived glioblastoma cells (PDC) cultured in defined media. These inducers and effectors of GBM cell motility are involved in pro-migratory autocrine signaling (TGFB1, SMAD3, THBS1 [START_REF] Joseph | TGF-beta is an inducer of ZEB1dependent mesenchymal transdifferentiation in glioblastoma that is associated with tumor invasion[END_REF] (Daubon et al 2019)), morphological rearrangements undertaken by motile cells and linking extracellular guidance signals with the remodeling of the cytoskeleton (ACTN4 (Sen et al 2009) (Ji et al 2019) (Tentler et al 2019), PTK2 (Xia et al 2016) (Frolov et al 2016), PXN (Frolov et al 2016) (Sun et al 2017) (Lopez-Colome et al 2017), TLN1 (Kang et al 2015), VCL (Xia et al 2016)) and remodeling of the matri-cellular environment (TNC (Xia et al 2016), SPARCL1 (Gagliardi et al 2017)). Each signature element is detected in ≥3 % malignant cells and a majority presents a statistically significant positive correlation of their expression across all cells (Fig. 1b, Table S2).

This clustering analysis resulted in six cell clusters (Fig. 1c) with differing expression profiles of the signature genes (Fig. 1d). Of note, taking into account UMAP components in addition to Principal components (PC) enhanced the homogeneity of transcriptomic profiles between cells of a given group, as shown by increased Silhouette index (Fig. S3B). The Normalized Mutual Information (NMI) value was then calculated to determine the contribution of different tumors to each cluster, as described (Saurty-Seerunghen et al 2019). A NMI value of 1 implies that clusters gather cells corresponding to a single tumor label, whereas a value of 0 denotes that each tumor contributes to each cluster. As indicated by the low NMI value (NMI = 0.11), each cluster was composed of cells coming from distinct tumors (Fig. 1e). Using the geometric mean of the expression values of each of the signature gene as a motility score, allowed to distinguish the two cell clusters with the highest (Mig HIGH ) and lowest (Mig LOW ) motile potential, respectively (Fig. 1f). Differential expression analysis between Mig HIGH and Mig LOW cell groups (Mann Whitney, BH-adjusted p-value < 0.01) highlighted the genes overexpressed in Mig HIGH (Table S3-Sheet 2).

Enrichment of motility-related genes in Mig HIGH cells

Next, we evaluated the coherence of our findings by confronting them with experimentally- acquired previous knowledge. Finding enrichment in motility-related terms upon gene ontology (GO) analyses performed on the genes overexpressed at least two-fold in Mig HIGH cells (log2 fold-change≥2) provided the first support of the relevance of the analytical strategy implemented (Fig. 1g, Table S4-Sheet2). To further probe this relevance, we applied the same analytical strategy to the other three selected datasets. THBS1 being detected in <1% malignant cells in the PA-10X dataset, was excluded from the motility signature for that dataset. Each of these analyses resulted in clusters composed of cells coming from different tumors (NMI =0.058-0.15) (Fig. S4A-C). Differential expression analysis between Mig HIGH and Mig LOW cell clusters (Mann Whitney, BH-adjusted p-value < 0.01) (Table S3-Sheets3-5), followed by GO analyses of the most overexpressed genes, again highlighted enrichment in motilityrelated terms (Fig. S5, Table S4-Sheets3-5). To determine the similarities of Mig HIGH cell transcriptional profiles of distinct datasets, we compared the different lists of overexpressed genes. The highest overlaps were found between lists of genes derived from datasets obtained with the same sequencing technique. 69.2% (4852/7010) of the genes overexpressed in Mig HIGH cells from N-S dataset were also overexpressed in Mig HIGH cells from D-S dataset, and 86.0% (6627/7707) of the genes overexpressed in Mig HIGH cells from N-10X dataset were also overexpressed in Mig HIGH cells from PA-10X dataset (Fig. S4D).

Comparing the four lists of genes overexpressed in Mig HIGH GBM cells, we found that 43-53% of genes overexpressed in Mig HIGH GBM cells from one dataset were also overexpressed in Mig HIGH cells from the three other datasets (Fig. S4D, Table S3-Sheet6). These results favor the idea that the implemented analytical strategy captures cells with related transcriptional profiles from distinct datasets.

To further determine the extent to which our analytical approach captures cells in a motile state, we selected three gene modules previously found to be associated with GBM cell motility (Table S3- Sheet7), and determined their enrichment in Mig HIGH cells and their correlation with the motility score. Most motile cells, including GBM cells, having been shown to undergo transcriptional reprograming akin to the so-called epithelial-mesenchymal transformation (EMT) process, we first used an EMT gene module of 14 genes associated with mesenchymal transformation of neural cells and GBM cell motility (Carro et al 2010) (Armento et al 2017). We also sought for enrichments in two other gene modules recently shown to be mobilized in invasive GBM cells. One was highlighted by Bhaduri and colleagues who identified an outer radial glia (oRG)-like malignant cell population with increased invasive behavior in GBM, and delineated 36 oRG-like associated genes with highly correlated expression in their GBM scRNA-seq dataset (Bhaduri et al 2020). The other is involved in the regulation of cell motility in several cancers and driven by TEAD transcription factors (Thompson 2020). TEAD1 role in the regulation of GBM cell motility was demonstrated by Tome-Garcia and colleagues who identified 32 TEAD1-regulated genes involved in GBM cell motility (Tome-Garcia et al 2018). Each of these three gene modules was found to be enriched in Mig HIGH cells from each analyzed dataset (Fig. 1h). We also observed a positive linear relationship between the motility score and the EMT, oRG and TEAD scores (Fig. 1i, Fig. S6), and higher scores for each of these modules in Mig HIGH cells compared to Mig LOW cells from each dataset (Fig. 1j, Fig. S6).

Altogether, these results support the relevance of the analytical strategy implemented for capturing cells with high motile potential from distinct datasets.

Metabolic characteristics of cells with high motile potential

We went on to identify metabolic pathways deregulated in GBM cells with high motile potential. The genes coding for enzymatic components of the metabolic pathways overexpressed in Mig HIGH cells from each dataset were identified using the 2019 KEGG list of 1586 metabolic elements (Kanehisa et al 2017) (Table S5-Sheet2). KEGG pathway analyses highlighted enrichment in metabolic pathways related to extracellular matrix (ECM) modeling (e.g. glycosaminoglycan synthesis/degradation) and membrane composition rearrangements (e.g. unsaturated fatty acids (FA), FA elongation, GPI-anchor biosynthesis) (Fig. 2a, Table S5-Sheets3-6). In addition to these metabolic pathways expected to be mobilized in motile cells, we also observed enrichment in pathways involved in energy production (Fig. 2a, Table S5-Sheets3-6). Upregulation of genes involved in glycolysis, amino acid catabolism, fatty acid degradation, protein recycling (lysosome), TCA cycle, oxidative phosphorylation, suggests that GBM cells with high motile potential have higher energetic needs than cells with low motile potential. Enrichment in oxidative phosphorylation, a term gathering elements of the electron transport chain (ETC), associated with enrichment in TCA cycle components, suggests in addition a higher mitochondrial load in Mig HIGH cells. Finally, we observed enrichment in pathways involved in anti-oxidative processes, such as the pentose phosphate pathway (PPP), and the cysteine, sulfur and glutathione metabolisms (Fig. 2a, Table S5-Sheet3-6). Enrichment in pathways or organelles known to be involved in ROS scavenging and/or production of anti-oxidants (porphyrin and selenocompound metabolisms, peroxisome), as well as enrichment in purine/pyrimidine metabolism, known to be required for DNA/RNA repair, are also compatible with enhanced responses to oxidative stress. Overall, the results of this analysis point to an enhanced mobilization of anti-oxidative pathways counteracting the deleterious actions of oxidative compounds derived from the over-mobilization of energy production pathways.

Trajectory modeling was used next to delineate metabolic components potentially crucial for acquisition of a motile cell state. Trajectory reconstruction based on expression of the motility signature genes was achieved with the STREAM python package using the N-S dataset which has the highest sequencing depth. The inferred trajectory corresponded to a branched path, highlighting an intermediary branch (S3-S1 branch) linking cells with low motile potential to cells with high motile potential (Fig. 2b). Among the metabolic enzymes overexpressed in S3-S1 branch compared to each other branch, characterizing therefore the switch between low and high motile states (Table S6), two were also overexpressed in Mig HIGH cells of all datasets considered. NIT2 encodes an omega-amidase involved in anaplerosis through shuttling of α-ketoglutarate and oxaloacetate from the glutamine and asparagine metabolisms into the TCA cycle (Cooper et al 2016). MPST (Mercaptopyruvate Sulfurtransferase) encodes an enzyme involved in cysteine catabolism and is part of anti-oxidative cell defenses. MPST converts 3-mercaptopyruvate, derived from cysteine transamination, to pyruvate while transferring a sulfur to a thiophilic acceptor, thereby forming a persulfide (Filipovic et al 2018). MPST is overexpressed in GBM tissues compared to normal brain tissues (Fig. 2c). MPST was selected for further studies based on the following rationale. Mitochondria-located GOT2, an aminotransferase which ensures the formation of the MPST substrate 3-mercaptopyruvate, is overexpressed in GBM cells with high motile potential (Table S3). MPST participates to anti-oxidant defenses, notably through persulfidation (also known as S-sulfhydration) of proteins (P-SSH), a post-translational modification of reactive cysteine residues that protects proteins from irreversible cysteine hyperoxidation (Filipovic et al 2018) (Pedre & Dick 2021). Persulfidated MPST (MPST-SSnH) has also been proposed as a key intermediate in production of GSH-persulfides (GSSH) (Kimura et al 2017). Although MPST participation in the control of GBM cell properties is unknown, evidences have been provided for its stimulatory effect on the migration of endothelial cells (Tao et al 2017) (Abdollahi Govar et al 2020) and murine colon cancer cells (Augsburger et al 2020). These data and the results from our modeling studies led us to postulate that enhanced oxidative stress coupled with MPST mobilization plays a key role in GBM cell motility (Fig. 2d).

Enhanced ROS production in cells with high motile potential

The predictive values of the results from the computational analyses were experimentally evaluated using GBM patient-derived cells (PDC) exhibiting different relative migratory and invasive properties. Spheroid-on-Matrigel migration assays, invasion-across-collagen assays and Matrigel-coated transwell assays classified 5706**-and 6240*-PDC as the most and least motile, respectively, whereas R633-and P3-PDC exhibited an intermediary motility (Fig. S7). Of note, these different motility profiles determined in vitro are faithful to the histological aspect of the tumors developing upon intra-cerebral xenografts, the most motile cells forming the most diffuse tumors (Bogeas et al 2018) (Daubon et al 2019).

We first determined whether Reactive Oxygen Species (ROS) affect GBM cell motility. Intracellular ROS levels were evaluated with CellROX ® Green or Deep red fluorogenic probes. Cell exposure to the anti-oxidant N-Acetyl Cysteine (NAC) or the ROS generator menadione resulted in decreased or increased fluorescent signals detected by FACS, thereby confirming CellROX ® suitability for detecting ROS in GBM PDC (Fig. S8A). Microscopic observation of spheroids first seeded on Matrigel and then labeled with CellROX ® Green probe showed a higher ROS signal in GBM cells migrating out of spheroids formed by all PDC tested including the ones with the lowest motility (6240**-PDC) compared to their static counterparts (Fig. 3a).

Considering this result with our prior observation of enriched expression in genes coding for the TCA cycle and the ETC, we measured mitochondrial mass using the MitoTracker Green reagent. An increased mitochondrial load was detected in cells migrating out from cell spheroids or from cultured patients' GBM tissue fragments (GBO) (Fig. S8B). FACS was then used to isolate PDC according to their CellROX or Mitrotracker signals. Cells with high ROS production exhibited a 2.1-2.6 times higher mitochondrial mass than cells with low ROS production (Fig. 3b). Conversely, cells with high mitochondrial mass exhibited higher ROS production than cells with low mitochondrial mass (1.9-3.2 times more) (Fig. 3c). We next assessed the invasive potential of GBM PDC with high ROS production. To do so, we FACS-sorted GBM PDC according to their basal ROS production detected with CellROX ® probe (Fig. 3d). Evaluation of the ability of each cell population to invade through Matrigel-coated transwells showed a 1.4-7.7-fold increase in invasiveness of GBM cells with high ROS levels compared to cells with low ROS levels (Fig. 3d). Finally, we examined the effects of decreasing ROS levels on cell motility by exposing to NAC the most motile PDC (5706** and R633, Fig. S7). Decreasing ROS levels with NAC treatment resulted in decreased GBM cells' ability to invade through Matrigel-coated transwells (Fig. 3e). These results indicate that enhanced ROS production resulting from active oxidative phosphorylation is instrumental in GBM cell motility, as predicted by the results of our computational analyses.

MPST enzymatic activity is required for GBM cell motility

To determine if MPST is necessary for GBM cell motility, we knocked down MPST expression in the three GBM PDC with high and intermediate motile potential, using lentiviral transduction of small hairpin (sh) RNA. Cell transduction with shMPST lentiviral constructs led to decreased MPST mRNA and protein levels (Fig. 4a-b) with no major cytotoxic effect (Fig. 4c). Similar results were obtained with another shMPST construct (Fig. S9). Knocking down MPST expression resulted in a sharp decrease in cell migration (51-86%) as determined with spheroid-on-Matrigel assays (Fig. 4d). Monitoring cells' ability to move across 4.5µm-high channels in a microfluidic chip further showed that MPST knockdown also impaired GBM cells' ability to change their shape so as to migrate through narrow and tortuous spaces (Fig. 4e). Cell invasion was also decreased by 44-95% upon MPST knockdown as observed using spheroid-incollagen and Matrigel-coated transwells invasion assays (Fig. 4f-g). The observed decrease in cell motility can be due to an inhibition of MPST catalytic activity or to defects in additional effects mediated by MPST independently from its enzymatic activity. To distinguish between these two possibilities, we treated cells or GBOs with 200 µM I3-MT-3, a selective pharmacological MPST inhibitor (Hanaoka et al 2017). Spheroid-on-Matrigel assays showed that I3-MT-3 decreased PDC migration compared to vehicletreated cells (Fig. 4h). Alike PDC, cell migration out of GBO was also inhibited by I3-MT-3 (Fig. 4h). These results indicated that MPST catalytic activity was required for GBM cell motility. Given that MPST participates to the cells' anti-oxidative defenses through its role in protein persulfidation, we further determined protein persulfidation levels in GBM PDC with decreased MPST expression or activity using the recently published Dimedone-Switch method (Zivanovic et al 2019). Decrease in the global protein persulfidation levels upon MPST knockdown (Fig. 4i), and upon MPST enzymatic inhibition (Fig. 4j), showed that MPST participates to protein persulfidation in GBM cells.

Finally, we evaluated in vivo the consequences of MPST knockdown on tumor development using orthotopic xenografts of PDC stably expressing luciferase and either shControl or shMPST. Bioluminescence imaging showed cell engraftment and initial tumor development whether the cells expressed shControl or shMPST. However, a delayed tumor formation was observed in xenografts performed with 2 of the 3 shMPST-PDC (Fig. 5a). In all cases, we observed a very robust decreased bioluminescent signal in mice engrafted with cells expressing shMPST as compared to shControl, indicating a reduced tumor burden in mice grafted with shMPST-PDC, compared to mice grafted with shControl-PDC (Fig. 5b). Mice survival was monitored for two of the xenograft series, one performed with cells driving either slow -within 6 months -or rapid -within 3 months -tumor development. In both cases, we observed a significant improvement in the survival of the mice grafted with PDC expressing shMPST. Mice xenografted with shMPST-PDC survived on average twice longer than mice xenografted with shControl-PDC (Fig. 5c).

Taken together, these results show that enhanced oxidative stress coupled with MPST mobilization plays a key role in GBM cell motility.

Discussion

The diversity of cancer cells populating GBM is a major challenge in treating these brain cancers. Understanding the molecular basis of the varying functional cell states that co-exist within GBM and that are commonly shared across tumors and as observed at the time of patients' diagnosis, opens a path to overcome this challenge. A way to access to such information is coupling experimental validations with computational modeling of transcriptomes of single cells sorted from patient GBM, which offer the closest possible map of the patients' tumor cells. ) and 45h (GBO) postseeding on Matrigel. Mean ± SD, n = 5-8 independent biological samples, * p <0.05, Mann-Whitney test. i. Decreased protein persulfidation levels in shMPST-PDC vs shControl-PDC. R633 and P3 PDC. Ingel detection of protein persulfidation (P-SSH) using the dimedone switch method with Cy7.5 as a reporting molecule and fluorescent signal of the NBF-adducts (488nm). P-SSH levels calculated as a ratio of Cy7.5/488 fluorescent signals. Mean ± SD, n= 3 independent biological samples, * p <0.05, unpaired t-test. j. Decreased protein persulfidation levels upon cell treatment with the MPST inhibitor I3-MT-3. R633 PDC. In-gel detection and calculation of P-SSH levels as in i. Mean ± SD, n= 3 independent biological samples, * p <0.05, unpaired t-test. metabolic characteristics of motile cells across thirty patients' GBM, and to identify MPST as an antioxidative enzyme crucial for GBM cell motility.

Our results suggest that GBM cells with high motile potential are characterized by higher energetic needs than cells with low motile potential, as shown by the enrichment of several energetic pathways. This is in coherence with motile cells needing more energy to colonize surrounding brain tissue and adapt to new microenvironments with unique nutrient and oxygen availability, and hence mobilizing different pathways to meet their high energetic needs (Garcia et al 2021). Of note, this overall energetic mobilization appears to converge towards the TCA and the ETC, a modeling result supported by the experimental demonstration of enhanced ROS production and higher mitochondrial mass load in motile GBM cells. These results are in coherence with previous studies on epithelial cancer cells showing that an ETC overload with preserved mitochondrial functions and an increased mitochondrial superoxide production promote tumor cell migration and invasion (Porporato et al 2014). Results from studies using GBM cell lines cultured in serum also highlight a role for mitochondria. An example is the trafficking of energetically active mitochondria to membrane protrusions accompanying enhanced invasive properties of the LN229 GBM cell line (Caino et al 2015). Another is the inhibitory effect of ROS scavenging on EGFstimulation of T98G cell migration (Pudelek et al 2020) counteract oxidative stress. Upregulation of the pentose phosphate pathway (PPP) implicates enhanced NADPH production necessary for the function of several antioxidant proteins, whereas enhanced glycine, glutamate and cysteine metabolism together with enhanced glutathione metabolism implicates enhanced availability of the major anti-oxidant compound GSH. Cysteine metabolism involvement in anti-oxidative process extends beyond supporting GSH synthesis. Protein cysteine residues are known targets of peroxide (H2O2), a major ROS species, and undergo H2O2-dependent oxidation or antioxidantdependent reduction. The CSE (Cystathionine Gamma-Lyase also known as CGL, encoded by CTH), CBS (Cystathionine Beta-Synthase) and MPST (also known as MST) enzymes of the cysteine catabolism participate to protein protection from dismal hyperoxidation of their reactive cysteine residues (Zivanovic et al 2019).

MPST was the only metabolic enzyme endowed with anti-oxidant properties predicted by trajectory modeling to be at the crossroad of the path leading GBM cells from low to high motility. Our results demonstrate the biological validity of this prediction. MPST knockdown resulted in a robust inhibition of GBM cell motility. The cells exhibit a decreased ability to move around in a freely available space, as shown with spheroid-on-Matrigel migration assays and microfluidics assays, the latter showing in addition a sharp reduction in the cells' ability to change their shape so as to migrate through narrow spaces. MPST facilitation of cell migration might extend to cell types other than GBM cells, MPST inhibition having previously been involved in the migration of a mouse carcinoma cell line (Augsburger et al 2020), and endothelial cells (Tao et al 2017) (Abdollahi Govar et al 2020). We also observed a robust decrease in GBM cell invasiveness upon MPST knockdown, GBM cells exhibiting an impaired ability to remodel the ECM so that they can move into their surroundings as shown with spheroid-in-collagen and Matrigel-coated transwell assays. In coherence with these in vitro results, alteration of GBM cell properties by MPST knockdown translates in vivo into reduced tumor burden, and a robust increase in mice survival, albeit the cells retain their tumor-initiating properties. MPST enzymatic activity is required for its pro-motile effects, the MPST pharmacological inhibitor I3-MT-3 (Hanaoka et al 2017) having the same consequences as MPST knockdown on GBM cell migration. MPST catalyzes stepwise the desulfuration of 3-mercaptopyruvate that generates an enzyme-bound persulfide, and the transfer of the persulfide's outer sulfur atom to proteins or small molecule acceptors (Pedre & Dick 2021). MPST knockdown was accompanied with reduced persulfidated protein levels in GBM cells, indicating that protein persulfidation is part of the mechanisms by which MPST affects GBM cell motility. Whether the advantage conferred by this post-translational modification stems from a mere overall protection of proteins from hyperoxidative damage or from more specific consequences remains to be determined. The few surveys of persulfidated proteins performed without exogenous dihydrogene sulfide (H2S) supply showed that persulfidation affected proteins involved notably in energetic pathways (Ida et al 2014) (Gao et al 2015) and migration and associated signaling pathways (Zivanovic et al 2019) (Murphy et al 2019) (Zuhra et al 2021). For instance, proteins of the glycolysis pathway were reported to be persulfidated in a human pancreatic beta cell line (Gao et al 2015). The resulting increase in glycolytic metabolic flux could stem from enhanced catalytic activity of glycolytic enzymes such as GAPDH or LDHA under their persulfidated form (Mustafa et al 2009) (Gao et al 2015) (Untereiner et al 2017). Persulfidation studied in HEK293 cells has also been reported to enhance actin polymerization and alter actin-dependent cytoskeletal rearrangements (Mustafa et al 2009). Interestingly, several proteins pertaining to migration, including the focal protein kinase encoded by PTKN2, have also been found to undergo persulfidation in HeLa cells treated with EGF, a known pro-migratory signaling factor (Zivanovic et al 2019). Most of these reported increases in protein persulfidation are presumably mediated by CBS and/or CSE activity, as supported by studies implementing inactivation of these enzymes and/or by their conspicuous expression in the biological systems under scrutiny (Mustafa et al 2009) (Gao et al 2015) (Untereiner et al 2017) (Zivanovic et al 2019) (Zuhra et al 2021). MPST-mediated protein persulfidation has so far been reported for two proteins, the sulfurtransferase Mocs3 and thioredoxin (Trx). Through Mocs3, it is involved in a cascade leading to protein urmylation and tRNA thiolation, and through Trx in H2S generation (Pedre & Dick 2021). The overall decrease in protein persulfidation levels we observed in GBM cell upon MPST knockdown indicates that MPST affects a wider range of proteins.

Sulfur metabolism, considered from the perspective of H2S levels detected in tumor tissue homogenates, has been reported to be lower in GBM than in normal brain tissues, and controlled H2S production appears to be necessary for tumor fitness insofar as increasing H2S levels with H2S donors decreases GBM cell viability, and tumor development from subcutaneous mouse GBM cell grafts (Silver et al 2021). Finding decreased H2S production upon inhibition of CBS expression or activity in mouse GBM cells (Silver et al 2021), and in a serum-cultured cell line (Takano et al 2014), point to CBS as the mediator of H2S production in GBM as in other cancer cells studied so far. In this context, our results unravel a novel role for sulfur metabolism in the control of GBM development through MPST. We believe that disclosure of the participation of this enzyme in the control of GBM cell motility was made possible by taking into account GBM heterogeneous populations through the consideration of the cell functional state in link with its metabolic underliers. Only single-cell data allow the access to analysis at this level of granularity, avoiding confusion brought up by averaging different cell types and functional states when considering a whole tissue or cell culture. However, it remains unknown whether our modeling results pertain to a specific mode of cell motility and/or path of migration. Both individual and collective migrations have been identified in GBM (Vollmann-Zwerenz et al 2020) (Volovetz et al 2020). Considering that oRG move in an individual manner in the developing cortex is in favor of our capturing of cells undergoing solitary displacements. However, it is yet to be demonstrated that gene markers for oRG-like GBM cells are effectively supporting individual migration rather than collective migration. With respect to the path of migration, we found that MPST knockdown impaired cell invasiveness not only through collagen, a major component of the tumor ECM, but also through Matrigel, whose composition mixes its major laminin component to collagen IV and heparan sulfate proteoglycans, similar to the blood vessels' basement membrane found in GBM (Ljubimova et al 2006). These results suggest that MPST could be necessary for GBM cells moving along blood vessels as well as through the brain parenchyma. In the different datasets considered in this study, the edge or core origin of the cells is most often unknown, though they are most likely coming from dense tumor cell regions. Removal of the surrounding parenchyma is usually undertaken with great care by neurosurgeons in order to spare at the best motor and cognitive functions. Although Darmanis and colleagues distinguished malignant cells from the tumor core and the periphery, the cancer cells isolated from the latter are in too small numbers to result in robust computational analyses. There is therefore a high probability that the results of our analyses relate to dense tumor areas' cells, moving either along the dense blood vessel network that characterize these highly angiogenic tumors, or across the tumor parenchyma enriched in components favoring cell displacement. This thesis work aimed at exploring, at the single-cell level, cell functioning states as present in patients' GBM at the time of diagnosis, and to identify their eventual metabolic weaknesses. It brought several contributions to the computational analysis of cancer cell transcriptomes as well as to the understanding of GBM pathophysiology.

We demonstrated that unsupervised without a priori clustering analyses are not adapted for identifying cell functioning states from single-cell transcriptomic data. Indeed, these analyses result in cancer cell grouping predominantly by tumor of origin, and in normal cell grouping predominantly driven by lineage subtypes. Although these observations were not novel, they were not previously considered as an issue to be addressed. We showed that the identity of each patient tumor is encoded by information dispersed throughout the cell transcript repertoires, which masks any other potentially relevant biological information. We showed that this holds true for several cancers, including those with varying founding mutations (e.g. GBM, melanoma) as well as those with more homogeneous ones (e.g. IDH mut gliomas). Furthermore, the results from our computational analyses of normal cells coming from GBM tumors show that an identity card is also encrypted throughout their transcriptome, this time corresponding to their lineage identity. These identity cards appear as a reflection of the cells' evolution, either along their lineage specification for normal cells or along tumor development for cancer cells. Extracting information other than cell identities from the cells' transcriptomes demands therefore performing computational analyses on a priori selected gene subsets. We thus developed a data-driven clustering strategy based on a signature of a functionally coherent set of genes. This strategy allowed to unmask information related to cell functioning states and to reconstruct the associated metabolic pathways at play from single-cell transcriptomes. We identified GBM cell populations with distinct tumorigenic (part 1 of Results section) or motile (part 3 of Results section) potential, uncovered their associated metabolic pathways and experimentally validated the role of a metabolic enzyme crucial for each functioning state. Firstly, we found that cells with high tumorigenic potential are characterized by an enrichment in lipid and amino acid metabolic pathways. We demonstrated that ELOVL2, involved in the synthesis of very-long PUFA, is required for GBM cell tumorigenicity and that it may act through the formation/release of extracellular vesicles mediating cell-cell communication, known to affect GBM growth. Secondly, we found that GBM cells with high motile potential are characterized by high energetic needs and high oxidative stress, including an enrichment in cysteine metabolism. We demonstrated that MPST, involved in cysteine catabolism, is essential for GBM cell motility. Results obtained suggest that MPST acts by protecting protein cysteine residues from hyperoxidation. Altogether, we demonstrated that our methodological strategy is adapted to the discovery of key druggable elements of metabolic modules whose therapeutic targeting is likely to overcome the inter-and intra-tumor variability of genomic anomalies.

To the best of our knowledge, our studies are the first ones to take advantage of single-cell transcriptomic data from patient tumors for identifying GBM cell functioning states in articulation with their associated metabolic pathways. Till date, analysis of single-cell transcriptomic GBM data led mainly to the identification of GBM cell populations resembling normal cell lineages or to deciphering the role of the tumor microenvironment on GBM growth. The originality of our approach lies in the use of molecular signatures to achieve cell groupings according to their potential functional state from unicellular transcriptomes. Molecular signatures have previously been used, once cell clustering has been achieved, for biologically interpreting and annotating identified clusters.

For example, they have been used for assigning cell lineage identities (Tirosh et al 2016b), for assessing functioning states (e.g. stem [START_REF] Patel | Single-cell RNA-seq highlights intratumoral heterogeneity in primary glioblastoma[END_REF] (Tirosh et al 2016b)) or "activities" of biological processes or pathways (e.g. cell cycle (Tirosh et al 2016a) (Tirosh et al 2016b), T-cell activity (Tirosh et al 2016a)). Our results demonstrate that the already-available data can be further leveraged to explore other aspects crucial for understanding GBM physio-pathology, highlighting the high richness of the available single-cell transcriptomic data, and the possibility to use these data to explore different and novel perspectives. It was not so obvious that single-cell transcriptomes could unveil information related to metabolic activities associated to cell functional states for several reasons. Single-cell transcriptomes provide a temporal and spatial snapshot of mRNA abundance in a sample, and may thus not reflect the whole intra-sample heterogeneity, the dynamics of gene expression and the full complexity of cell functioning states. Besides, they do not provide information on RNA turnover dynamics (i.e. transcription and degradation rates) or on post-transcriptional RNA modifications. The latter affects RNA stability, processing, splicing and translation, and has been shown to play pivotal roles in controlling GBM cell functioning [START_REF] Dong | The Emerging Roles of RNA Modifications in Glioblastoma[END_REF]. Therefore, RNA modifications may contribute to cancer cell heterogeneity in GBM, as observed in breast cancers [START_REF] Marcel | Ribosomal RNA 2'Omethylation as a novel layer of inter-tumour heterogeneity in breast cancer[END_REF] or diffuse large B-cell lymphoma [START_REF] Krogh | Profiling of ribose methylations in ribosomal RNA from diffuse large B-cell lymphoma patients for evaluation of ribosomes as drug targets[END_REF]. Most importantly, transcriptomic data only provide an indirect readout for protein levels in a given cell. They do not either provide information on protein stability. For instance, the mRNA levels and the protein levels may not correlate when considering a protein with long half-life. Moreover, they do not capture the numerous possible post-translational modifications affecting protein activity. For instance, a protein can remain inactive as long as it is not phosphorylated. Furthermore, sequencing techniques relying on oligo (dT) priming, such as SMART-seq2 or 10X Genomics, do not capture non-polyadenylated RNAs (e.g. rRNAs) which also regulate protein levels in a cell. For instance, rRNA composition and modifications shape different subpopulations of ribosomes with distinct functions which differentially control translation [START_REF] Monaco | 2'-O-Methylation of Ribosomal RNA: Towards an Epitranscriptomic Control of Translation?[END_REF] [START_REF] Venezia | Emerging Role of Eukaryote Ribosomes in Translational Control[END_REF] and thus, may influence cell functioning.

When analyzing single-cell transcriptomic data, it is also important to remember to take into account the technical biases involved and use the right data for the right question. For instance, currently we have to make a compromise between sequencing depth and number of cells profiled: deeply sequencing a small number of cells (e.g. SMART-seq2 data), or profiling a large number of cells at a shallower sequencing depth (e.g. 10X Genomics data). Data with higher sequencing depth may be more appropriate when studying transcription factors for instance, which are usually expressed at low levels and thus have a higher probability of being missed using shallower techniques. However, they would be less appropriate for identifying rare cell populations given the relatively small number of cells that can be profiled. For example, cells endowed with drug-resistance mechanisms may constitute a minor cell population difficult to identify among few profiled cells. In this case, having more cells to analyze would be recommended. Likewise, one should remember the limits of analytical tools used. A flaw of most clustering methods is that they will cluster data, regardless of whether or not there are any groups present, and whether or not the identified groups are biologically meaningful. For example, well-defined groups may not be present in scRNA-seq data when analyzing a continuum as in differentiation processes. In this case, trajectory analyses may be more adapted than clustering analyses. In this thesis work, before analyses, we verified the clustering tendency of the data using the Hopkins statistic which tests whether or not the data points are uniformly randomly distributed [START_REF] Lawson | New index for clustering tendency and its application to chemical problems[END_REF]. Randomly-distributed data indeed implies that there is no meaningful cluster. Even if clusters are identifiable in a dataset, it is important to assess their biological relevance. In this work, we determined whether or not genes expected to be associated to the cell functioning state under scrutiny were among the genes characterizing the clusters. To this aim, we reviewed the literature and also performed enrichment analyses. Another limit of analytical methods is that they may perform differently depending on the quality of the data analyzed. Though data from both SMART-seq2 and 10X Genomics sequencing techniques resulted in identification of cells with high motile potential with similar transcriptional profiles and of MPST being a potential metabolic vulnerability, we observed that the homogeneity of identified clusters differed between sequencing techniques. 10X Genomics data yielded less homogeneous clusters and more discrete motility scores than SMART-seq2 data (part 3 of Results section). We also observed variations in fold change range across independent datasets when applying the same analytical workflow on four GBM scRNA-seq datasets obtained with two different sequencing techniques (Appendix 5). These variations may stem from differences in capture between sequencing techniques, differences in quality of sequencing, as well as real differences in gene expression levels between samples. The observed variability in fold change range reinforces the idea that thresholds should be set carefully and be adapted to the datasets considered. Moreover, trajectory inference (TI) analyses may also be affected by data quality. TI analytical tools have successfully been used for modeling mainly differentiation trajectories, even with data obtained using shallow sequencing techniques. Since changes in differentiation states involve rather discrete changes in gene expression, shallow sequencing techniques may still be able to capture the gene expression variations. However, when modeling transitions between cell functioning states, which involve combinations of subtle gene expression changes, shallow data may be less appropriate.

The methodological strategies implemented proved fruitful in identifying metabolic vulnerabilities of glioblastoma cells, common to numerous patient tumors irrespective of the interand intra-tumor genomic heterogeneity. A weakness of our strategy is that it requires the definition of a biologically-relevant molecular signature with strong experimental evidences, and hence relies on background knowledge of genes crucial for a given cell functioning. However, these may not be well characterized or may be difficult to find in the literature. Moreover, the defined signatures may be not exhaustive and may be restricted to only some aspects of the cell functioning state considered, which can also bias the results. For example, defining a motility signature using only genes required for collective migration would have allowed to explore only this aspect of cell motility. Finally, the number of genes constituting the signature should be limited in order to avoid predominant cancer cell grouping according to its tumor of origin while being sufficient to be used as a basis for the computational approaches implemented. We have already evolved our strategy in the course of our work to obtain better circumscribed cell groups and to integrate trajectory inference analyses to identify metabolic enzymes required for the switch from one cell functioning state to another. We now envisage to refine existing molecular signatures and to define novel ones for yet poorly explored functional states using correlation-based analyses on current datasets, and machine learning approaches applied on data obtained from discrete biological samples in a homogenous state. The refined signatures would also allow us to test whether subgroups of cells in a given cell functioning states can be identified (e.g. subgroups of cells with high tumorigenic potential) and whether these subgroups are characterized by distinct networks of transcription factors or specific metabolic activities.

In its present state, our analytical strategy allowed to obtain an overview of the metabolic pathways at play in link with GBM cell functioning as present in the patient tumors at the time of diagnosis. It also led to circumscribe potential key metabolic players. The biological consequences of their genetic and/or pharmacological manipulation supported the pertinence and robustness of our methodological approach, ELOVL2 being found to be required for GBM cell tumorigenicity and MPST for GBM cell motility. The role of ELOVL2 in inter-cellular communication has to be further investigated. As for MPST, the identification of the proteins losing their persulfidation mark in response to MPST knockdown should bring further insight on the mechanism at stake in motile cells. This implies the implementation of specific precipitation of persulfidated proteins coupled with mass spectrometry analyses. Determining whether conditional expression of MPST will suffice to push a static GBM cell into moving will help to further verify the biological soundness of the predictions resulting from the trajectory inference modeling. The goal of this research work was to identify metabolic vulnerabilities that could interfere with the development of fully-established tumors, as seen in patients. To demonstrate that targeting a given metabolic enzyme will indeed affect the cell functioning state under scrutiny in fully-established patient tumors, one should reduce the expression of the enzyme or decrease its activity, once tumors are fully grown. In an attempt to do so, we are currently developing a new tool to induce the conditional knockdown of genes of interest in PDC. Of note, the in silico analyses pointed towards several additional potential metabolic vulnerabilities previously unsuspected within the context of GBM. We are currently experimentally evaluating the most promising ones.

Our methodological strategy proved fruitful for identifying metabolic modules at play in glioblastoma cell in a tumorigenic or motile state. We envisage to apply it to other functioning states crucial for GBM growth, notably angiogenesis and therapy resistance. GBM cells may adopt a proangiogenic state to favor nutrient supply within the tumor. The therapy resistant state prevents efficiency of current therapy which includes chemotherapy using temozolomide and irradiation (Stupp protocol). Thus, for this state, we will define two distinct signatures: chemo-resistant and irradiation-resistant. The drug-resistant state in primary GBM could be more challenging to identify in single-cell transcriptomic data than the tumorigenic or motile state since it may constitute a rare (if not inexistent) cell population. This might be because of a low number of cells actually in that state in a primary untreated tumor, because this cell population may be present in tumors at a specific timepoint of tumor evolution (not necessarily present at the time of diagnosis) or because these cells may be present in specific regions of the tumor (not the region from which the analyzed cells come). In addition, if cells in a drug-resistant state in primary GBM are dormant cells with low transcriptional activity, they may be wrongly discarded as low-complexity cells during the cell quality control step.

Finally, analysis of mRNA levels alone is not sufficient to characterize and understand biological systems. This stands true also when investigating cell functioning states. A cell functioning state arises as a result of a combination of transcriptional, proteomic and epigenomic changes in a cell. Thus, analyzing mRNA levels alone provides a partial view of the cell's reality. In this thesis work, we used only single-cell transcriptomic data because of their availability. More recently, scATAC-seq data from GBM patient tumors have become available (Wang et al 2019b) [START_REF] Guilhamon | Single-cell chromatin accessibility profiling of glioblastoma identifies an invasive cancer stem cell population associated with lower survival[END_REF]) [START_REF] Park | A single-cell based precision medicine approach using glioblastoma patient-specific models[END_REF]. These data could be analyzed from the same perspective as we did for single-cell transcriptomic data to have a better understanding of GBM cell functional heterogeneity and plasticity. We could leverage them to identify the epigenetic regulators crucial for establishment and/or maintenance of specific cell functioning states in GBM. We could also determine if some cells are primed to adopt certain functioning states (e.g. a drug-resistant state), despite their plasticity.

Moreover, analysis of single-cell epigenomic data could help determine how each tumor imprints and maintains its identity card within the whole transcriptome (part 2 of Results section). As omics technologies progress, simultaneous profiling of the genome, transcriptome, proteome and/or epigenome in the same cell will become more accessible, and the availability of multi-omics data will allow obtaining a more integrated view of GBM cells. Their appropriate leveraging is bound to allow important progresses in our understanding of this deadly cancer with an unprecedented analytical finesse that will hopefully translate in improved therapies for patients.

Figure 1 :Figure 2 :Figure 3 :Figure 4 :Figure 5 :Figure 6 :Figure 7 :Figure 8 :

 12345678 Figure 1 : Primary brain and other central nervous system (CNS) tumors ........................................... 18 Figure 2 : Glioma classification based on histopathological features ................................................... 19 Figure 3: Classification of diffuse gliomas based on histological and genetic features ........................ 20 Figure 4: Subgroups of glioblastoma .................................................................................................... 22 Figure 5: Glioblastoma microenvironment ........................................................................................... 24 Figure 6: Schematic representation of mechanisms involved in GBM cell motility. ............................ 28 Figure 7: Main genomic alterations in GBM ......................................................................................... 31 Figure 8: Models of GBM cell heterogeneity ........................................................................................ 33 Figure 9: "Warburg effect" in proliferating tissues and cancer cells .................................................... 36 Figure 10: Analyzing molecular profiles of cell populations using bulk and single-cell approaches .... 42 Figure 11: General scRNA-seq experimental workflow ........................................................................ 43 Figure 12: General analytical workflow for analyzing scRNA-seq data from complex tissues. ............ 46 Figure 13: Schematic representation of data integration strategies. ................................................... 50 Figure 14: Illustrating the curse of dimensionality ............................................................................... 51 Figure 15: Schematic representation of clustering analysis ................................................................. 53 Figure 16: Schematic representation of main trajectory inference steps ............................................ 55 Figure 17: Schematic representation of gene expression comparison ................................................ 56 Figure 18: Schematic representation of a gene regulatory network reconstruction ........................... 60 Figure 19: Information dispersed throughout cells' transcriptome encodes the tumor of origin for cancer cells and the lineage subtype for normal cells. ......................................................................... 97 Figure 20: Information dispersed throughout cells' transcriptome encodes the tumor of origin for cancer cells coming from gliomas ......................................................................................................... 98 Figure 21: Cancer cells from tumor core and periphery group together. ............................................ 99 Figure 22: Information dispersed throughout cells' transcriptome encodes the tumor of origin for cancer cells coming from non-brain tumors. ...................................................................................... 100

FFABP7

  Main clinical and molecular characteristics distinguishing IDH wt from IDH mut GBM ............... 21 Table 2: Examples of single-cell high-throughput technologies ........................................................... 40

Figure 1 :

 1 Figure 1 : Primary brain and other central nervous system (CNS) tumors Figure adapted from (Ostrom et al 2020). (A) Average annual age-adjusted incidence rates with 95% confidence intervals (CI) of the top eight highest incidence cancers. Age-adjusted to the 2000 United States (US) standard population, $: female only, *: male only. Red bars correspond to primary brain and other CNS tumors. (B) Average annual age-adjusted mortality rates with 95% CI of primary brain and other CNS tumors (red bar) compared to the most common causes of cancer death (white bar) and most common non-cancer causes of death (grey bar). Ageadjusted to the 2000 US standard population, $: female only, *: male only, COPD: Chronic Obstructive Pulmonary Disease. (C) Distribution of primary brain and other CNS tumors. Percentages may not add up to 100% due to rounding.

Figure

  Figure 2 : Glioma classification based on histopathological features World Health Organization (WHO) 2007 classification. The designation of tumors as astrocytoma or oligodendroglioma is based on their histological resemblance to astrocytes or oligodendrocytes, respectively. See text for further details. GBM: Glioblastoma. Figure adapted from (Tilak et al 2021).

Figure 3 :

 3 Figure 3: Classification of diffuse gliomas based on histological and genetic features World Health Organization (WHO) 2016 classification. NOS (Not Otherwise Specified) refers to tumors for which full Isocitrate dehydrogenase (IDH) evaluation cannot be performed. Figure reproduced from (Louis et al 2016).

  four GBM subtypes (mesenchymal or ME, classical or CL, proneural or PN and neural or NE) by analyzing the global gene expression profiles of 200 GBM tissues from The Cancer Genome et al 2010). The CL subtype is

  et al 2010) (Brennan et al 2013). Proneural GBMs with IDH mutations (IDH mut -GBM) confer favorable clinical outcomes compared with IDH wt -GBM (Brennan et al 2013) (Ceccarelli et al 2016). These classifications were based on histomorphology and molecular (genetic alterations, transcriptional profiles) characteristics of GBM tissue fragments. Technological advances have allowed to analyze the transcriptomes of individual GBM cells, which revealed the coexistence of different genetic clones and transcriptomic subtypes within a given tumor. These studies further highlighted the extensive heterogeneity between GBM patient tumors but also within the same GBM (see § I.5).

Figure 4 :

 4 Figure 4: Subgroups of glioblastoma Mutations of the Isocitrate dehydrogenase 1/2 (IDH1/2) gene distinguish secondary glioblastoma (GBM) from primary GBM. Secondary GBM (first column) account for 10% of GBM, evolve from low-grade gliomas and are of better prognosis than primary GBM. Primary GBM has been classified into proneural, mesenchymal and classical subgroups based on gene expression profiles. Each subgroup is characterized by specific mutational or cytogenetic profiles. Figure adapted from (Sturm et al 2012).

  GBM are heterogeneous tumors composed of different malignant cell populations intermingled with different non-malignant cell populations (endothelial cells, immune cells, neural cells) and with the extracellular matrix (ECM) (Figure 5), creating a complex and unique tumor microenvironment (TME) in each GBM. The TME can vary across different anatomic locations, with distinct TME niches differing in terms of oxygen pressure, blood vessel density, growth factors, nutrients, metabolites and cellular composition among others (Hambardzumyan & Bergers 2015) (Schiffer et al 2018). For instance, in vivo measurements of metabolite levels in the core and infiltration zone of patient-derived xenograft GBM models revealed that distinct microenvironments have different metabolic profiles (Cudalbu et al 2021) (Mishkovsky et al 2021). The tumor core enriched in cancer cells was characterized by low levels of metabolites related to neuronal metabolism (low N-acetyl aspartate and neurotransmitters) and high levels of choline compounds, glutamine and lipids whereas the infiltration zone showed the opposite metabolic profile. The TME may also evolve over time, in response to therapeutic agents. For example, Wang and colleagues observed changes in TME cellular composition between primary and recurrent IDH wt -GBM, with a monocyte depletion and an increase in T-cell number (Wang et al 2017a). GBM and surrounding cells communicate in several ways. They use not only cell-secreted soluble factors, but also non-secretable factors including second messengers (e.g. Ca 2+ ), transcription Box 2: In vitro and in vivo models to study GBM

Figure

  al 2015)[START_REF] Valdebenito | Tunneling nanotubes, TNT, communicate glioblastoma with surrounding non-tumor astrocytes to adapt them to hypoxic and metabolic tumor conditions[END_REF] [START_REF] Pinto | Patient-derived glioblastoma stem cells transfer mitochondria through tunneling nanotubes in tumor organoids[END_REF]. The communication between GBM cells and surrounding cells in the TME is bi-directional, with surrounding cells influencing GBM cell functioning and malignant GBM cells manipulating non-malignant cells in the TME so as to sustain tumor cell growth and infiltration.

Figure 5 :

 5 Figure 5: Glioblastoma microenvironment Glioblastoma (GBM) are composed of malignant cells, as well as extracellular matrix (ECM), endothelial cells (green box), immune cells (red boxes) and neural cells (blue boxes). +: pro-tumor function, -: antitumor function, ±: mixed pro-tumor and anti-tumor functions. Figure adapted from (Broekman et al 2018).

and

  Notch signaling in GBM cells, and leukocyte chemotaxis and migration, T-cell dysfunction and cell-cell adhesion in immune cells. Pombo Antunes and colleagues found that tumor-associated macrophages in GBM have pro-angiogenic, poor antigen presenting and high T-cell suppressive capacities, hence favoring GBM cell growth (Pombo Antunes et al 2021).

  cell dissemination makes complete surgical resection of GBM impossible. Moreover, studies showed that irradiation and chemotherapy could promote the invasive phenotype of GBM cells(Lu et al 2012) (Huber et al 2013). The "spared" cells drive inevitable recurrence (de Gooijer et al 2018). Hence, better understanding of the molecular mechanisms crucial for GBM cell motility is required for developing more efficient therapies for GBM patients.

Figure 6 :

 6 Figure 6: Schematic representation of mechanisms involved in GBM cell motility. Figure reproduced from (Armento et al 2017). See text for further details.

  [START_REF] Chan | Roles of the Rac1 and Rac3 GTPases in human tumor cell invasion[END_REF] [START_REF] Tran | Increased fibroblast growth factor-inducible 14 expression levels promote glioma cell invasion via Rac1 and nuclear factor-kappaB and correlate with poor patient outcome[END_REF] [START_REF] Kwiatkowska | The small GTPase RhoG mediates glioblastoma cell invasion[END_REF] [START_REF] Gont | PREX1 integrates G proteincoupled receptor and phosphoinositide 3-kinase signaling to promote glioblastoma invasion[END_REF], Cdc42[START_REF] Fortin | Cdc42 and the guanine nucleotide exchange factors Ect2 and trio mediate Fn14-induced migration and invasion of glioblastoma cells[END_REF] [START_REF] Okura | A role for activated Cdc42 in glioblastoma multiforme invasion[END_REF] and RhoG[START_REF] Kwiatkowska | The small GTPase RhoG mediates glioblastoma cell invasion[END_REF] has been shown to enhance GBM cell migration and invasion. Expression of these genes may vary depending on the route of invasion -GBM cells invading the brain parenchyma have lower RhoA but higher Rac1 and Cdc42 activity than cells invading the perivascular area[START_REF] Hirata | In vivo fluorescence resonance energy transfer imaging reveals differential activation of Rho-family GTPases in glioblastoma cell invasion[END_REF]. In addition, downstream effectors of Rho GTPases have also been implicated in GBM cell motility (e.g. Cdc42 effector N-WASP[START_REF] Tang | c-Src and neural Wiskott-Aldrich syndrome protein (N-WASP) promote low oxygen-induced accelerated brain invasion by gliomas[END_REF]). Moreover, motile cells are characterized by a polarized cell morphology and an active machinery for actin filament. In GBM cells, it has been reported that cofilin, involved in de-polymerization and polymerization of actin filaments, enhances cell motility[START_REF] Yap | The motility of glioblastoma tumour cells is modulated by intracellular cofilin expression in a concentration-dependent manner[END_REF] [START_REF] Park | Transcriptional profiling of GBM invasion genes identifies effective inhibitors of the LIM kinase-Cofilin pathway[END_REF] [START_REF] Jin | Nogo-A inhibits the migration and invasion of human malignant glioma U87MG cells[END_REF]. During cell movement, rearranged actin cytoskeleton is connected to the ECM by formation of focal adhesion complexes (FAC). FACs are composed of paxillin and talin which are intracellularly bound to integrins. Upon ligand binding, integrins form clusters and recruit focal adhesion kinase (FAK) and vinculin(Armento et al 2017). Expression of focal adhesion proteins, like talin, have been shown to be required for GBM cell invasiveness (Sen et al 2012)(Kang et al 2015). In addition to changes in cell morphology, motile cells also have to alter their cell volume to migrate through narrow and tortuous extracellular spaces. Ion channels known to regulate intracellular volume have been shown to promote GBM cell invasion such as NKCC1 (Haas & Sontheimer 2010) (Garzon-Muvdi et al 2012) and CLCN3 (Lui et al 2010) (Song et al 2020). Several signaling pathways have been acknowledged as important regulators of GBM cell motility. Different axonal guidance molecules, which mediate cell-to-cell communication as well as cell migration and axon guidance in the brain (Roig-Puiggros et al 2020), have been shown to regulate GBM cell motility. Netrin expression enhances GBM cell invasiveness via Notch signaling activation[START_REF] Ylivinkka | Motility of glioblastoma cells is driven by netrin-1 induced gain of stemness[END_REF]. The erythropoietin-producing human hepatocellular carcinoma (Eph) receptors and ephrins also affect GBM cell motility. EphA2[START_REF] Miao | EphA2 promotes infiltrative invasion of glioma stem cells in vivo through cross-talk with Akt and regulates stem cell properties[END_REF], EphB2 (Wang et al 2012)[START_REF] Nakada | The phosphorylation of EphB2 receptor regulates migration and invasion of human glioma cells[END_REF] and EphrinB2[START_REF] Nakada | The phosphorylation of ephrin-B2 ligand promotes glioma cell migration and invasion[END_REF] [START_REF] Krusche | EphrinB2 drives perivascular invasion and proliferation of glioblastoma stem-like cells[END_REF] were found to promote GBM cell invasiveness while overactivation of several other members of the Eph/ephrin family were detected in invasive GBM cells compared to non-invasive ones[START_REF] Nakada | Ephrin-B3 ligand promotes glioma invasion through activation of Rac1[END_REF])[START_REF] Nakada | The phosphorylation of ephrin-B2 ligand promotes glioma cell migration and invasion[END_REF]. Growth factor signaling pathways have also been shown to be crucial for GBM cell motility such as the TGF-β (Joseph et al 2014)(Daubon et al 2019) and the EGFR signaling pathways[START_REF] Nakada | Receptor Tyrosine Kinases: Principles and Functions in Glioma Invasion[END_REF].Various proteins involved in the epithelial to mesenchymal transition (EMT) process, observed during normal development, also regulate GBM cell migration. For example, SNAI1/Snail, SNAI2/Slug, ZEB1 and ZEB2 exert a pro-migratory effect on GBM cell migration (Yang et al 2010)(Carro et al 2010) [START_REF] Qi | ZEB2 mediates multiple pathways regulating cell proliferation, migration, invasion, and apoptosis in glioma[END_REF] [START_REF] Myung | Snail plays an oncogenic role in glioblastoma by promoting epithelial mesenchymal transition[END_REF] [START_REF] Zhang | SHP-2-upregulated ZEB1 is important for PDGFRalpha-driven glioma epithelial-mesenchymal transition and invasion in mice and humans[END_REF].Changes in the microenvironment are also important for cell motility. The ECM has to be remodeled or destroyed to allow cell invasion. Its main components include proteoglycans, hyaluronan, link-proteins like TNC, and others. Several studies have highlighted deregulated ECM dynamics in GBM with an upregulation of fibrous proteins and laminin as well as ECM degrading and remodeling enzymes (e.g. MMP2[START_REF] Wang | The expression of matrix metalloproteinase-2 and -9 in human gliomas of different pathological grades[END_REF], MMP9[START_REF] Wang | The expression of matrix metalloproteinase-2 and -9 in human gliomas of different pathological grades[END_REF] [START_REF] Park | Role of hyaluronan in glioma invasion[END_REF] [START_REF] Zhou | Increased expression of MMP-2 and MMP-9 indicates poor prognosis in glioma recurrence[END_REF], uPA[START_REF] Chandrasekar | Downregulation of uPA inhibits migration and PI3k/Akt signaling in glioblastoma cells[END_REF] [START_REF] Brat | Pseudopalisades in glioblastoma are hypoxic, express extracellular matrix proteases, and are formed by an actively migrating cell population[END_REF] [START_REF] Kunigal | SPARC-induced migration of glioblastoma cell lines via uPA-uPAR signaling and activation of small GTPase RhoA[END_REF] (Hu et al 2015b),

  al 2016). GBM cells can also adopt different cell functioning states (functional heterogeneity).

Figure 7 :

 7 Figure 7: Main genomic alterations in GBM Main alterations in GBM involve canonical signal transduction and tumor suppressor pathways. Summary done using exome sequencing and DNA copy-number data from 251 GBM. % represents the overall alteration rate. Figure reproduced from (Brennan et al 2013).

Figure 8 :

 8 Figure 8: Models of GBM cell heterogeneityIn the clonal evolution model, cells divide and acquire mutations, generating genetic heterogeneity. Upon selective pressure, only the most fit clone(s) will survive. In the Cancer Stem Cell (CSC) model, CSCs divide asymmetrically to generate CSCs and more differentiated cells. Both display heterogeneity. Upon selective pressure, resistant CSCs will survive. In these two models, the surviving cells will then expand and cause tumor relapse. In the last model, heterogeneity occurs as a result of cellular plasticity in response to microenvironment cues. This model implies that cells adapting to selective pressure are the source of tumor relapse. Figure adapted from (Inda et al 2014).

Figure 9 :

 9 Figure 9: "Warburg effect" in proliferating tissues and cancer cells Schematic representation of glucose metabolism. Glucose is first metabolized to pyruvate via glycolysis. In the presence of oxygen, non-proliferating tissues (differentiated) completely oxidize most of that pyruvate during oxidative phosphorylation (OXPHOS). In the absence of oxygen, they redirect pyruvate away from mitochondrial OXPHOS by generating lactate (anaerobic glycolysis). Cancer cells and normal proliferative tissues tend to convert most glucose to lactate despite presence of oxygen (aerobic glycolysis). Figure reproduced from (Vander Heiden et al 2009).

  -Minh et al showed that fast-cycling cells depend more on anaerobic glycolysis, while slowcycling cells rely on OXPHOS and lipid oxidation (Hoang-Minh et al 2018). Kathagen-Buhmann et al showed that the pentose phosphate pathway (PPP) is mainly used during GBM cell proliferation, and glycolysis during migration (Kathagen-Buhmann et al 2016). These findings indicate that there exists a metabolic heterogeneity in link with cell functional heterogeneity. Altogether, these rewiring of metabolic pathways favor several of the catabolic pathways used by cancer cells for energy production while supporting anabolic pathways (synthesis of nucleotides, amino acids, and lipids) (Marin-Valencia et al 2012) in link with cell functioning such as proliferation and invasion. In addition to generating ATP, aerobic glycolysis provides glucose-derived carbon skeletons for cellular building blocks to support cancer cell growth and migration (Vander Heiden et al 2009) (Maher et al 2012) (Garcia et al 2021). The tight regulation of the activity of glycolytic enzymes including phosphofructokinase 1 (PFK1) and pyruvate kinase M2 (PKM2) in cancer cells can result in a slower glycolysis rate with an accumulation of metabolites upstream of pyruvate, and thus in their diversion into other biosynthetic pathways (Vander Heiden et al 2009). For instance, the pentose phosphate pathway (PPP) uses the glycolytic intermediate, glucose 6-phosphate (G6P) to

  Several variations in lipid metabolism have also been shown to drive changes in GBM cell functioning. Inhibition of fatty acid binding protein 7 (FABP7), a lipid chaperone mediating fatty acid uptake and subsequent oxidation, delays invasiveness of GBM models in vitro and in vivo (De Rosa et al 2012)[START_REF] Hoang-Minh | Infiltrative and drugresistant slow-cycling cells support metabolic heterogeneity in glioblastoma[END_REF]. Geng and colleagues showed that sterol O-acyltransferase (SOAT1) inhibition blocked cholesterol esterification by preventing SREBP-1-regulated fatty acid synthesis, resulting in decreased GBM cell tumorigenicity[START_REF] Geng | Inhibition of SOAT1 Suppresses Glioblastoma Growth via Blocking SREBP-1-Mediated Lipogenesis[END_REF]. Wang et al showed that inhibition of the mevalonate pathway, which controls cholesterol production and post-translational modifications of Rho-GTPases, decreased GBM cell proliferation and self-renewal(Wang et al 2017b).Polyunsaturated fatty acids (PUFA) have recently been highlighted as key lipids for GBM cell functioning. We demonstrated that variations in expression of the ELOVL2 enzyme, involved in synthesis of PUFA, controlled GBM cell tumorigenicity in patients' tumors ((Saurty-Seerunghen et al 2019), see Results section for more details). Of note, ELOVL2 has been found as one of the genes associated with superenhancers in GBM(Gimple et al 2019b).Evidences for a driver role of amino acids and their derivatives have also been provided by our team. Decreased activity of the mitochondrial enzyme SSADH, involved in GABA metabolism, was demonstrated to drive GBM cell conversion into a less aggressive functioning state. Mechanistically, decreased SSADH activity or expression results in enhanced production of the GABA by-product GHB that interfers with DNA epigenetic regulation(El-Habr et al 2017). Noteworthily, this coupling between metabolic reprograming and cell functioning was identified in the patients' tumors (El-Habr et al 2017). Variations in the levels of another metabolic intermediate linking metabolism and the epigenome, the acetyl donor for acetylation reactions acetyl-CoA, have been demonstrated to promote GBM cell migration and adhesion to ECM(Lee et al 2018b).

  -CELL RNA-SEQUENCING DATA ANALYSIS Single-cell technologies are evolving with an outstanding rapidity. Since the first published single-cell RNA sequencing (scRNA-seq) method in 2009 (Tang et al 2009), technological advances have allowed to increase throughput in terms of number of genes quantified and of number of cells analyzed while reducing cost (see § II.1). Today, scRNA-seq is the most widely used high-throughput single-cell approach. It has allowed tremendous progress in understanding intra-sample heterogeneity (e.g. diversity of immune cell populations in healthy (Jaitin et al 2014) (Szabo et al 2019) and diseased (Szabo et al 2019) (Pombo Antunes et al 2021) states), in identifying rare cell populations (e.g. cells involved in drug resistance and cancer relapse (Ho et al 2018) (Hong et al 2019)) and in studying development[START_REF] Trapnell | The dynamics and regulators of cell fate decisions are revealed by pseudotemporal ordering of single cells[END_REF] [START_REF] Petropoulos | Single-Cell RNA-Seq Reveals Lineage and X Chromosome Dynamics in Human Preimplantation Embryos[END_REF]. Despite the popularity of scRNA-seq, analysis of mRNA levels alone is not sufficient to characterize and understand biological systems. scRNA-seq methods provide an indirect readout for protein levels in a given cell -they do not capture post-translational modifications and information on protein stability and activity, hence the need for direct quantification of the proteome in single cells. High-throughput single-cell protein profiles begin to be accessible using liquid chromatography-mass spectrometry (LC-MS) based methods (Table2) (e.g. Single Cell ProtEomics by Mass Spectrometry or SCoPE-MS[START_REF] Budnik | SCoPE-MS: mass spectrometry of single mammalian cells quantifies proteome heterogeneity during cell differentiation[END_REF], improved boosting to amplify signal with isobaric labeling or iBASIL[START_REF] Tsai | An Improved Boosting to Amplify Signal with Isobaric Labeling (iBASIL) Strategy for Precise Quantitative Single-cell Proteomics[END_REF], ScoPE2[START_REF] Specht | Single-cell proteomic and transcriptomic analysis of macrophage heterogeneity using SCoPE2[END_REF]). The last eight years also saw the advent of single-cell epigenomic approaches (Table2) for profiling chromatin accessibility (e.g. scATAC-seq (Cusanovich et al 2015)[START_REF] Buenrostro | Single-cell chromatin accessibility reveals principles of regulatory variation[END_REF]), chromosome conformation (e.g. single-cell Hi-C (Nagano et al 2013)), DNAor chromatin-binding proteins (e.g. single-cell ChIP-seq[START_REF] Rotem | Single-cell ChIP-seq reveals cell subpopulations defined by chromatin state[END_REF]) and DNA modifications (e.g. single-cell methylome analysis technique based on RRBS (scRRBS)[START_REF] Guo | Single-cell methylome landscapes of mouse embryonic stem cells and early embryos analyzed using reduced representation bisulfite sequencing[END_REF]). Nowadays,

  oligonucleotide-labeled antibodies were published: Cellular Indexing of Transcriptomes and Epitopes by sequencing or CITE-seq[START_REF] Stoeckius | Simultaneous epitope and transcriptome measurement in single cells[END_REF] and RNA Expression And Protein sequencing or REAP-seq[START_REF] Peterson | Multiplexed quantification of proteins and transcripts in single cells[END_REF]. Later, in 2019, another assay, single-cell RNA and Immunodetection or RAID, was published which combines intracellular immunostaining with Antibody RNA-barcode Conjugates (ARCs) with single-cell mRNA sequencing, hence allowing joint measurement of the transcriptome and intracellular (phospho-)proteins from fixed single cells[START_REF] Gerlach | Combined quantification of intracellular (phospho-)proteins and transcriptomics from fixed single cells[END_REF].However, the throughput of the proteomics data is still much lower than the associated transcriptomics data. It is also possible to jointly capture RNA and DNA from the same cells, for example using gDNA-mRNA sequencing (DR-Seq)[START_REF] Dey | Integrated genome and transcriptome sequencing of the same cell[END_REF], G&T-seq[START_REF] Macaulay | G&T-seq: parallel sequencing of single-cell genomes and transcriptomes[END_REF] orTARGET-seq (Rodriguez-Meira et al 2019). Moreover, single-cell methods have been developed to simultaneously capture the transcriptomic and epigenomic profiles (e.g. RNA and DNA methylation using single-cell methylome and transcriptome sequencing or scMT-seq[START_REF] Hu | Simultaneous profiling of transcriptome and DNA methylome from a single cell[END_REF]). Recently, assays allowing joint measurement of chromatin accessibility and DNA or proteins have been Table 2: Examples of single-cell high-throughput technologies This table gathers the examples described in the text. The list is not exhaustive due to the constant novel development of the single-cell technologies.

  several techniques: limiting dilution, microscopy-guided isolation, flow cytometry/fluorescenceactivated cell sorting (FACS), laser capture microdissection and microfluidics (Hwang et al 2018). Dissociation and single-cell isolation should be done as quickly as possible to ensure that cell viability is not altered and to minimize changes in gene expression. Capture of all cell types present in a tissue is however not warranted. In glioma, not all normal cell types are identified -the two main ones being microglia/macrophages and oligodendrocytes although these tumors also contain astrocytes, all cells composing blood vessels, neurons, and a variety of immune cells (Venteicher et al 2017) (Filbin et al 2018) (Neftel et al 2019) (Johnson et al 2020). The potential difficulty to spontaneously recover these normal cell types may be due to their rarity and/or greater sensitivity to cell dissociation procedures (e.g. neurons) or due to inappropriate dissociation protocols which favor cancer cell separation (Tirosh & Suva 2018). It is possible to enrich for specific cell types by using immunopanning techniques for example (Darmanis et al 2017).

Figure 10 :

 10 Figure 10: Analyzing molecular profiles of cell populations using bulk and single-cell approaches Figure reproduced from (Mossner et al 2021). (A) Illustrating the coexistence of heterogeneous cell populations within a sample. (B) Analyzing multicellular bulk pieces "A" and "B" reveals differences in the global RNA expression profiles which could be wrongly attributed as intrinsic, tumor cell-specific changes. (C) In contrast, analyzing individual cells allows to associate transcriptome profiles with specific cell populations (such as tumor or normal cells) and to explore cell heterogeneity within a sample.

  Figure 11: General scRNA-seq experimental workflow Figure adapted from (Haque et al 2017). See text for details.

  al 2019), Bowtie[START_REF] Langmead | Ultrafast and memory-efficient alignment of short DNA sequences to the human genome[END_REF], BWA[START_REF] Li | Fast and accurate short read alignment with Burrows-Wheeler transform[END_REF] and STAR[START_REF] Dobin | STAR: ultrafast universal RNA-seq aligner[END_REF]. For UMI-based sequencing data, UMI sequences should be trimmed before aligning reads. The quality of aligned reads should also be assessed by evaluating the total number of reads, the proportion of uniquely mapped reads, the fraction of reads mapping to annotated exonic regions and the coverage uniformity[START_REF] Hwang | Single-cell RNA sequencing technologies and bioinformatics pipelines[END_REF]. Tools like RNA-seQC (DeLuca et al 2012) allows to do so. Cells with unusually high numbers of low-

Figure 12 :

 12 Figure 12: General analytical workflow for analyzing scRNA-seq data from complex tissues. See text for details. In red are steps applicable only when analyzing cancer scRNA-seq data.

  e. Data integration. With the growing number of scRNA-seq datasets available today, it is common practice to integrate multiple datasets collected across different conditions or time points, from different laboratories, and obtained with different sequencing techniques. This step is performed so as to have a more integrated view (e.g. when studying different development stages) or to increase number of cells under study (to increase biological signal). However, one should consider the systematic technical variations involved (e.g. differences in single-cell capture, sequencing depth, read length, plates or flow cells, protocols, sample acquisition and handling, sample composition, sampling time) that lead to batch effects in downstream analyses.

Figure 13 :

 13 Figure 13: Schematic representation of data integration strategies A. In the merging strategy, common cell groups between different datasets are identified and merged into a single group. Figure reproduced from (Haghverdi et al 2018). B. In the projection strategy, similarities between individual cells from a query sample and a reference sample are sought for. The cells from the query sample are then projected onto the most resembling cell populations of the reference sample (scmap-cell, top panel). Alternatively, cell groups from a query sample are projected onto most similar cell groups of the reference sample (scmap-cluster, bottom panel). Colors illustrate different cell populations in the reference sample. Figure adapted from (Kiselev et al 2018).

Figure 14 :

 14 Figure 14: Illustrating the curse of dimensionality As dimensions increase, the volume of the space increases, causing data to become sparse and all objects in that space to appear dissimilar to one another. Figure reproduced from (Rhys 2020).

  frequently used. In supervised approaches, an identity is assigned to each cell after learning the cell identities from annotated training data (e.g. a reference atlas)[START_REF] Abdelaal | A comparison of automatic cell identification methods for single-cell RNA sequencing data[END_REF]. While being robust to batch effects and highly reproducible, such approaches rely on accurate reference data for predicting cell identities[START_REF] Ranjan | scConsensus: combining supervised and unsupervised clustering for cell type identification in single-cell RNA sequencing data[END_REF].Examples of supervised approaches for scRNA-seq data (also known as supervised classifiers) include ACTINN based on neural networks[START_REF] Ma | ACTINN: automated identification of cell types in single cell RNA sequencing[END_REF] and SingleCellNet based on random forest algorithms[START_REF] Tan | SingleCellNet: A Computational Tool to Classify Single Cell RNA-Seq Data Across Platforms and Across Species[END_REF].Unsupervised approaches first involve grouping cells based on their transcriptomic similarity, and then annotating cell clusters, using marker genes for example, to assign labels to each cluster[START_REF] Ranjan | scConsensus: combining supervised and unsupervised clustering for cell type identification in single-cell RNA sequencing data[END_REF]. Two categories of unsupervised approaches can be used to generate cell clusters from similarity matrices: clustering algorithms and community detection methods. Clustering algorithms frequently used for scRNA-seq data include hierarchical clustering (e.g. SINCERA[START_REF] Guo | SINCERA: A Pipeline for Single-Cell RNA-Seq Profiling Analysis[END_REF], CIDR(Lin et al 2017b)) and centroid-based clustering (e.g. RaceID[START_REF] Grun | Single-cell messenger RNA sequencing reveals rare intestinal cell types[END_REF], SC3[START_REF] Kiselev | SC3: consensus clustering of singlecell RNA-seq data[END_REF]). Hierarchical clustering can be either divisive or agglomerative. Agglomerative clustering initially considers each sample as a cluster. The clusters are then progressively obtained by merging neighboring clusters together. The process continues until there is only one cluster left. This results in a dendrogram, a tree that represents the progressive clusters merging. Divisive clustering works in the opposite direction -it starts with one cluster containing all samples and then identifies subgroups. Cell clusters are obtained by selecting a number of clusters and cutting the grouping tree accordingly. Hierarchical clustering is especially useful when aiming to arrange clusters into a natural hierarchy. Because it needs to compute and store an n*n distance matrix, it can be rather slow and computationally intensive for very large datasets such as scRNA-seq datasets. Centroid-based clustering methods are suited for very large datasets. They group data by selecting cluster centroids and assigning cells to the closest cluster centroid. Each cluster center is then representative of the cluster's properties. The most popular centroid-based clustering method is the K-means algorithm (e.g. implemented in RaceID[START_REF] Grun | Single-cell messenger RNA sequencing reveals rare intestinal cell types[END_REF] and SC3[START_REF] Kiselev | SC3: consensus clustering of singlecell RNA-seq data[END_REF]. A main disadvantage of centroid-based methods is that the number of clusters expected is required as input. This information being usually unknown, one has to proceed by inference. The Hierarchical Clustering on Principal Components (HCPC) approach, which combines hierarchical and K-means clustering(Husson et al 2010), exploits the advantages of the two clustering methods. Hierarchical clustering is first performed on the data, and the resulting dendrogram analyzed. The number of clusters (N) can be determined either graphically based on the dendrogram or using indices that compare intra-cluster

  Figure 15: Schematic representation of clustering analysis Objects are clustered by similarity, i.e. by shape or color in this example.

  al 2017)(Neftel et al 2019) (Bhaduri et al 2020) [START_REF] Darmanis | A survey of human brain transcriptome diversity at the single cell level[END_REF] [START_REF] Grun | Single-cell messenger RNA sequencing reveals rare intestinal cell types[END_REF]

Figure 16 :

 16 Figure 16: Schematic representation of main trajectory inference stepsThe data are first projected into a lowdimensional space (A). A path going through the data points is then identified (B1), such that it minimizes transcriptional changes between neighboring cells. Finally, the trajectory is constructed (B2). Figure adapted from[START_REF] Street | Slingshot: cell lineage and pseudotime inference for single-cell transcriptomics[END_REF].

Figure 17 :

 17 Figure 17: Schematic representation of gene expression comparison Figure adapted from (DGEworkshop).

  specific cell functioning states for each cancer type. Secondly, cancers like GBM can be very heterogeneous tissues at the molecular level. Within the same tumor, malignant cells may use different mechanisms. For example, these mechanisms can vary depending on the cells' genomic background or according to the cells' microenvironment. In that case, different gene combinations can mark the same functioning state.II.9 Gene signature scoresCalculating gene signature scores allows to have a global view on the expression levels of a biologically-relevant set of genes. It has been used for assigning cell lineage identities (Tirosh et al 2016b) (Caruso et al 2020) and for assessing functioning states (e.g. stem-like (Patel et al 2014) (Tirosh et al 2016b), tumorigenic (Saurty-Seerunghen et al 2019) or the "activity" of biological processes or pathways (e.g. cell cycle (Tirosh et al 2016a) (Tirosh et al 2016b) (Couturier et al 2020), T-cell activity (Tirosh et al 2016a), developmental and injury response programs (Richards et al 2021)).

  the geometric mean of each signature gene for calculating a signature score (Saurty-Seerunghen et al 2019). Geometric means, instead of arithmetic means, avoid giving too much weight to zero values, which are very abundant in scRNA-seq data. In addition, this score is independent of dataset composition and different scores can be compared since gene set size is taken into account. Very recently, use of Mann-Whitney U statistic has been proposed for signature scoring (UCell tool (Andreatta & Carmona 2021b)). UCell scores have the advantage of depending only on relative gene expression in individual cells and are also unaffected by dataset composition. . They are based on measurements of gene co-expression (Figure 18) such as correlation as described by Iacono and coll. (Iacono et al 2019), mutual information (e.g. Multivariate Information-based Inductive Causation or MIIC (Verny et al 2017), Algorithm for the Reconstruction of Accurate Cellular Networks or ARACNE (Margolin et al 2006), Context Likelihood of Relatedness or CLR (Faith et al 2007), Partial Information Decomposition and Context or PIDC (Chan et al 2017)), or via regression models as used by GEne Network Inference with Ensemble of Trees or GENIE3 (Huynh-Thu et al 2010), Single-Cell rEgulatory Network Inference and Clustering or SCENIC (Aibar et al 2017) and SCODE[START_REF] Matsumoto | SCODE: an efficient regulatory network inference algorithm from single-cell RNA-Seq during differentiation[END_REF]. The idea behind all these algorithms is that two genes sharing a regulatory relationship are expected to be co-expressed. Usually, the strengths of the co-expression indexes are ranked to model gene relationships. This ranking is often used to infer the direct or indirect nature of the relationship between two genes. It should however be noted that the nature of the relationships can be modeled with more confidence when the final network is built from iterative analysis of different combinations of genes as done in MIIC[START_REF] Verny | Learning causal networks with latent variables from multivariate information in genomic data[END_REF].A recent comparison of five commonly-used GRN inference methods originally developed for bulk sample data (Partial correlation (Pcorr), Bayesian Networks (BN), GENIE3, ARACNE and CLR), and three single cell-specific methods (SCENIC, SCODE, PIDC) showed that most methods assessed perform poorly with single-cell expression data, even if they are specifically developed for this type of data[START_REF] Chen | Evaluating methods of inferring gene regulatory networks highlights their lack of performance for single cell gene expression data[END_REF], urging caution in the use and interpretation of such analyses.

Figure 18 :

 18 Figure 18: Schematic representation of a gene regulatory network reconstruction First, co-expression relationships between each possible gene pair recorded in the expression dataset are defined based on correlation, mutual information or regression models. These individual relationships can then be represented as a network, each node corresponding to a gene, and the edge thickness corresponding to the strength of the co-expression relationship. Figure reproduced from (van Dam et al 2018).

Fig. S1 (

 S1 Fig. S1 (related to Fig. 1). Removing cells with low-complexity transcriptomes and unsupervised grouping analysis of GBM and normal cells simultaneously. A-B. Filtering out cells with low-complexity transcriptomes. Cells with more than 90000 transcripts and more than 1700 genes were selected for further analyses. 1033 GBM (A) and 2417 normal cells (B) were retained. Left panels: selected cells colored in blue, rejected cells colored in gray. Right panels: cells colored by tumor. C. Unsupervised grouping analysis of the mixed set of GBM and normal cells distinguishes cancer cells from normal cells. Each dot represents a cell. Normal cells colored in light blue, cancer cells colored in violet. Left panel: PCA visualization. Right panel: tSNE visualization. Saurty-Seerunghen MS et al. Capture at the single cell level of metabolic modules distinguishing aggressive and indolent glioblastoma cells.

Fig. S7 (

 S7 Fig. S7 (related to Fig. 5) Increased ELOVL2 expression is associated with increased tumor burden. A. Down-regulated ELOVL2 expression in patient-derived GBM cells deprived of tumorigenic properties, compared to their tumorigenic counterparts. Mean±SD, Mann-Whitney test. Lee microarray dataset GEO ID: GSE4536. B. Higher ELOVL2 expression in GBM tissues compared to normal brain tissues. Mean±SD, Mann-Whitney test. TCGA tissue transcriptome dataset (microarrays) of 528 primary GBM and 10 normal brain tissues. C. Significantly higher ELOVL2 expression in GBM cells compared to normal cells. Mean±SD, Mann-Whitney test. scRNA-seq dataset from Darmanis and colleagues [14]. D-E. ELOVL2 expression is higher in GBM bearing a wild-type form of IDH1 (D) and with EGFR gene amplification (EGFR AMP ) (E) than in IDH1-mutant and non-amplified EGFR, respectively. French tissue transcriptome dataset (microarrays). IDH1 WT , n=95. IDH1 MUT n=33. EGFR WT , n=46. EGFR AMP , n=32. Mean±SD, Mann-Whitney test. F. High ELOVL2 expression is associated with a poorer survival for patients. TCGA GBM tissue transcriptomes (RNA-seq) of 153 primary GBM. Log-rank test.

  al 2010) (Snuderl et al 2011) (Sottoriva et al 2013) (Francis et al 2014) (Lee et al 2017) (Puchalski et al 2018) (Yu et al 2020). These results led us to postulate that information dispersed throughout each whole cell transcriptome encodes the tumor identity. Our observation that tumor-driven cancer cell grouping

  GBM(Neftel et al 2019) [START_REF] Yuan | Single-cell transcriptome analysis of lineage diversity in high-grade glioma[END_REF] (Pombo Antunes et al 2021), from other adult gliomas including oligodendroglioma(Tirosh et al 2016b) and diffuse astrocytoma, IDH mutant[START_REF] Venteicher | Decoupling genetics, lineages, and microenvironment in IDH-mutant gliomas by single-cell RNA-seq[END_REF], from pediatric gliomas including pediatric GBM(Neftel et al 2019) and H3K27M-glioma[START_REF] Filbin | Developmental and oncogenic programs in H3K27M gliomas dissected by single-cell RNA-seq[END_REF]), and from non-brain tumors including melanoma(Tirosh et al 2016a), head and neck cancers[START_REF] Puram | Single-Cell Transcriptomic Analysis of Primary and Metastatic Tumor Ecosystems in Head and Neck Cancer[END_REF], breast cancers[START_REF] Chung | Single-cell RNA-seq enables comprehensive tumour and immune cell profiling in primary breast cancer[END_REF] and acute myeloid leukemia (van Galen et al 2019).

Figure 19 :

 19 Figure 19: Information dispersed throughout cells' transcriptome encodes the tumor of origin for cancer cells and the lineage subtype for normal cells. Down-sampling gene numbers decreases the influence of the tumor of origin on cancer cell grouping (A) whereas it decreases the influence of the lineage subtype on normal cell grouping (B).Normalized Mutual Information (NMI) scores, relative to NMI score of the grouping analysis performed with all detected genes, are depicted. SMART-seq2 dataset from Darmanis et al 2017. 10-20 independent analyses performed with randomly-selected genes for each gene number analyzed. Mean ± SD. One-sample t-test. * p <0.05 compared to the NMI score of the grouping analysis performed with all genes detected in GBM cells. NS: non-significant.

  environments, one dominated by normal cells, the periphery, and the other by cancer cells, the tumor core. The biological relevance of this dispersed information is supported by our finding that normal cell lineage identity is also encoded by information dispersed throughout cells' transcriptome. This tumor identity does not appear to be simply linked to genomic specificities, as we previously demonstrated by considering circumscribed sets of genes reflecting the specific genomic background of each tumor (Saurty-Seerunghen et al 2019). It is also likely to integrate additional information related for example to epigenomic alterations inherited from the cell of origin. The source and the physical support of this information shared by all cancer cells from a given tumor remain to be unraveled. The founding genomic anomalies are a likely source. They may translate into a common global transcriptional background for all daughter cells of the cell of origin. As the tumor evolves, the transcriptional profile of each cell changes but cells may still retain a common ancestor transcriptional background (fingerprint) that would correspond to the tumor identity card. Another possible way by which the tumor identity card is imprinted in the cell transcriptome may be through epigenetic memory of the founding cell. Epigenetic memory corresponds to the transmission of specific gene expression patterns (selfpropagating molecular signatures) without irreversible changes of the genetic information (Bonasio et al 2010). Tumors may use this epigenetic memory to imprint and maintain their identity card in all their cells. Thus, each tumor may be characterized by different combinations of self-propagating transcriptional states which are maintained through feedback loops and networks of transcription factors (trans epigenetic signals), or by cis epigenetic signals (e.g. DNA methylation, histone modifications) that are physically associated and inherited along with the chromosome on which they act.

Figure 21 :

 21 Figure 21: Cancer cells from tumor core and periphery group together. UMAP representation. Cells colored by tumor of origin (left panel) or tumor region (right panel).

  These genes are involved in morphological rearrangements undertaken by motile cells (ACTN4 (Sen et al 2009) (Ji et al 2019) (Tentler et al 2019), PTK2 (Xia et al 2016) (Frolov et al 2016), PXN (Frolov et al 2016) (Sun et al 2017) (Lopez-Colome et al 2017), TLN1 (Kang et al 2015), VCL (Xia et al 2016)), remodeling of the matri-cellular environment (TNC (Xia et al 2016), SPARCL1 (Gagliardi et al 2017)), and pro-migratory signaling (TGF-β1 (Joseph et al 2014) (Daubon et al 2019), SMAD3 (Daubon et al 2019) , THBS1 (Daubon et al 2019)).

  al 2011) (Sottoriva et al 2013) (Lee et al 2017) (Puchalski et al 2018) (Yu et al 2020) (Neftel et al 2019) (Couturier et al 2020) (Bhaduri et al 2020) (Wang et al 2020)

  al 2019) (Darmanis et al 2017), and two datasets obtained with the 10X Genomics technique (Neftel et al 2019) (Pombo Antunes et al 2021) (Table

Fig. 1 (

 1 Fig.1 (See legend on next page)

(

  See figure on previous page) Fig.1 Capturing glioblastoma cells with high and low motile potential a. Schematic outline of the computational analytical strategy. See text for details. b. Correlations between expression of the ten genes of the motility signature. Only significant correlations are shown (p < 0.01). c. Malignant cell clustering based on the motility signature. UMAP representation. d. Different combinatorial expressions of the signature genes characterize each cluster. Heatmap representation of the relative expression of each signature element per cluster. e. Contribution of each tumor to identified clusters. Each cell cluster contains cells coming from distinct tumors (NMI score = 0.11). f. Identification of cell groups with the highest and lowest mean motility scores (C6: Mig HIGH , C2: Mig LOW ). *: Clusters with mean motility score statistically different from each of the other clusters, p <0.01, one-way ANOVA, Tukey's multiple comparisons test. NS: non-significant. Mig: Migratory/Invasive potential of Glioblastoma cells. g. Motility-related terms highlighted by ontology analysis of genes overexpressed in Mig HIGH versus Mig LOW cells. Dot plot representation. Genes overexpressed with fold change ≥2. Benjamini-Hochberg (BH)-adjusted p-value <0.05. KEGG: Kyoto Encyclopedia of Genes and Genomes, BP: Biological processes, CC: Cellular components, MF: Molecular functions. h. Enrichment in gene modules previously associated with GBM cell motility in Mig HIGH cells. EMT: genes associated with epithelio-mesenchymal transformation and GBM cell motility. oRG: genes signing for outer Radial Glia (oRG)-like malignant cell population with increased invasive behavior in GBM. TEAD: TEAD1-regulated genes involved in GBM cell motility. p <0.001, hypergeometric test. i. Linear regression model between motility score and EMT, oRG and TEAD scores. p <0.0001. j. Higher EMT, oRG and TEAD scores in Mig HIGH versus Mig LOW cells. * p <0.0001, Mann-Whitney test. b-g, and i-j: Results from N-S dataset analysis.

Fig. 2

 2 Fig.2 Metabolic characteristics of glioblastoma cells with high motile potential a. KEGG pathway analysis of metabolism genes overexpressed in Mig HIGH cells compared to Mig LOW cells highlights enrichment in metabolic pathways involved in energy production, oxidative stress response, extracellular matrix (ECM) modeling and membrane composition rearrangements. BH-adjusted pvalue <0.05. Results from N-S dataset analysis. b. Inferred trajectory of glioblastoma cells from low to high motility. Cells colored by motility clusters (top panel) and motility score (bottom panel). MPST (Mercaptopyruvate Sulfurtransferase) marks the path crossroad between low and high motile potential. c. MPST overexpression in GBM tissues compared to normal brain tissues. Expression in 163 glioblastoma (GBM) versus 207 normal brain (Normal) tissues. Boxplot representation from GEPIA2 website. TCGA RNA-seq dataset. One-way ANOVA test, *: BH-adjusted p-value <0.05. d. Schematic representation of the metabolic pathways over-mobilized in Mig HIGH cells compared to Mig LOW cells, as inferred from the KEGG pathway and trajectory modeling analyses.

Fig. 3

 3 Fig.3 Motile GBM cells exhibit enhanced ROS production and higher mitochondrial mass than their static counterparts a. Higher ROS production detected in cells migrating out of spheroids. ROS levels assessed using CellROX Green reagent. R633, P3 and 6240** PDC. Brightfield (top panel) and CellROX fluorescence signal (488nm, bottom panel) imaging. Scale bars = 50µm (R633) and 200µm (P3, 6240**). b. Higher mitochondrial mass in cells with high ROS production compared to cells with low ROS production. R633 and P3 PDC. FACS analysis of ROS levels using CellROX Deep Red reagent and mitochondrial mass using MitoTracker Green reagent. Cells classified based on their basal ROS levels. Mean ± SD, n = 3 independent biological samples, * p <0.05, unpaired t-test with Welch's correction. c. Cells with high mitochondrial mass have higher ROS production than cells with low mitochondrial mass. R633 and P3 PDC. FACS analysis of ROS levels using CellROX Deep Red reagent and mitochondrial mass using MitoTracker Green reagent. Cells classified based on their mitochondrial mass. Mean ± SD, n = 3 independent biological samples, * p <0.05, unpaired t-test with Welch's correction. d. Cells with high ROS production have higher invasive properties than cells with low ROS production. 5706**, R633 and P3 PDC. Left panel: Example of FACS-sorting of PDC into ROS LOW and ROS HIGH fractions. Right panel: Cell invasion across Matrigel-coated transwells. Mean ± SD, n = 3-4 independent biological samples, * p <0.05, unpaired t-test with Welch's correction. e. Decreasing ROS levels in GBM-PDC using 1mM NAC decreases cell invasion. 5706** and R633 PDC. Cell invasion across Matrigel-coated transwells. Mean ± SD, n = 4-7 independent biological samples, * p <0.05, unpaired t-test with Welch's correction.

Fig. 4 (

 4 Fig.4 (See legend on next page)

  Following this blueprint led to circumscribe the overall (See figure on previous page) Fig.4 MPST knockdown inhibits GBM cell migration and invasion a. Decreased MPST mRNA levels in PDC expressing shMPST. Mean ± SD, n = 5-6 independent biological samples. * p <0.05, Mann-Whitney test. b. Decreased MPST protein levels in shMPST expressing PDC. Western Blot analysis. MPST MW: 33/35 kDa. MPST levels normalized by amount of proteins loaded. Mean ± SD, n = 1-3 independent biological samples, * p <0.05, one-sample t-test. c. No major impact of MPST knockdown on cell viability. Cell viability assessed using Trypan Blue exclusion test. Mean ± SD, n = 5-12 independent biological samples, * p <0.05, Mann-Whitney test. 5706**, R633 and P3 PDC. d. MPST knockdown decreases cell migration on Matrigel. 5706**, R633 and P3 PDC. Cell migration assessed on Matrigel after 7h (5706**) or 24h (R633 and P3). Microphotographs depict examples of cell migration assays. Scale bars = 200µm. A solid and dotted line delineate the spheroid core and the migration area, respectively. The dot plot depicts the quantification of the migration. Mean ± SD, n = 10-41 independent biological samples, * p <0.05, Mann-Whitney test. e. MPST knockdown decreases cell migration into microfluidics chips. Microphotographs illustrate PDC migration in microfluidic chips loaded with equal numbers of PDC expressing a shMPST together with GFP and PDC-expressing a shControl into microfluidics chips. Note that only rare shMPST-PDC (green cells) cross through the narrow microchip channels, contrary to shControl-PDC (white cells). Scale bar = 200µm. The graph depicts the mean distances travelled by shControl-and shMPST-PDC. Mean ± SD, n = 32-40 from 4 independent biological samples, * p <0.05, Mann-Whitney test. P3 PDC, 24h migration time-length. f. MPST knockdown decreases cell invasion into collagen. 5706**, R633 and P3 PDC. Cell invasion after 24h (R633 and P3) or 40h (5706**). Microphotographs depict examples of invasion assays. Scale bars = 200µm. A solid and dotted line delineate the spheroid core and the invasion area, respectively. The dot plot depicts the quantification of the invasion. Mean ± SD, n = 4-8 independent biological samples, * p <0.05, Mann-Whitney test. g. MPST knockdown decreases cell invasion across Matrigel-coated transwells. 5706**, R633 and P3 PDC. Cell invasion assessed after 24h. Mean ± SD, n = 4-8 independent biological samples, * p <0.05, Mann-Whitney test. h. Inhibiting MPST enzymatic activity using the pharmacological inhibitor I3-MT-3 decreases cell migration. Cells treated with 200µM I3-MT-3 or vehicle for 72h (R633, P3) or 45h (GBO). Cell migration assessed 48h (R633, P3

Fig. 5

 5 Fig.5 MPST knockdown delays glioblastoma growth a. MPST knockdown delays tumor development. Bioluminescent analyses of tumor growth initiated by grafting PDC transduced with a luciferase construct and either shControl (shCTL) or shMPST constructs. 5706**, R633 and P3 PDC. n = 5 (R633) or 6 (5706**, P3) mice per group. b. MPST knockdown decreases tumor burden as shown by quantification of the tumor bioluminescent signals. 5706**, R633 and P3 PDC. DPG: days post-graft. Mean ± SD, n = 5 (R633) or 6 (5706**, P3) mice per group, * p <0.01, Mann-Whitney test. c. Kaplan-Meier survival curves demonstrating a significant survival benefit of mice grafted with PDC expressing shMPST compared to mice grafted with PDC expressing shControl. 5706** and P3 PDC. n = 6 mice per group. Log-rank (Mantel-Cox) test.

Fig.

  Fig.S1 related to Methods. Identification of malignant and normal cells in N-10X dataset A. Cell clustering based on GMM-based CNV predictions. Heatmap representation of CNV (copy number variations) predictions (A1). Potential malignancy status assigned, following cell clustering based on CNV predictions, highlighted on UMAP representation (A2). B. CNV predictions at canonical glioblastoma loci (Chr7 and 10). UMAP representation. C. Expression of marker genes of normal cell types. Pan-immune cells (PTPRC), macrophages (ITGAM, FCGR3A, CD14), microglia (CSF1R, TMEM119), T-cells (CD2, CD3D) and oligodendrocytes (MOG, MAG). UMAP representation. D. Cell clusters identified based on their repartition on UMAP plot. Hierarchical clustering followed by Kmeans clustering on UMAP components. E. Cell malignancy status assigned based on CNV status prediction, marker gene expression and clustering results. UMAP representation.

Fig.

  Fig.S2 related to Methods. Identification of malignant and normal cells in PA-10X dataset A. Cell clustering based on GMM-based CNV predictions. Heatmap representation of CNV (copy number variations) predictions (A1). Potential malignancy status assigned, following cell clustering based on CNV predictions, highlighted on UMAP representation (A2). B. CNV predictions at canonical glioblastoma loci (Chr7 and 10). UMAP representation. C. Expression of marker genes of normal cell types. Pan-immune cells (PTPRC), macrophages (ITGAM, FCGR3A, CD14), microglia (CSF1R, TMEM119), T-cells (CD2, CD3D) and oligodendrocytes (MOG, MAG). UMAP representation. D. Cell clusters identified based on their repartition on UMAP plot. Hierarchical clustering followed by Kmeans clustering on UMAP components. E. Cell malignancy status assigned based on CNV status prediction, marker gene expression and clustering results. UMAP representation.

Fig

  Fig.S3 related to Fig.1. Grouping strategy A. Schematic representation of HCPC clustering strategy integrating UMAP components. comp: components, HCPC: Hierarchical clustering on principal components, PCA: Principal component analysis, PC: Principal components, UMAP: Uniform Manifold Approximation and Projection. B. Integrating UMAP components results in more homogeneous clusters, as shown by elimination of cell outliers on the UMAP representations (compare upper left and right panels) and the increase in the Silhouette width index (compare lower left and right panels). Clustering based on expression of the ten elements of the motility signature, and performed using 4916 glioblastoma cells from N-S dataset.

Fig.

  Fig.S4 related to Fig.1. Motility signature captures cells with similar profiles in independent datasets. A-C. Clustering malignant cells from D-S (A), N-10X (B) and PA-10X (C) datasets based on the motility signature genes. From left to right subpanels: UMAP representation; identification of cell groups with the highest and lowest mean motility scores (Mig HIGH , Mig LOW ); Heatmap representation of the relative expression of each signature element per cluster, *: Clusters with mean motility score statistically different from each of the other clusters, p <0.01, one-way ANOVA, Tukey's multiple comparisons test, NS: non-significant. D. Overlaps between lists of genes overexpressed (OEG) in Mig HIGH compared to Mig LOW cell groups from distinct datasets. High overlap between lists from datasets obtained with the same sequencing techniques, with 69.2% of OEG in Mig HIGH cells from N-S dataset also overexpressed in Mig HIGH cells from D-S dataset (4852/7010), and 86% of OEG in Mig HIGH cells from N-10X dataset also overexpressed in Mig HIGH cells from PA-10X dataset (6627/7707). When comparing N-S, D-S, N-10X and PA-10X to the 3 other datasets, lists of OEG in Mig HIGH cells overlap by 56.6% (3590/6348), 53.2% (3590/6750), 48.4% (3590/7418) and 43.0% (3590/8348), respectively. Due to graphical constraints, two overlaps are not shown on the third venn diagram: 119 genes were identified as overexpressed in Mig HIGH cells in the two SMART-seq2 datasets only, and 241 overexpressed in Mig HIGH cells from N-S, D-S and N-10X datasets (not from PA-10X).

Fig.

  Fig.S5 related to Fig.1. Motility-related terms highlighted by ontology analysis of genes overexpressed in Mig HIGH versus Mig LOW cells from independent datasets Dot plot representation of enriched terms are related to motility. Analyses performed using genes overexpressed with fold change >2 (A, D-S dataset, 1486 genes) and (B, N-10X dataset, 600 genes), and fold change >1.5 (C, PA-10X dataset, 338 genes). Benjamini-Hochberg (BH)-adjusted p-value < 0.05. KEGG: Kyoto Encyclopedia of Genes and Genomes, BP: Biological processes, CC: Cellular components, MF: Molecular functions.

Fig

  Fig.S6 related to Fig.1. Enrichment in EMT (A), oRG (B) and TEAD (C) gene modules previously associated with GBM cell motility in Mig HIGH cells. EMT: genes associated with epitheliomesenchymal transformation and GBM cell motility. oRG: genes signing for outer Radial Glia (oRG)-like malignant cell population with increased invasive behavior in GBM. TEAD: TEAD1regulated genes involved in GBM cell motility. Upper subpanels: Linear regression models between motility score and EMT, oRG and TEAD scores. p <0.0001. Lower subpanels: Higher EMT, oRG and TEAD scores in Mig HIGH versus Mig LOW cells. * p <0.0001, Mann-Whitney test.

Fig.

  Fig.S7 related to Fig.3. Comparative migratory and invasive properties between GBM PDC A. Cell migration on Matrigel assessed after 24h (5706**, R633 and P3) and 4 days (6240**). The microphotographs illustrate examples for each PDC. Scale bar = 200µm. A solid and dotted line delineate the spheroid core and the migration area, respectively. The dot plot depicts the quantification of the migration. Mean ± SD, n = 6-12 independent biological samples. B. Cell invasion into collagen assessed after 40h (5706**) and 23-25h (R633, P3). Microphotographs provide examples of a cell invasion assay for each PDC. Scale bars = 200µm. A solid and dotted line delineate the spheroid core and the invasion area, respectively. The bar graph depicts the quantification of the invasion. Mean ± SD, n = 4-5 independent biological samples. C. Quantification of cell invasion across Matrigel-coated transwells after 24 hours. 5706**, R633, P3 and 6240** PDC. Mean ± SD, n = 3-5 independent biological samples.

Fig.

  Fig.S8 related to Fig.3. ROS levels measured using CellROX reagent A. FACS analysis of CellROX signal in basal conditions, and following cell treatment with the anti-oxidant N-Acetyl Cysteine (NAC, 1mM, 1h) or the ROS generator menadione (Mena, 0.1mM, 30min). 5706** PDC. Note the decrease in fluorescent signal in NAC-treated PDC and the increased fluorescent signal in Menadione-treated PDC. Mean ± SD, n = 3 independent biological samples, * p <0.05, unpaired t-test with Welch's correction. B. Higher mitochondrial mass detected in cells migrating out of spheroids or GBO. Mitochondrial mass assessed using MitoTracker Green reagent. R633 and P3 spheroids and GBOs seeded on Matrigel for 16h (P3 and GBO) or 24h (R633). Phase contrast (top panel) and MitoTracker fluorescence signal (bottom panel). Scale bars = 100µm.

Fig.

  Fig.S9 related to Fig.4 MPST knockdown using a second lentivirus, shMPST-142, also decreases cell migration and invasion Left panel: Decreased MPST mRNA levels in shMPST142-PDC compared to shControl-PDC. Q-PCR assay. Mean ± SD, n = 2-3 independent biological samples. Middle panel: No effect on cell viability assessed using Trypan Blue exclusion test. Mean ± SD, n = 3-5 independent biological samples, unpaired ttest with Welch correction. Right panel: Decreased cell migration upon MPST knockdown. Cell migration assessed on Matrigel 7h post-seeding. Mean ± SD, n = 8-9 independent biological samples, * p <0.05, Mann-Whitney test.

Table 1 : Main clinical and molecular characteristics distinguishing IDH wt from IDH mut GBM

 1

Table reproduced from (Louis et al 2016)

  different (for e.g. when analyzing DNA and RNA levels in the same cell), this step is not recommended when analyzing a single variable such as RNA levels, because it shrinks features

containing strong signals and inflates those with no signal

[START_REF] Nguyen | Ten quick tips for effective dimensionality reduction[END_REF]

. Normalizing feature variances in transcriptomic data gives the same weight to all genes and implies that gene expression variations observed between cells are not informative, but are due to technical biases rather than biological variations. In our own analyses, we did not normalize feature variances and retained these variances as potentially relevant biological information. A data transformation necessary prior analyzing transcriptomic data, and implemented by all, is sequencing depth normalization. Different samples will have different library sizes (sequencing depths), which will artificially differentiate samples. To make samples comparable to each other, samples are normalized by dividing each measurement by a corresponding sample size factor, as done when CPM (Counts Per Million) or TPM (Transcripts Per Million) values are calculated. CPM or TPM correspond to the number of gene-mapped reads (Counts or Transcripts) divided by the total number of mapped reads per cell, the resulting number being then divided by one million (Per Million). This sequencing depth normalization is required for between-sample comparisons (i.e. when comparing gene expression between cells). TPM adds another level of transformation prior to sequencing depth normalization.

(

  Tirosh et al 2016b) [START_REF] Venteicher | Decoupling genetics, lineages, and microenvironment in IDH-mutant gliomas by single-cell RNA-seq[END_REF] [START_REF] Filbin | Developmental and oncogenic programs in H3K27M gliomas dissected by single-cell RNA-seq[END_REF]. As expected, analyses selecting these features disclose first cell characteristics governed by marked differences in gene expression, i.e. lineage identities or similarities and genomic anomalies. Other cell characteristics depending on more subtle and progressive changes in gene expression, such as cell functioning states, are not highlighted. Our analyses seeking for cell functioning states, we therefore chose not to filter genes based on these criteria, and to keep all information contained in the datasets.

	There are four main groups of DR methods: linear (e.g. principal component analysis or PCA,
	multidimensional scaling or MDS), non-linear (e.g. t-Stochastic Neighbour Embedding or t-SNE,
	Uniform manifold approximation and projection or UMAP), neural network-based (e.g. variational
	autoencoder or VAE) and model-based (e.g. Zero Inflated Factor Analysis or ZIFA), the most widely
	used methods being cited into parenthesis. In a comparative study between 18 different DR

methods, Sun and colleagues showed that PCA and UMAP performed overall better in single-cell data analysis than the abovementioned methods in terms of their accuracy and effectiveness in downstream analyses, their robustness and stability, and their performance in large-scale data applications

[START_REF] Sun | Accuracy, robustness and scalability of dimensionality reduction methods for single-cell RNA-seq analysis[END_REF]

. In PCA, new variables corresponding to linear combinations of the initial variables are constructed (principal components or PCs). Only the PCs capturing the maximum variance of the data are retained for further analysis

[START_REF] Lever | Principal component analysis[END_REF]

. Being a linear DR method, it preserves global data structure very well

[START_REF] Nguyen | Ten quick tips for effective dimensionality reduction[END_REF]

. However, PCA doesn't perform that well in highly heterogeneous data

[START_REF] Xiang | A Comparison for Dimensionality Reduction Methods of Single-Cell RNA-seq Data[END_REF]

. On the other hand, UMAP represents the best approximation of underlying topology of the data. It presents the following advantages: (1) it yields meaningful representations, (2) it is able to resolve subtly differing cell populations, (3) it preserves both local and global structure, (4) it performs quite well in highly heterogeneous data, (5) it is quite rapid, and (6) it can be used for large number of cells

[START_REF] Becht | Dimensionality reduction for visualizing single-cell data using UMAP[END_REF]

[START_REF] Xiang | A Comparison for Dimensionality Reduction Methods of Single-Cell RNA-seq Data[END_REF]

.

Acknowledgements

LIST OF FIGURES AND TABLES viii

List of Figures

RESULTS

GENERAL DISCUSSION

Appendix 2: Glossary

Cell functioning/functional state: The condition of a given cell which translates its physiological function in the whole tissue (e.g. migratory or static, tumorigenic (able to initiate tumors) or nontumorigenic)

Cell plasticity: Cell's ability to reversibly change its phenotype in response to environmental cues, without genetic mutations [START_REF] Shen | Cell plasticity in cancer cell populations[END_REF] Cell type: Cell's identity based on its differentiation level (e.g. stem cells, OPC and oligodendrocytes) Dropout An event in which a transcript is not present in the sequencing data due to inefficient mRNA capture, stochastic mRNA expression and/or amplification failure [START_REF] Haque | A practical guide to single-cell RNA-sequencing for biomedical research and clinical applications[END_REF] Invasion: Cell's ability to remodel its microenvironment, to penetrate and move into its surroundings (de Gooijer et al 2018)

Migration: Cell's ability to move around in a space that is freely available (de Gooijer et al 2018)

Pseudo