Integration of single-cell RNA sequencing and spatial transcriptomics to reveal the glioblastoma heterogeneity

Glioblastoma (GBM), a deadly brain tumor, is still one of a few lasting challenges of contemporary oncology. Current therapies fail to significantly improve patient survival due to GBM tremendous genetic, transcriptomic, immunological, and sex-dependent heterogeneity. Over the years, clinical differences between males and females were characterized. For instance, higher incidence of GBM in males or distinct responses to cancer chemotherapy and immunotherapy between males and females have been noted. Despite the introduction of single-cell RNA sequencing and spatial transcriptomics, these differences were not further investigated as studies were focused only on revealing the general picture of GBM heterogeneity. Hence, in this mini-review, we summarized the current state of knowledge on GBM heterogeneity revealed by single-cell RNA sequencing and spatial transcriptomics with regard to genetics, immunology, and sex-dependent differences. Additionally, we highlighted future research directions which would fill the gap of knowledge on the impact of patient’s sex on the disease outcome.


Introduction
Glioblastoma (GBM) is one of the deadliest human tumors, with a 14-month median survival length and five-year overall survival (OS) of approximately 6.8% (https://seer.cancer.gov/). Despite intensive research and the introduction of novel therapeutic regimens, the survival rate has not been improved in the last decades. 1 The disappointing results are mainly due to ineffective surgical resection and rapid local progression. 2 Moreover, there are no useful biomarkers to detect the emergence of GBM, and the early course of the disease is often asymptomatic. 3 Nowadays, GBM is treated with surgery, temozolomide-based (TMZ) chemotherapy, and radiotherapy. 2 The failure of conventional and targeted therapies is most likely due to intratumoral heterogeneity, intrinsic mechanisms of cell death resistance (due to a high frequency of TP53 and PTEN mutations) and redundant prosurvival signaling pathways. 4 Additionally, in accordance with recent reports, patient sex may have a major impact on GBM therapeutic outcomes and prognosis. The response rate to conventional therapies is higher in females, whereas immunotherapy works better in males due to higher molecular and cellular heterogeneity of glioma cells. 5 Moreover, a mutation in the isocitrate dehydrogenase 1 encoding gene, IDH1, contributes to better chemotherapy outcomes and prolonged OS in males only. 6-8 On the other hand, hypermethylation of the promoter of the MGMT gene coding for O6-methylguanine-DNA methyltransferase enhances the effect of TMZ chemotherapy and prolongs OS in females. 8,9 Also, there are several negative prognostic biomarkers which are listed in the Figure 1. 8,10-13 Figure 1. Major sex-dependent clinicopathological and gene expression differences affecting therapeutic outcomes, prognosis, genetics, and immunology of glioblastoma patients. Legend: Btg2 -B cell translocation gene 2; CDKcyclin-dependent kinase; FZD7frizzled class receptor; Hmg2high mobility group box 2; IDH1isocitrate dehydrogenase 1; ILinterleukin; MGMT -O6-methylguanine-DNA methyltransferase; MHCmajor histocompatibility complex; p63transformation-related protein 63; RB1retinoblastoma protein 1; Shhsonic hedgehog human; TNF-αtumor necrosis factor alpha.

REVISED Amendments from Version 1
We have corrected Figure 1, and the legend of Figure 3 according to the reviewers' comments. Table 1 was also removed.
Any further responses from the reviewers can be found at the end of the article Over ten years ago, the first GBM transcriptional subtypes were identified with a partial enrichment of PDGFRA or EGFR alterations. 14,15 However, with the recent introduction of single-cell RNA sequencing (scRNA-seq) and spatial transcriptomics (ST), new cellular states were revealed. Interestingly, distinct states may co-exist within the tumor with variate frequency and are associated with genetic alterations. 4,16 Moreover, GBM heterogeneity manifests in unique developmental states of GBM cells in the tumor. GBM mimics mechanisms of neural development, thus it contains GBM stem cells (GSCs) which are thought to play a major role in possessing tumor-propagating potential and exhibiting preferential resistance to radiotherapy and chemotherapy. 17,18 Finally, GBM has a heterogeneous and highly immunosuppressive TME composed of normal brain residents such as neurons, astrocytes, oligodendrocytes, and microglia, immune system infiltrating cells including mostly monocytes/macrophages, as well as endothelial and mesenchymal cells. 19,20 Glioma-associated microglia and macrophages (GAM) are the most abundant population of immune cells, constituting up to 30% of the tumor mass. 21-23 Despite the high content of glioma-associated microglia and macrophages (GAM), GBM is considered immunologically "cold" due to relatively poor infiltration of activated T cells and anergic state of those present in TME. 24 Moreover, the broad range of cell-to-cell interactions between cancer cells and components of the TME affects the biological status of the tumor with an increase in its evasion capacity and resistance to treatment. 25 In summary, a better understanding of multilayer GBM heterogeneity, including genetics, epigenetics, developmental stages, TME, and immunology, is required to establish effective therapies. 16 Thus, the purpose of this mini-review was to discuss recent advances in exposing the GBM heterogeneity and highlight future research directions.
Integration of single-cell RNA sequencing and spatial transcriptomics Although scRNA-seq should expose a detailed functional characterization of transcriptomics of a single cell and allows integration of in vivo states with in vitro models, in principle, it provides an indirect inference of cellular interactions. 26 However, with the recent introduction of spatial transcriptomics (ST), it is believed that spatial and functional organization are strictly related, especially in the context of a neuronal tissue. 27 ST is a revolutionary method that enables the characterization of cellular interactions and spatial organization of examined tissues. 28 Given the fact that ST is not yet at single-cell resolution, methods for the integration of scRNA-seq and ST are vital to understand the heterogeneity of GBM ( Figure 2). With the rapid development of complex algorithms and machine learning technologies, a variety of tools to integrate scRNA-seq and ST were recently developed, which enabled a precise characterization of GBM and its TME. 29-36 Figure 2. Integration (deconvolution) strategy of single-cell RNA sequencing and spatial transcriptomics data (based on the 10x Genomics protocol). Single-cell RNA sequencing provides high-throughput and high-resolution profiling of gene expression. However, it lacks spatial information due to tissue dissociation. Conversely, ST offers a spatial context without single-cell resolution. Currently, the spot diameter of ST and Visium platforms is 100 μm and 55 μm, respectively, capturing from one to 30 cells. Thus, to gain a spatial single-cell resolution, it is necessary to integrate both methods.
The first layer of glioblastoma heterogeneity: bulk-RNA sequencing Over the years, with the constant introduction of more advanced sequencing methods and the launching of The Cancer Genome Atlas (TCGA) in 2008, a few layers of transcriptional GBM heterogeneity were discovered ( Figure 3) providing a classification of GBM subtypes. 14,16,26 In 2010, Verhaak et al., using bulk-RNA sequencing data of 200 GBM patients from TCGA, identified four distinct transcriptional subtypes of GBM: 1) classical, 2) mesenchymal, 3) proneural, and 4) neural. The results were further validated in a combined cohort of 260 GBM patients from previous studies and confirmed the association of each subtype with certain genetic events. 37-40 The classical subtype was highly associated with chromosome 7 amplification paired with chromosome 10 loss, EGFR gene alterations, and disruption of RB, Notch, and Sonic hedgehog signaling. The mesenchymal subtype was characterized by alterations in NF1 and PTEN genes affecting the AKT pathway and expression of mesenchymal markers and tumor necrosis factor superfamily and NF-κB proteins. The proneural subtype was associated with various gene alterations, including PDGFRA, IDH1, TP53, PIK3CA/PIK3R1, and higher expression of OLIG2, SOX, DCX, DLL3, ASCL1, and TCF4. Finally, the neural subtype was associated with the expression of neuron markers NEFL, GABRA1, SYT1, and SLC12A5. However, it was later confirmed that it is not a tumor-specific subtype of GBM with a lack of gene abnormalities but resected fragments with a high contribution of normal tissues. 15 Next, identified subtypes were assigned to neural cell types using transcriptomic data gene sets. 41 The classical, mesenchymal, and proneural subtypes were associated with murine astrocytic, astroglia, and oligodendrocytic signatures, respectively. Interestingly, the frequency of each subtype can vary within the same tumor as multiple subtypes can co-exist or change over time and as a response to therapy. 4,42 Interestingly, the prevalence of mesenchymal, neural, and proneural subtypes is higher in males than females, while the classical subtype occurs with the same frequency. 10,43 The second layer of heterogeneity: single-cell RNA sequencing In 2019, Neftel et al. showed the relationships between genetic subtypes and cellular states by deconvolution of scRNAseq and TCGA bulk data on GBMs with lineage tracing in GBM murine models. 16 The researchers depicted four cellular states: 1) neural progenitor-like (NPC-like), 2) oligodendrocyte-progenitor-like (OPC-like), 3) astrocyte-like (AC-like), and 4) mesenchymal-like (MES-like) corresponding to previously established TCGA signatures ( Figure 3). 14 As previously mentioned, these states may co-exist within the same tumor with different frequencies influenced by genetic alterations in CDK4, PDGFRA, EGFR, and NF1, which favor a particular state, respectively. Based on distinct gene expression patterns signatures can be divided into mesenchymal (MES1-like [hypoxia-independent], MES2like [hypoxia-dependent]) and neuro-developmental (NPC1-like, NPC2-like, OPC-like, AC-like) states. 44,45 In general, GBM cells correspond primarily to one of the four states, however, each of the tumors contains at least two cellular states, with most tumors containing all four states. The most frequent hybrid states are AC-like/MES-like, NPC-like/OPC-like, and AC-like/OPC-like.
On the other side, scRNA-seq provided novel insights into GBM immunology (Figure 4), especially on the localization of GAMs within the tumor. Microglia tend to reside in the tumor periphery with the adjacent brain parenchyma, while tissue-invading monocyte-derived macrophages (MDM) are most abundant within the tumor core. [46][47][48] Moreover, the expression of immune checkpoint receptor ligands differs in myeloid cells between tumor core and peritumoral tissue. 46 Comparison of microglial cells from human IDH wild-type GBM and age-matched controls revealed substantially downregulated expression of microglia core genes and upregulated expression of inflammatory-(IFI27, IFITM3), metabolic-(LPL, APOE, TREM2), and hypoxia-associated (HIF1A, VEGFA) genes. 49 In terms of sex-dependent differences, Ochocka et al. showed that microglia MHC class II-associated genes were significantly upregulated and more reactive in males than in females. 50 Furthermore, Pombo et al. investigated the evolution of functional GAM profiles across disease stages by sequencing samples from newly diagnosed and recurrent GBMs. Microglia-derived GAMs were predominant in newly diagnosed tumors, but were surpassed by more heterogeneous MDMs in the recurrent ones, especially in the hypoxic tumor environment. 51 Regarding T cells, Mathewson et al. showed that the inhibitory natural killer (NK) cells receptor CD161 is expressed in tumor-infiltrating lymphocytes, but absent in T regulatory cells (Tregs) or patient-matched peripheral blood mononuclear cells (PBMCs). Moreover, CLEC2D (CD161 ligand) was primarily expressed by malignant and myeloid cells, revealing similarities with the PD-1/PDL-1 (programmed death-1/ programmed death-1 ligand) system. The third layer of heterogeneity: spatial transcriptomics Local cellular interactions between tumor and cells located in TME play a major role in the adaptation of GBM and facilitate growth, infiltration, and therapy resistance, contributing to unique spatial signatures in GBM. [52][53][54] In 2022, Ravi et al. published an atlas of spatially resolved transcriptomics of 28 specimens (20 patients) and complemented it with spatially resolved metabolomics and proteomics. 26 The researchers described five spatially distinct transcriptional programs of GBM: 1) radial glia, 2) reactive-immune, 3) neural development, 4) spatial OPC, and 5) reactive-hypoxia ( Figure 3). The first two were associated with high expression of astrocyte -related genes (GFAP, AQP4, VIM, CD44). Specifically, radial glia program had an increased expression of radial-glia-associated genes (HOPX, PTPRZ1) and reactive-immune program had a functional enrichment of inflammation-associated genes (HLA-DRA, C3, CCL4, CCL3) and interferon-γsignaling. The next two programs were associated with neuronal lineages (neurons or oligodendrocytes) and were named accordingly. The last program was associated with hypoxia-response (VEGFR, HMOX1, GAPDH) and glycolytic (LDHA, PGK1) genes. In order to integrate novel programs with already established bulk and single-cell classifications, spatial-weighted regression and bilateral integration of the top-scoring gene signatures were carried out, confirming an overlap between radial glia, spatial OPC, neuronal, reactive-hypoxia, and AC-, OPC-, NPC-, and MES2like (hypoxia-dependent) states, respectively.

Reactive-immune program
The hybrid meta-modules described by Neftel et al. were associated with reactive-immune program suggesting a close functional relationship between AC-and MES-like states. 16 In further analysis, an imaging mass cytometry-based singlecell profiling showed significant enrichment of myeloid and lymphoid cells across hybrid regions. A substantial enrichment of tumor-associated myeloid cells and T cells was found among the reactive immune program. Moreover, the mean PD-1 protein level on T cells was increased, suggesting locally enhanced immunosuppression. These findings were supported by the enrichment of CD163+ myeloid cells, which support phagocytosis and immunosuppression. Furthermore, a previous study by Ravi et al. in 2022 focused on T-cell dysfunction and indicated that exhausted T cells are preferentially located within regions with mesenchymal transcriptional programs. The study revealed that the spatial and functional interaction between the myeloid and lymphoid compartment leads to an interleukin-10 mediated T cell exhaustion. 55 This was further confirmed that the enrichment of memory and exhausted T cells occur in the reactiveimmune and reactive-hypoxia areas. 26

Reactive hypoxia program
The reactive hypoxia program is associated with histologically confirmed areas of necrosis and a high prevalence of copy-number alterations (CNAs), including focal amplification of oncogenes or losses of tumor suppressors. Spatially resolved metabolomics of reactive-hypoxia regions revealed a significant enrichment of the pentose phosphate pathway, phosphoadenylate metabolism, glycolysis, and amino sugar metabolism. 26 Ravi et al. showed that hypoxia and oxidative stress highly contribute to genomic instabilities, including various chromosomal alterations, that are driving forces in GBM resistance to therapies. Moreover, the reactive-hypoxia program is prevalent among non-cycling cells resulting in a S-phase arrest contributing to the genomic instability. 56 Finally, based on migratory gene-expression signatures, the effect of oxidative stress on cellular migration was explored, revealing the opposing drivers of genomic diversity, resulting in clonal evolution in GBM.

Conclusions
In this mini-review, we showed a comprehensive overview of the GBM heterogeneity which was revealed by scRNA-seq and ST, and the urgent need to integrate both methods in future research. We discussed distinct layers of heterogeneity with regard to GBM genetics, transcriptomics, immunology, and patients' sex. Worryingly, despite a growing number of clinical reports about the impact of patients' sex on GBM prevalence, therapeutic outcomes, and prognosis, there are a few studies considering patients' sex while using scRNA-seq or ST. Hence, in the future, it is crucial to distinguish patients' sex in order to characterize potential differences which could have a major impact on the development of therapeutic agents and overcoming GBM treatment resistance.

Data availability
There are no underlying data associated with this article.

Open Peer Review
transcriptomics, immunogenomics, and sex. In the form, there are a couple of typos and grammar mistakes as well as redundant "words" that should be corrected like "on the other side" on pages 3 and 6 and the repetition of "on the other hand" in different contexts.

○
The authors should also correct the label "in this study" on page 4, as this is a mini-review.

○
For the content, all is well-covered but it would be better if sex is integrated into the other 3 pillars of heterogeneity.

○
The mini-review can be Approved after these minor corrections.
Is the topic of the review discussed comprehensively in the context of the current literature? Yes

Is the review written in accessible language? Partly
Are the conclusions drawn appropriate in the context of the current research literature? Yes are only a few studies investigating sex-dependent differences. Studies on which the heterogeneity classifications were made did not take into account this matter. Hence, the purpose of this mini-review was to point out and address this gap in the field.

Theo Mantamadiotis
The University of Melbourne, Parkville, Australia The mini-review explores the very interesting observations that male patient with GBM exhibit a higher incidence and fair worse than female patients with GBM. The authors discuss many valid biological phenomena related to the tumor microenvironment and sex-specific differences among tumor infiltrating immune cells and other cell types found in GBM.
Although the sex-specific biology is discussed, the review is divided into sections based on technologies, rather than biology. This approach distracts from the interesting biology. Table 1 seems more appropriate for a technical or bioinformatic focused review, unless the authors can demonstrate that each of these approaches has discovered something relevant to the biology of sex-specific differences.
Other parts of the manuscript which require improvement are Figure 1 legend text. The diagram and legend could be better presented. A clear explanation of what Figure 1 shows is required. The manuscript needs editing to correct language errors. For example, in the Abstract, the words "However, despite.." are used together. "However" is redundant here, so should be deleted. The word "exposed" as used in the abstract and with the body text, is wrong -the authors mean "revealed". Another example is in the Conclusions section. Near the end of Conclusions, "it is warrant" is used but it's unclear what the context is. Reviewer Expertise: Cancer biology, brain cancer, molecular biology, gene expression, animal models, developmental biology I confirm that I have read this submission and believe that I have an appropriate level of expertise to confirm that it is of an acceptable scientific standard, however I have significant reservations, as outlined above.
Author Response 22 Apr 2023 Jakub Mieczkowski Table 1 seems more appropriate for a technical or bioinformatic focused review, unless the authors can demonstrate that each of these approaches has discovered something relevant to the biology of sex-specific differences. Table 1 was removed from the main text. There were no sex-specific differences mentioned in these articles.
Other parts of the manuscript which require improvement are Figure   The title and description of the figure were changed accordingly. Names of specific subtypes were cited from primary publications. The goal of this article was to show