Tumor Heterogeneity in Glioblastomas: From Light Microscopy to Molecular Pathology

Simple Summary Glioblastomas (GBMs) are the most frequent and aggressive malignant tumors arising in the human brain. One of the main reasons for GBM aggressiveness is its diverse cellular composition, comprised by differentiated tumor cells, tumor stem cells, cells from the blood vessels, and inflammatory cells, which simultaneously affect multiple cellular functions involved in cancer development. “Tumor Heterogeneity” usually encompasses both inter-tumor heterogeneity, differences observed at population level; and intra-tumor heterogeneity, differences among cells within individual tumors, which directly affect outcomes and response to treatment. In this review, we briefly describe the evolution of GBM classification yielded from inter-tumor heterogeneity studies and discuss how the technological development allows for the characterization of intra-tumor heterogeneity, beginning with differences based on histopathological features of GBM until the molecular alterations in DNA, RNA, and proteins observed at individual cells. Abstract One of the main reasons for the aggressive behavior of glioblastoma (GBM) is its intrinsic intra-tumor heterogeneity, characterized by the presence of clonal and subclonal differentiated tumor cell populations, glioma stem cells, and components of the tumor microenvironment, which affect multiple hallmark cellular functions in cancer. “Tumor Heterogeneity” usually encompasses both inter-tumor heterogeneity (population-level differences); and intra-tumor heterogeneity (differences within individual tumors). Tumor heterogeneity may be assessed in a single time point (spatial heterogeneity) or along the clinical evolution of GBM (longitudinal heterogeneity). Molecular methods may detect clonal and subclonal alterations to describe tumor evolution, even when samples from multiple areas are collected in the same time point (spatial-temporal heterogeneity). In GBM, although the inter-tumor mutational landscape is relatively homogeneous, intra-tumor heterogeneity is a striking feature of this tumor. In this review, we will address briefly the inter-tumor heterogeneity of the CNS tumors that yielded the current glioma classification. Next, we will take a deeper dive in the intra-tumor heterogeneity of GBMs, which directly affects prognosis and response to treatment. Our approach aims to follow technological developments, allowing for characterization of intra-tumor heterogeneity, beginning with differences on histomorphology of GBM and ending with molecular alterations observed at single-cell level.


Introduction
Gliomas are the most frequent primary malignant brain tumors, with glioblastoma (GBM) alone accounting for 15% of all central nervous system (CNS) tumors [1]. The optimum treatment with maximum safe surgical resection, followed by radiation with concurrent and adjuvant temozolomide (TMZ), significantly increases the overall survival (OS) of GBM patients compared to radiation monotherapy [2][3][4]. More recently, the addition were included in a database containing information about the journals, type of heterogeneity, molecular methods used, main field of knowledge, and type of tumors included in the study. Although this is not intended to be a systematic review, we thought this analysis would help the readers have an overview of the evolution of the field in tumor heterogeneity in GBM.
Two authors (A.P.B. and B.E.S.) individually filtered abstracts deemed more pertinent to the theme of this review, and there was consensus on the inclusion of 70 articles. After careful reading of the articles, we excluded the ones that were not within the scope of the review. In course of preparing the current manuscript, we added references that were not identified or had not been published during the initial search process. Statistics and plot designs were performed with SPSS 25 (IBM) and R statistics [10].
Ninety papers were deemed the most relevant for the theme, while the remaining references were indirectly associated with the theme of GBM tumor heterogeneity. The majority of the articles specifically addressed GBMs (51 out of 90 articles). However, articles including all diffuse gliomas (WHO grades 2-4), high-grade gliomas (HGGs-WHO grades 3 and 4), and lower grade gliomas (LGGs-WHO grades 2 and 3) were included because histopathological grade alone has decreasing importance as determinant for GBM diagnosis [11]. Some of the selected papers addressed pediatric HGGs, which are entities molecularly different from adult GBMs [11]; however, those references were included due to their importance in the field of neuro-oncology. Figure 1 depicts the overall profile of the papers. Remarkably, the exponential advances in molecular pathology over the last decade have resulted in multiple comprehensive papers, encompassing complex multi-platform analyses that are more likely to be accepted by premier journals (Figure 2). This may be explained both by the technological development and by the mounting interest in exploring tumor heterogeneity in GBMs to support advances in image exams and therapeutic Remarkably, the exponential advances in molecular pathology over the last decade have resulted in multiple comprehensive papers, encompassing complex multi-platform analyses that are more likely to be accepted by premier journals (Figure 2). This may be explained both by the technological development and by the mounting interest in exploring tumor heterogeneity in GBMs to support advances in image exams and therapeutic approaches, discovery/validation of therapeutic targets and of prognostic/predictive biomarkers. Bioinformatics tools and algorithms have supported the analysis of the enormous amount of data produced by multiplatform high-throughput profiling studies, but are beyond our scope in this review. Timeline of methods by year of publication based on the articles selected in this review. "Histopathology" = H&E and IHC stains (the 1st WHO classification is from 1979). "Multiomics" = studies with independent results from two or more high-throughput methods; studies with multiple methods for determination of groups in the study were not included in this count.

Inter-Tumor Heterogeneity: Classification and Grading of Gliomas
The glioma classification arose from neurosurgeons' need to predict survival and propose therapeutic approaches for patients with brain tumors. In the 1920s, Bailey and Cushing classified brain tumors according to the embryological origin as astrocytomas (the most numerous group), ganglioneuroma, oligodendroglioma, ependymoma, pinealoma, and papilloma choroideum. GBM was identified as a less differentiated astrocytoma, then named "spongioblastoma multiforme" [12].
Multiple histopathology-based grading systems sought to refine the prognosis of astrocytomas. The first World Health Organization (WHO) grading system in 1979 consisted of three histopathological grades, and did not include GBM as a diagnosis [13]. Later, the Saint Anne-Mayo system increased reproducibility of astrocytoma grading, by systematically evaluating four histopathological features: nuclear atypia, mitotic count, endothelial proliferation, and necrosis ( Figure 3) [14].
The omics revolution of the past two decades shifted glioma classification from subjective histopathological criteria toward molecular profiling based on the genome, transcriptome, and methylome of individual tumors [15][16][17][18]. In this setting, as studies of brain tumors have always been at the forefront of biomedical research, GBM was one the first tumors to have its genome and transcriptome described by The Cancer Genome Atlas (TCGA) initiative [19,20]. Consequently, brain tumors were the first solid tumors to incorporate molecular features in their nomenclature. The current WHO 2016 classification still uses the histopathological grading system, but now molecular markers such as isocitrate dehydrogenase (IDH) 1/2 mutations and co-deletions in chromosomes 1p and 19q complete the integrated diagnosis of gliomas [21][22][23]. This resulted in three main glioma classes: IDH mutant, 1p19q codeleted (oligodendrogliomas), IDH mutant, 1p19q intact (astrocytomas), and IDH wild-type gliomas. Specific alterations are associated with each group, such as mutations on CIC, FUBP1, and TERT promoter; on ATRX and TP53; and on EGFR, PDGFRA, and TERT promoter, respectively [21,23]. In 2020, the consensus C-IMPACT-NOW 6 proposed an updated nomenclature for astrocytomas (oligodendrogliomas were not included) [11,24,25]. For IDH mutant astrocytomas, any combination of microvascular proliferation, necrosis (characteristics of anaplastic astrocytomas) or CDKN2A/B homozygous deletion now allows for the diagnosis of Astrocytoma, IDH-mutant, WHO grade 4-note the use of Arabic is now suggested, rather than Roman, numbers [11]. The diagnosis of GBM is now exclusive for IDH wild-type tumors with specific molecular features (gains in chromosome 7 and losses in chromosome 10, EGFR amplifications, and TERT pro- moter mutations), independent of the histopathological grade [11,26]. These changes are recommendations for the upcoming WHO classification predicted to be released in 2021.    Several other classifications based on the molecular profile of GBMs were proposed even before the WHO 2016 classification was released. Gene expression signatures categorize GBMs as proneural (enriched in IDH mutant tumors), neural, classical and mesenchymal, each with different prognosis, and histogenesis linked to, respectively, oligodendrocytes, neurons, astrocytes, and astroglial cells [18]; however, the possibility of the neural subtype corresponding to normal entrapped neurons was raised by subsequent authors [19,27]. DNA methylation pattern identified a group of GBMs and LGGs with hypermethylation at a large number of loci, called CpG Island Methylator Phenotype (G-CIMP), and favorable outcomes, which are virtually always related to the presence of IDH mutations [16]. Finally, a comprehensive analysis of these classifications defined three subgroups of IDH mutant gliomas (G-CIMP high-LGm1, G-CIMP low-LGm 2, and 1p19q codeleted-LGm 3) and four subgroups of IDH wild-type gliomas (classicallike-LGm4, mesenchymal-like-LGm5, and two rare subtypes-LGm6-GBM and PA-like), closely associated with transcriptional subtypes [15].
Not surprisingly, the molecular subtypes are closely associated with the main prognostic factors in GBM, IDH mutations and methylation of the promoter region of O 6methylguanine-DNA methyltransferase (MGMT) [21]. Virtually all proneural GBMs harbor IDH mutations and usually present G-CIMP, therefore reducing the prognostic/predictive power of MGMT promoter methylation in this group of tumors. Furthermore, MGMT promoter methylation is considered a predictive marker only in classical GBM, but not in the mesenchymal subtype [19].
Because of the possibility of establishing GBM prognosis by other methods, despite its importance, the transcriptomic classification of GBMs has been used more in the research field than in the clinical setting. Recently, many authors have tried to develop and validate GBM sub-classifications using histopathological features and immunohistochemistry (IHC) approaches [27][28][29]. Initially, only moderate correlation of molecular features with the histopathological grades was observed [15]. However, in further analyses, small cell morphology and microvascular proliferation ( Figure 3c) were associated with the classical subtype; oligodendroglial features with proneural subtype; pleomorphism and epithelioid cells, as well as inflammatory infiltration (Figure 4b,c and Figure 6a-c), with the mesenchymal subtype [27]. Moreover, a strong correlation between the GBM transcriptional subtypes and IHC profile was demonstrated, with higher expression of ASCL1, OLIG2, and PDGFRα in proneural GBMs; EGFR in classical GBM; and p53, pN-DRG1, YKL40, and MET in mesenchymal GBM, which were able to predict the molecular subtypes with high efficacy using machine learning tools [27]. Other authors characterized a classical/proliferating profile (overexpression of EGFR and Olig2, and high proliferative activity), and a mesenchymal/microglial profile (overexpression of ALDH1A3 in tumor cells and Iba-1 in microglia) by IHC, which were closely related to the transcriptional profiles with similar names [29]. A subset of GBMs contained areas of both profiles ("subtype heterogeneous") and had worse OS. Their results also confirmed the inflammatory profile of the mesenchymal subtype, which predominates in recurrent GBMs [29]. With the growing need for molecular tests to establish the final diagnosis of CNS tumors [11,21], the use of an IHC panel makes the subclassification of GBMs more speedy and affordable for oncology services worldwide.    In summary, glioma grading and classification is a complex, evolving field. Intertumor heterogeneity studies are the basis of the WHO classification and its regular updates. These studies need a large number of samples, but allow for the observation of common histopathological and/or molecular characteristics, as well as the discovery and validation of prognostic/predictive biomarkers. Nevertheless, noticeable variability in outcomes is still seen among individuals of a single group of tumors. Therefore, closer analysis of individual tumors is necessary to assess how intra-tumor heterogeneity influences the outcomes and responses to treatment of GBM patients. In the next sections, we will review how different methods have helped advance the understanding of tumor biology and how intra-tumor heterogeneity studies are paving the way to a more individualized care for patients with gliomas.

Spatial Heterogeneity and Spatial-Temporal Histomorphology-Based Methods
Although histopathology was the foundation for gliomas grading and classification, hematoxylin and eosin (H&E) analyses are intrinsically subjective and highly dependent on the representation of the tumor. As an example, several features used for astrocytoma grading are prone to subjectivity, resulting in up to 20% inter-observer disagreement [30]. One of the reasons for this discrepancy is that 62% of all gliomas retain areas histologically representative of grades 2, 3, and 4 within a single tumor, which enhances the chance of grade underestimation in small tissue fragments (e.g., stereotactic biopsies) (Figure 4a) [30,31]. Despite that, histopathology has guided GBM research, and several methods of molecular evaluation, such as fluorescence in situ hybridization (FISH) and IHC, require correlation with histopathological features. These methods have greatly benefitted from the development of tissue microarrays (TMAs) [32]. TMAs enable the simultaneous study of dozens of formalin-fixed, paraffin embedded (FFPE) specimens in a single block, supporting back-toback comparisons of multiple tumors, multiple areas from a single tumor, and of matched pairs of primary and recurrence tumors.
Despite the limitations, several histopathological features were associated with response to different treatments. Tumors with high cellularity were negatively associated with survival in patients with no adjuvant treatment. Not surprisingly, high cellularity presented a positive correlation in patients treated with radiation plus 1,3 bis (2-chloroethyl)-1nitrosourea (BCNU), reflecting the susceptibility of tumor cells to chemo-radiation. Necrosis, microvascular proliferation, and infiltration of inflammatory cells were not associated with survival or response to treatment [33]. Correlative studies with autopsy material and computerized tomography (CT) images, indicated that small anaplastic and small fibrillary cell populations (not the pleomorphic cells- Figure 4b) enhanced GBM invasiveness and mass effect, but also improved response to treatment. The authors also described different cell densities in the necrotic center, the contrast-enhancing rim, and in the perilesional area [34,35], which are radiological regions of interest in all image exams in the context of neuro-oncology. Although the radiological aspects of glioma heterogeneity are beyond the scope of this review, it is important to note that advances in neuroimaging techniques such as magnetic resonance image (MRI) sequences improved the guidance for stereotactic biopsies, and predictions based on texture analysis and machine learning algorithms are advancing radiomics to a new field of radiogenomics [36][37][38].
Unfortunately, some high-throughput molecular techniques demand tissue dissociation, which results in the loss of important information on the cellular and subcellular expression of altered genes. Morphology-based studies with IHC may detect the consequences of multiple molecular alterations in tumor cells. For example, gene mutations with gain of function and amplifications result in increased protein expression, while reduced expression may follow loss-of-function mutations and DNA methylation. Several authors have extensively demonstrated the expression of GBM markers, such as EGFR, MGMT, IDH1 R132H, and ATRX, in tumor cells using IHC, which is much more applicable for neuropathology routine [39][40][41][42][43][44][45][46]. With systematic evaluation of IHC stains, several authors showed that while endothelial cells do not express EGFR, regional differences in EGFR staining in GBM tumor cells may be explained by different levels of EGFR amplification [45,46], with higher EGFR amplification being associated with high proliferative activity, higher mutational burden, and shorter OS [45]. Heterogeneous MGMT expression in gliomas including positive stain in endothelial cells and lymphocytes has been described for decades [40,42,44]. Although MGMT IHC is not the gold standard prognostic test for GBMs, it is suggested that MGMT expression assessed by IHC can refine the predictive power of DNA methylation tests used in the clinical setting [41,43]. Our group assessed MGMT and GFAP co-expression with fluorescence IHC (Automated QUantitative Analysis-AQUA) and demonstrated that IHC may be a better prognostic test than the DNA methylation tests for patient stratification in GBM clinical trials [39]. On the other hand, because TP53, ATRX and IDH mutations are early events in gliomas [23], protein expression assessed by IHC is either diffusely lost (ATRX) or diffusely present (protein IDH1 R132H, p53) in tumor cells ( Figure 5). ATRX has shown heterogeneous expression in only about 20% of GBMs, associated with higher frequency of EGFR amplifications [47]. Whereas p53 IHC staining is not sensitive nor specific for TP53 mutations [21], the diffuse staining reflects its early occurrence in gliomagenesis.

Chromosomal and DNA Ploidy Analyses
Early studies associated differences in karyotypes of tumor cell clones with specific morphology and growing rate. Fast-growing fibroblast-like cells showed hyperploid karyotypes, while slow-growing astrocyte-like cells presented near-diploid karyotypes [48]. More recently, the combination of DNA ploidy and stem cell markers at single-cell level provided new insights on the evaluation of tumor clones. Considering genomic instability is a hallmark of GBM [49], DNA ploidy identified a subset of GBMs that displayed both pseudodiploid and aneuploid tumor cell clones, characterizing the clonal evolution process towards aneuploidy. Although both pseudodiploid and aneuploidy clones were able to form spheroids in vitro, the aneuploid clones instigated more aggressive tumors in vivo, with shorter OS [50]. To investigate stemness of these tumor cell clones, the authors evaluated the expression of several cancer stem-cells (CSCs) cell-membrane markers. Both diploid and aneuploid clones expressed CD56, CD90, and CD29 in similar levels; A2B5 expression was stronger in the diploid compared with the aneuploid fraction. However, CD133 and CD15 were heterogeneous, enriched either in the aneuploid or in diploid fraction tumor clones in different tumors. These findings suggested that in polygenomic GBMs the tumorigenic potential depends on DNA ploidy of CSCs [50].
At chromosomal level, early cytogenetic analyses of multiple samples obtained from one tumor diagnosed as "low grade oligoastrocytoma" demonstrated a ubiquitous tumor cell clone with gains in chromosome 7 (47XY, +7), while cells with concurrent chromosome 7 gains and chromosome 9 losses (47,XY,+7/47,idem,+del(9)(q21.2)/47,idem -9,+ 16) were deemed a subclonal population [51]. Intriguingly, the region with both chromosomal alterations exhibited higher pleomorphism and mitotic activity [51]. C-IMPACT-NOW-6 recently incorporated these early chromosomal alterations to the updated glioma classification independent of histopathology [11]. Therefore, although the chromosomal changes were focal, that tumor could be now classified as grade 4 astrocytoma or GBM (depending on IDH status), even in the absence of microvascular proliferation or necrosis [11,24,51]. Several further studies with comparative genomic hybridization (CGH) and laser cytometry showed that DNA aneuploidy, chromosome 7 gains, and losses in chromosomes 9, 10, 13, and 22 were diffuse in grade 3 and 4 gliomas, suggesting that genetic instability and cytogenetic heterogeneity, constitute early alterations in glioma progression [26,52,53]. Yet those findings were not associated with any histopathologic features [52]. Remarkably, Cancers 2021, 13, 761 9 of 22 cytogenetic evaluation of TME components, specifically endothelial cells showed that intra-tumor endothelial proliferation displayed unique chromosomal changes, with modest overlapping with the alterations observed in the tumor cells. This is believed to be due to tumor endothelial cells being generated by transdifferentiation of glioblastoma cells [54]. Only one of the tumors in the cohort exhibited the characteristic chromosome 7 gains observed in GBMs [55]. In single-cell studies, chromosome 7 gains and chromosome 10 losses were common to all 430 analyzed cells from a small GBM series. Other chromosomal alterations were observed in the cohort; however, this method may underestimate the presence of focal alterations occurring in very small percentage of tumor cells [56].
Given that alterations in receptor tyrosine kinase (RTKs) are a hallmark of 50% of GBMs [19,20], multiple groups have investigated the frequency and distribution of RTK amplifications with array comparative genomic hybridization (aCGH) with single-nucleotide polymorphism (SNP) analysis of specific gene loci in chromosomal regions [56][57][58]. EGFR amplifications were more frequent in GBMs than in LGGs, while PDGFRA amplifications were more homogenous among different tumor grades [57]. In addition to previously reported intra-tumor heterogeneity of EGFR amplifications [45,46], a subset of GBMs displayed mosaic amplifications in multiple RTKs, mainly involving EGFR, PDGFRA, and MET, with EGFR and PDGFRA co-amplification being the most common combinations, followed by EGFR and MET co-amplification. Rather than co-occurring in single cells, those RTKs were amplified in different cell populations, which overexpressed exclusively the correspondent protein [57,58]. Importantly, EGFR amplifications were associated with aggressive phenotype, represented by poor vascularization and small cell morphology, the more invasive cell type as previously described [34,35]. Conversely, PDGFRA-amplified cells were significantly associated with fibrillary morphology and more vascularized areas [57]. Post-mortem analysis of one GBM with a mosaic pattern of RTK amplifications showed that tumor cells with PDGFRA amplifications were limited to the bulk of the tumor, while EGFR amplifications were observed in the cell population present at areas of infiltration [58]. Single-cell RNA sequencing has also confirmed the possibility of mosaic expression of multiple RTK amplifications in EGFR, PDGFRA, PDGFA, FGFR1, FGF1, NOTCH2, and JAG1 in GBMs [56].
More recently, gene amplifications were linked to extrachromosomal DNA (ecDNA) elements, which likely give origin to different tumor cells sub-clones. The presence of ecDNA was more frequent in IDH-wild-type, compared to IDH-mutant GBMs, and targeted several oncogenes, namely MET, EGFR, CDK4, PDGFRA, and MYC, were maintained in recurrent tumors after adjuvant treatment [59]. The therapeutic implication of ecDNA elements as a mechanism of gene amplification was demonstrated in vivo. Tumors with MET amplification associated with high-frequency ecDNA elements had worse response to treatment with a MET inhibitor than the tumors with MET amplification without ecDNA elements, even though both presented MET overexpression, thus suggesting that ecDNA is a significant mechanism of treatment resistance to be considered in further studies [59].

Genomic and Transcriptomic Analyses
As the molecular classification of GBMs relies on both genomic and transcriptomic data, virtually all current GBM studies perform DNA sequencing and/or gene expression data for cohort characterization. Therefore, whole genome sequencing (WGS), whole exome sequencing (WES), RNA sequencing, gene expression arrays, and other methods are ubiquitous in glioma research [23,[59][60][61][62][63][64][65]. Importantly, the development of cell sorting through bar coding techniques enabled the evaluation of gene expression at single cell level [56,62,66,67].
In an evolutionary tree, truncal mutations are observed in every sequenced cancer cell, while branches represent subclonal mutations, present in only a subset of cancer cells. Although those are concepts closely associated with the definition of clonal and subclonal mutations, respectively, the distinction between truly clonal from pseudoclonal mutations is dependent on the number of tumor regions sequenced [7]. Suzuki et al. described regional and temporal heterogeneity based in the genomic profile of multiple LGGs areas, showing that approximately 60% of mutations were branch (subclonal) mutations, while only 10% were defined as truncal (clonal) mutations [23]. Based on these findings, this group was one of the first to describe the putative order of molecular alterations in gliomas incorporated in the WHO 2016 classification with IDH point mutations, TERT promoter mutations, and 1p19q co-deletion as early events. These events are followed by frequent mutations in ATRX/TP53 (in 1p19q intact tumors) or FUBP1/CIC (in 1p19q co-deleted tumors), which affect different tumor cell sub-clones, reflecting on the concept that clonal evolution of tumors occurs in both space and time [23].
Even before the genomic landscape of LGGs was described, GBM spatial-longitudinal evolution was inferred based on the frequency of the alterations between the samples by assessing DNA somatic copy number alterations (CNAs), single-nucleotide variants (SNVs), and gene expression profiles in multiple tumor areas [68]. This was one of the first articles to observe individual clones with different transcriptional subtypes coexisting within the tumors. In addition, early CNAs on chromosomes 7, 9, and 10, loci of several GBM drivers (EGFR, CDK6, CDKN2A/B, MET, and PTEN) were diffusely seen in the tumors. Subsequent alterations affected genomic regions containing PDGFRA (chromosome 4), PTEN (chromosome 10), and TP53 (chromosome 17), and were observed in limited/unique areas. With that, the authors achieved temporal and spatial reconstruction of GBM ontogeny [68]. In fact, the genetic distance between tumor areas was observed by other authors, and is thought to be correlated with the physical (Euclidian) distance between the tumor areas [64].
Vinci et al. further confirmed the clonal nature of H3K27M, ATRX, and NF1 mutations, as well as the subclonal expansion of cells bearing mutations in TP53, BRAF, PDGFRA, among others, in pediatric GBM (pGBM) and DIPG [61]. The authors identified differences in cell morphology, growth, migration, and invasion directly associated with genomic profile of the tumor sub-clones, even when arising in rare sub-clone variants. For example, rare mutations in the histone H4 methyltransferase KMT5B, enhanced invasiveness, and sensitivity of tumor cells to PARP inhibitors in vitro and in vivo. Co-culture of KMT5B mutant and KMT5B wild-type cells and monoculture using conditioned medium from the parental cells (bulk cells) showed significant enhancement of migration and invasion compared with single clone culture, both in vitro and in vivo, and suggested the importance of cell-cell interaction and paracrine signaling as mechanisms in tumor development. This suggests that in DIPG, multiple sub-clones cooperate to enhance tumorigenic phenotypes [61]. The concept of 'cooperative invasion' had previously been described in other cancers and in GBM [70].
In adult GBM, multiple studies have demonstrated heterogeneous transcriptome profiles within single GBMs using single-cell RNA sequencing [56,62,71]. Patel et al. defined four meta-signatures in tumor cells enriched for genes associated with cell cycle, hypoxia, and complement/immune response that correlated with the molecular subtype of GBM in single cells, which did not always correspond to that assigned in the bulk tumor [56,62]. The gene meta-signature associated with cell cycle was expressed in a small proportion of cells in each tumor and was negatively correlated with stemness and with hypoxia gene signatures. This stemness meta-signature was observed in cells of proneural and classical subtypes, but not in the mesenchymal subtype. The variability of the cell cycle meta-signature among tumor cells may be associated with variations in the TME [56]. Finally, multiple groups have reported the occurrence of hybrid tumor cells expressing both classical and mesenchymal or proneural and classical signatures that are associated with worse OS [56,71]. Intra-tumor heterogeneity defined in single-cell analyses may present prognostic implications. Upon the evaluation of bulk proneural GBMs, the authors observed the proportion of tumor cells of alternate subtypes influences the clinical outcome, with increased heterogeneity being associated with worse OS [56]. Meyer at al. demonstrated the presence of TMZ-resistant clones in a GBM not submitted to previous treatment, suggesting that resistant clones do not arise solely from therapeutic selection. Similarly, a TMZ-sensitive subclone was isolated from an otherwise recurrent GBM, which presented several TMZ-resistant clones with significant difference in the sensitivity to multiple target therapies [71]. Recently, Akgul et al. also reported that multiple molecular differences in subclones reflected different responses both to the standard GBM treatment (TMZ + radiation) and to multiple targeted therapies [62]. These findings justify the concurrent use of different combination therapies, especially in recurrent GBMs [62,71].
Because advanced technologies may detect subtle differences between tumor cells, using multiple samples from single tumors can allow for the study of spatio-longitudinal heterogeneity [23,26,50,53,67,72,73]. In this setting, studies on multicentric/multifocal GBM (M-GBM) and M-LGG have also shed light on the clonal evolution of these tumors [74][75][76]. M-GBMs and single GBMs share common molecular alterations. Similar to other reports, in M-GBMs, chromosomal alterations are early events, as shown by shared CNAs present in all areas of the tumors. On the other hand, specific gene mutations unique to different foci point toward a parallel genetic evolution from a common tumor precursor clone (early separation of cell clones) [74]. Hayes et al. described different IDH1 point mutations (R132H and R132C) in different tumor areas of one M-LGG. The patient had a germline mutation on TP53 present in only in the tumor foci with the non-canonical IDH1 mutation. In another patient, two foci of M-LGG presented both grade 2 and 3 histopathologic grades. In this patient, multiple TP53 and ATRX point mutations were conserved among samples, but distinct between the two grades. Finally, two patients in this series presented M-LGG with divergent ATRX/1p19q status between the tumor foci, challenging the discouragement of the diagnosis of "oligoastrocytoma" in tumors with mixed histopathological aspects [75]. These findings confirm the divergent evolution in multifocal tumors and the need for careful evaluation of small samples that may misrepresent the whole tumor.

DNA Methylation Analyses
Epigenetic characterization of GBMs has provided great insight into inter-tumoral heterogeneity and glioma classification and prognosis [15,16]. DNA methylation of CpG islands in the gene promoter are associated with repression of gene expression [77]. These bear a close relationship with hypermethylated status conferred by IDH mutations in gliomas [15,16]. The methylome of tumors has been integrated in most heterogeneity papers using high-throughput methods, especially for temporal heterogeneity studies. DNA methylation has also advanced the knowledge in intra-tumor spatial heterogeneity, especially in the characterization of CSCs and inflammatory cells [78][79][80][81]. Notably, differences in DNA methylation between tumor areas as close as 5 mm apart were observed, both in IDH mutant and IDH wild-type gliomas [64]. For specific information about DNA methylation and its relationship with other epigenetic modifications in cancer, we recommend a recent review covering this topic [77].
Comparisons between the DNA methylome of secondary astrocytomas IDH mutant grade 4 (former IDH mutant GBMs), and their matched primary LGG found that the G-CIMP profile present in all primary tumors was maintained in the recurrences. However, the malignant progression of gliomas toward GBM showed loss of DNA methylation specifically in cell cycle-related genes [73]. The changes in the DNA methylation profile reflected the phylogenetic tree of gene mutations among tumor areas from a single tumor and between matched primary and recurrent tumors. However, in one patient with three LGG recurrences, the genomic changes occurred prior to (but in co-dependency with) epigenomic alteration in the glioma recurrences with grade progression. Therefore, the phyloepigenetic tree was slightly divergent from the phylogenetic tree of the same samples [73].
Klughammer et al. used reduced representation bisulfite sequencing (RRBS), a novel method that "allows for individual cell assessments, without the need for single-cell sequencing, and precludes the need for RNA extraction", to assess temporal and spatial heterogeneity in newly diagnosed and recurrent GBMs [72]. The authors reproduced the transcriptome subclassification of GBMs with high confidence, including the identification of multiple subtypes in tumor cell clones [56,71], and showed significantly hypermethylation of chromatin binding proteins (CTCF, EZH2, and KDM4A) and reduced methylation of stemness regulators (OCT4, NANOG, SOX2) in the mesenchymal subtype. Their pathway analysis identified enrichment of neural development and apoptosis signaling among genes whose promoters gained DNA methylation during disease progression, while genes whose promoters lost DNA methylation were enriched in the Wnt signaling pathway and T-cell activation [72]. In addition, using machine learning methods, specific DNA profiles were associated with GBM subtypes, the presence and density of inflammatory cells (macrophages and lymphocytes), histopathological features (necrotic area, cell morphology, and mitotic activity), and prognosis. For example, mesenchymal GBM was more associated with infiltration of inflammatory cells, lower cell density and extensive necrotic areas. These characteristics were also more significant in recurrent tumors, compared to primary tumors [72].
Pangeni et al. characterized the DNA methylation profile of GSCs, and observed higher frequencies of methylated CpGs than GBM bulk tumor cells, especially in mesenchymal GBMs. Specific DNA methylation signatures were observed in CSCs and in GBM bulk tumor, with only few overlappings, in both proneural and mesenchymal GBMs [82]. Peculiarly, these authors described predominance of global hypomethylation events in proneural GBMs, which are usually IDH-mutant, therefore hypermethylated tumors. However, the proneural GBMs of their series were reportedly IDH wild-type tumors, less likely to present high G-CIMP [16,82]. More specifically, because GSCs are pluripotent cells with self-renewal and proliferation properties, the expression of EZH2, a protein related to maintenance of stemness is of high importance, especially in H3K27-mutant tumors. Overexpression of EZH2 in HGGs is associated with worse OS. This gene is part of the polycomb repressive complexes involved in the epigenetic modulation of GSCs, the most studied epigenetic modulators in GSCs, and is therefore a potential therapeutic target [80].
Recently, Feng et al. proposed a signature based on the methylome profile composed by five candidate probes closely associated with the DNA damage repair function as an independent predictive marker for radiation therapy in gliomas [79]. Similarly, Wang et al. proposed a methylation-based classification of GBMs to predict response to immunotherapy [81]. They correlated DNA methylation and gene expression to define two clusters of GBMs with different prognosis, expression of immune checkpoint molecules, and response to immunotherapy, independent of other clinic prognostic factors. Six novel candidate genes (ANKRD10, BMP2, LOXL1, RPL39L, TMEM52, and VILL) were selected-based in survival analyses, to propose a reliable prognostic model, which resulted in an effective nomogram for survival prediction in GBM that was likely the first to incorporate global methylation signature [81].   The interactions between the TME components with the tumor cells are increasingly being explored and may determine response to immunotherapy and enable the development of new targeted therapies. Histologically, it is not difficult for the pathologist to identify and characterize, with H&E and IHC stains, the presence of microvascular proliferation and inflammatory infiltrate in GBMs. The characterization of macrophages relies in the expression of CD68, CD163 (marker of M2 polarized macrophages- Figure 6b), and IBA-1, among others. The difference between macrophage and resident microglia is based on high and low expression of CD45 (Figure 6a), respectively, and CD11b expression by microglia. The lymphocyte population is characterized by the expression of CD45, CD3, CD4, CD8, FOXP3, and more recently PD-1 [83,84]. The analysis of inflammatory cells has gained momentum with the development of cell sorting methods and single-cell analyses, improving the potential for use of immunotherapy and identification of new therapeutic targets. In addition, the analysis of cytokines and their effects in tumor initiation and progression has advanced with other methods not explored here, but is reviewed by other authors [85].
Several authors have identified that GBMs are 30-50% composed of microglia/macrophages (Figure 6b) [84,[86][87][88]. However, tumor-associated macrophages (TAMs) do not retain the ability to activate immune response through T-cells, because while expressing MHC class II, TAMs lack the expression of co-stimulatory factors (CD80, CD86, CD40) and most cytokines (mainly IL-1, IL-10, and TNF-α) needed for adaptive immunity [84]. Over the years, the more comprehensive understanding of macrophage functions led to subclassification according to their pro-inflammatory and anti-tumor properties as M1 and M2, respectively, observed in vitro [89,90]. However, these definitions do not translate well to the in vivo setting, and TAMs express marker characteristics of both M1 and M2 phenotypes, with the ability to transition between phenotypes, depending on the interactions with the TME [83]. Increasing evidence suggests that TAMs promote tumor invasion and growth (through several factors released from microglia, such as STI-1, EGF, CSF-1, CCL2, and TGF-β), and angiogenesis, through VEGF expression. The expression of TGF-β by TAMs is also associated with increased GSC invasiveness, whereas GSC recruit TAMs more efficiently than differentiated tumor cells. An in-depth review of the complexity of the interactions between the different TAMs phenotypes and between TAM and other components of gliomas may be seen in other reviews [83,89].
Tumor-infiltrating lymphocytes (TIL) usually correspond to up to 10% of glioma content, represented by T-cells, natural killer (NK) cells and B lymphocytes [83,87,88]. Nevertheless, their presence is so meaningful that it is the basis for another classification of solid tumors as "hot" (inflamed) and "cold" (non-inflamed) tumors based on an immunoscore initially developed for colorectal cancer that accounts for the presence of CD3 and CD8, markers of T-lymphocytes, in the tumor core and margins [91,92]. In general, GBM is considered a "cold" tumor, with important diagnostic/prognostic and therapeutic consequences [93]. Converting GBM from a "cold" to a "hot" phenotype is one of the goals in improving response to immunotherapies. For example, the metabolic product of mutant IDH1/2, 2-hydroxiglutarate (2-HG), suppresses antitumor activity of tumorinfiltrated cytotoxic T-cells in vivo [94], and therefore, IDH mutations are, per se, targets for new vaccines. In addition, activation of the PI3K pathway in GBM cells results in the expression of the immune checkpoint ligand PD-L1 ( Figure 6c); additionally, activation of the Ras-MAPK pathway induces the expression of IL6 and TGF-β, which enhances the recruitment of TAMs and inhibits both dendritic cells (DC) and lymphocytes infiltration into the tumor. This event is now seen as a potential avenue for development of multiple targeted interventions [93]. Single-agent neo-adjuvant immunotherapy with check point inhibitors (PD1 axis inhibitor, nivolumab) has failed to show clinical benefits in GBM [95]. However, neoadjuvant use of PD-L1 inhibitor resulted in significant modification in the TME towards a stronger anti-tumor immune response [95,96].
TILs may affect the outcomes of GBMs subtypes [67,72,78,82,97,98]. Wang et al. described the specific associations of transcriptomic subtypes of GBM with immune cells gene signatures [67]. Proneural GBMs were associated with a reduced CD4+ T-cells gene signatures, the classical GBMs were associated with an enhanced dendritic cell (DC) signature, and the mesenchymal GBMs were associated with enhanced M1 macrophage and neutrophil, but decreased NK cells gene signatures [15,67]. Other authors have reported increased infiltration of immune cells in mesenchymal type GBM compared to proneural and classical subtypes, using different assessment methods [72,78]. More specifically, mesenchymal tumors showed the highest amount of tumor-infiltrating CD3+ and CD8+ T-cells, while IDH mutant proneural tumors showed the lowest amount of CD3+ and CD8+ T-cell infiltration. Expression of CD68, FOXP3, and PD-1 did not show differences between the subtypes [78]. Importantly, these authors compared matched samples of patients before and after treatment with TMZ + radiation and an experimental DC vaccination. While CD3+ T-cells did not show differences over time, there was a reduction in CD4+ T-cells and an increase in CD8+ PD-1+ T-cells in all methylation types. Macrophage infiltration did not show differences between the subtypes with a general IHC marker (CD68) unable to show differences in M1/M2 macrophages. Although direct effects in OS were not detected, the findings confirm that the more aggressive mesenchymal GBMs presented higher T-cell infiltration, conferring a more immunogenic nature to those tumors [78]. To demonstrate the interaction between glioma tumor cells and immune cells, Zhai et al. showed that mRNA expression of Indoleamine 2,3-dioxygense 1 (IDO1), an immunossupressive enzyme, increases in tumor cells with histologic grade in gliomas. IDO1 showed increased expression in mesenchymal and classical GBMs compared with proneural GBMs, and was higher in IDH wild-type than in IDH mutant gliomas. Thus, high IDO1 expression was related to worse OS. Importantly, IDO1 expression in tumor cells increased with the increased expression of markers of TAMs (such as HLA-DRA and CD163) and neutrophils (for example, CD11b and CD16), and with multiple lymphocitic markers, including FOXP3, PD-1 and PD-L1, suggesting that T-cells directly regulate IDO1 expression in human GBM. Unfortunately, the IHC detection of IDO1 was minimal, and the evaluation of this important marker was only possible with mRNA measurements [98].
In pediatric HGGs, Bockmayr et al. identified four groups of tumors on the basis of the gene signatures of infiltrating immune cells-vascular (immune 1); monocytic and stromal (immune 2); monocytic and T-cell dominant (immune 3); and antigen-presenting cells (APC), NK cells, and T-cells (immune 4)-that were associated with the transcriptomic molecular subtypes. The vascular signature was associated with the classical subtype, while the mesenchymal subtype was associated with monocytes and T-cell gene signatures (immune 2 and 3). Interestingly, the immune 4 signature (innate immunity) was enriched in Pediatric HGGs with G34 mutations, a histone H3-mutated tumor with better prognosis than other pediatric HGGs [97]. This reflects the anti-tumor effects of the innate immune response.
Tumor immunology is an extremely complex field, with the same immune cells being able to perform anti-and/or pro-tumoral functions. Therefore, the development of combinatory therapies with novel drugs targeting multiple candidate molecules along the chain of events in the inflammatory response to tumors, rather than use of single-drug approaches, is an urgent need in GBMs.
A summary of the histopathological, immunohistochemical, and molecular features of grade 4 IDH mutant astrocytomas and different subtypes of GBMs is shown in Table 1.

Longitudinal Heterogeneity
Although maximum surgical resection has provided significant survival benefits for patients with recurrent GBM [99,100], some controversy over operating on patients with recurrent tumor still exists. With this in mind, studies on GBM longitudinal heterogeneity are needed to better understand the differences between primary and recurrent GBMs. Two recent comprehensive articles have begun to look at this by comparing patient's samples collected at multiple time points and have produced valuable findings regarding GBM evolution [53,101].
The Glioma Longitudinal Analysis (GLASS) consortium characterized the largest cohort of matched primary and recurrent gliomas, including all three WHO classes [101]. Overall, the mutational burden was low in primary gliomas of all classes, but increased in 70% of recurrent tumors. Hypermutational state was associated with treatment with alkylating agent, and significantly more frequent in IDH mutant-noncodel gliomas (47%), than in IDH mutant-codels (25%) and IDH wild-type gliomas (16%). The clonal structure at the primary tumors mostly persisted into recurrences, and the neutral-to-neutral trajectory (random mutations that are not subject to selection) was the most common evolutionary trajectory overall. Nevertheless, selection model was more frequent in IDH wild-type gliomas at any time point. Selection evolution at recurrence was associated with shorter OS than neutral evolution in all groups [101]. The order of IDH1, ATRX, and TP53 mutations followed previous reports [23]. On the other hand, in IDH wild-type tumors, their findings indicated chromosome 10 deletions as earlier events than chromosome 7 amplifications. A gross comparison of all shared mutations revealed that cancer cell fractions (CCF) increased in recurrent IDH-mutant-noncodel gliomas and decreased in recurrent IDH wild-type glioma, which suggests, respectively, reductions and increases in intra-tumor heterogeneity in recurrent tumors [101]. Concerning specific drivers of glioma progression, recurrent IDH mutant-noncodel gliomas presented an enrichment in CDKN2A homozygous deletions, with subsequent DNA aneuploidy and genetic instability that were associated with shorter OS [101].
Finally, because 80% of LGGs harbor IDH1 or IDH2 mono-allelic mutations [98], secondary GBMs (GBMs arising from previous LGGs) were often thought to be virtually always IDH mutant lower grade gliomas [21]. Korber et al. evaluated a cohort of matched primary and recurrent IDH wild-type gliomas [53], and overall, contrary to the previous report [101], mutational burden was similar in primary and recurrent tumors, and 2/3 of recurrent tumors originated from multiple clones of the primary tumor (oligoclonal origin), while 1/3 were monoclonal. Additionally, although the heterogeneity in recurrent tumors was inherited from the primary lesion, most sub-clones in recurrent tumors showed additional driver mutations not found in the primary lesion, indicating ongoing genetic evolution [53]. CNVs in chromosomes 7, 9, and 10 were early alterations and preceded SNVs, similar to several previous studies in spatial heterogeneity [23,26,51,52]. Among the SNVs, TERT promoter, TP53, PTEN and EGFR mutations were the most frequent clonal alterations. TERT mutations exhibited a survival advantage over other early mutations, such as CDKN2A/B and PTEN. More importantly, through high-level mathematical estimates, the authors were able to propose that an IDH wild-type GBM may grow for 2-7 years undetected [53]. Thus, early diagnosis and surgical approaches may potentially identify a primary IDH wild-type LGG that evolves in a secondary GBM.

Future Directions
Over the next decade, the investigation of intra-tumor heterogeneity in GBM will likely increase in scope and detail. Proteomics and metabolomics are areas with opportunities for expanded understanding of gliomas. Nevertheless, inter-and intra-tumor differences in proteome profile have been reported in small series of adult patients [102][103][104], and in pediatric neuro-oncology [105]. Multi-institutional studies in adult populations are needed, aiming to identify molecular targets not observed from RNAseq data, as gene expression data often does not correlate highly to changes in proteomic expression. When paired with genomic sequencing, proteomic evaluation can better tease apart protein-level effects of gene mutations in GBMs, and help identify new IHC markers to be implemented for routine neuropathological testing. In addition, multiregional sampling of GBMs may discriminate functional differences between non-necrotic, peri-necrotic, and necrotic tissue, and potentially yield insights into mechanisms underlying treatment resistance related to changes in the TME. Batch effects in mass spectrometry (MS), and the lower quality of FFPE samples compared to fresh-frozen tissue is still an issue in proteomic studies, particularly in phospho-proteomic studies. We have demonstrated quality improvement in MS analysis of FFPE samples, which enabled the identification of differential proteomic profiles between IDH wild-type and IDH mutant gliomas, as well as among multiple areas of LGGs and GBMs, reflecting different molecular pathways [106,107]. Similar to other methods, it is expected that proteomics will advance for analyses at single-cell level.
The rapidly expanding field of metabolomics is another area of study where opportunities exist to better understand glioma biology. Recent advances, such as desorption electrospray ionization mass spectrometry (DESI-MS), can allow for finely detailed spatial resolution of lipids and metabolites from biological samples [108,109]. In coming years, the 2D slices of DESI-MS will likely be expanded in order to recreate the 3D metabolic structure of tumors through combination of individual 2D slices. This may better show how gliomas secure nutrients and other metabolites crucial for continued growth, which could be used for the development of new markers for functional image exams. Given the phenotypic differences given by IDH mutations in gliomas and its product 2-HG, which interferes in the Kreb's Cycle, there is opportunity for exploration of the metabolic consequences of these mutations for GBMs.
The use of new methods to characterize tumor heterogeneity in GBM at protein and metabolite levels will fuel the next revolution in functional image examinations, which will trigger further improvements in early diagnosis, efficacy of surgical resection, histomolecular classification, and drug development, thereby substantially contributing to the advancement of personalized care for patients with GBM. Institutional Review Board Statement: Not applicable. All the microphotographs depicted in this review were captured from The OSU Tissue Archives (IRB 2013C0020).

Informed Consent Statement: Not applicable.
Data Availability Statement: Not applicable.