Multiregional Sequencing of IDH-WT Glioblastoma Reveals High Genetic Heterogeneity and a Dynamic Evolutionary History

Simple Summary Glioblastoma is the most common and aggressive primary brain malignancy in adults. In addition to extensive inter-patient heterogeneity, glioblastoma shows intra-tumor extensive cellular and molecular heterogeneity, both spatially and temporally. This heterogeneity is one of the main reasons for the poor prognosis and overall survival. Moreover, it raises the important question of whether the molecular characterization of a single biopsy sample, as performed in standard diagnostics, actually represents the entire lesion. In this study, we sequenced the whole exome of nine spatially different cancer regions of three primary glioblastomas. We characterized their mutational profiles and copy number alterations, with implications for our understanding of tumor biology in relation to clonal architecture and evolutionary dynamics, as well as therapeutically relevant alterations. Abstract Glioblastoma is one of the most common and lethal primary neoplasms of the brain. Patient survival has not improved significantly over the past three decades and the patient median survival is just over one year. Tumor heterogeneity is thought to be a major determinant of therapeutic failure and a major reason for poor overall survival. This work aims to comprehensively define intra- and inter-tumor heterogeneity by mapping the genomic and mutational landscape of multiple areas of three primary IDH wild-type (IDH-WT) glioblastomas. Using whole exome sequencing, we explored how copy number variation, chromosomal and single loci amplifications/deletions, and mutational burden are spatially distributed across nine different tumor regions. The results show that all tumors exhibit a different signature despite the same diagnosis. Above all, a high inter-tumor heterogeneity emerges. The evolutionary dynamics of all identified mutations within each region underline the questionable value of a single biopsy and thus the therapeutic approach for the patient. Multiregional collection and subsequent sequencing are essential to try to address the clinical challenge of precision medicine. Especially in glioblastoma, this approach could provide powerful support to pathologists and oncologists in evaluating the diagnosis and defining the best treatment option.


Introduction
Glioblastoma (GB) is the most common and aggressive primary brain malignancy in adults and one of the deadliest human cancers [1][2][3]. Mean survival rates are around 14 months with current and aggressive therapeutic modalities that include maximal safe, surgical resection, followed by radiotherapy and concomitant and adjuvant chemotherapy with the alkylating agent temozolomide [2][3][4][5]. Although GB is a rare cancer with an overall incidence of less than 10 per 100,000 people, its poor prognosis makes it a crucial public health problem [6]. Malignant gliomas are the cause of 2.5% of cancer deaths and are the third leading cause of cancer death in people between 15 and 34 years of age [7]. Glioblastomas are divided in the 2016 World Health Organization's WHO Classification of Tumors of the Central Nervous System [5,8] into two groups, based on genetic mutations in isocitrate dehydrogenase genes (IDH1 and IDH2): IDH-WT GB (90% of cases), defined as primary or de novo glioblastoma, and IDH-mutant GB, called secondary glioblastoma, with a history of previous lower-grade diffuse glioma. Overall survival (OS) in IDH1-mutant GB is more than three times higher than in IDH-WT GB [9].
The impossibility of a total tumor debulking and a poor distribution of drugs in the brain, due to the presence of the blood-brain barrier (BBB), contribute significantly to the lack of effective treatment options and a poor prognosis [1,3]. However, clinical results vary considerably between patients. Previous studies have shown marked differences between tumors at the genomic and transcriptomic level, which may be the basis for differences both in a patient's natural tumor history and in treatment responses [2]. In addition to inter-patient differences, GB evolution results in spatial and temporal intratumor heterogeneity [1,10].
Tumor heterogeneity is one of the main reasons for the poor prognosis and overall survival due to therapeutic failure and drug-resistance [1,3,[11][12][13]; therefore, it represents a challenge to achieve the goals of precision medicine. This raises the important question of whether the molecular characterization of a single portion of the tumor sufficiently represents the genomic landscape of GB in a biologically and clinically significant way for personalized-medicine approaches. Therapy selection based on the analysis from a single biopsy specimen may not be representative of the entire lesion and could result in treatment failure [2,14,15].
To face these challenges and fully define the mutational landscape and intra-tumor heterogeneity of GBs at the patient level, we sequenced the whole-exome of nine spatially different cancer regions of three primary IDH-WT GBs. Our results show that tumor heterogeneity displays a specific signature that reveals the evolutionary dynamics of GB at the individual patient level, highlighting the need for multiregional sequencing for precise therapy.

GB Patients
Three primary human GB surgical specimens (GB01, GB02 and GB03), diagnosed according to WHO diagnostic criteria [8], were collected from the Department of Neurosurgery (University of Pisa). Tumors were resected by the same surgeon and reviewed by the same pathologist. Patients received similar medical treatments. All cases had a diagnosis of GB with no previous history of any brain neoplasia and were not carrying R132 IDH1 or R172 IDH2 mutations or the 1p/19q codeletion. Patient clinical and demographic data are shown in Table 1. The study was approved by the Ethics Committee of the University Hospital of Pisa and all methods were performed in accordance with approved guidelines. Patients' data and samples were completely anonymized.

Sample Collection
All surgically resected tumor samples were chosen by the neurosurgeon during surgery. Spatially separated tumor regions, documented by photography, were cut with a scalpel and collected immediately after surgery under a biological hood. On the surface of each tumor, nine spatially organized samples (T1-T9) of approximately 3 mm 3 each, were collected with locations as shown in Figure 1.

GB Patients
Three primary human GB surgical specimens (GB01, GB02 and GB03), diagnosed according to WHO diagnostic criteria [8], were collected from the Department of Neurosurgery (University of Pisa). Tumors were resected by the same surgeon and reviewed by the same pathologist. Patients received similar medical treatments. All cases had a diagnosis of GB with no previous history of any brain neoplasia and were not carrying R132 IDH1 or R172 IDH2 mutations or the 1p/19q codeletion. Patient clinical and demographic data are shown in Table 1. The study was approved by the Ethics Committee of the University Hospital of Pisa and all methods were performed in accordance with approved guidelines. Patients' data and samples were completely anonymized.

Sample Collection
All surgically resected tumor samples were chosen by the neurosurgeon during surgery. Spatially separated tumor regions, documented by photography, were cut with a scalpel and collected immediately after surgery under a biological hood. On the surface of each tumor, nine spatially organized samples (T1-T9) of approximately 3 mm 3 each, were collected with locations as shown in Figure 1. Single portions were then snap-frozen in liquid nitrogen for subsequent DNA extraction.

Sample Processing
Each tumor region was used for genomic DNA extraction using a modification of the Maxwell 16 FFPE Tissue LEV DNA Purification Kit (Promega, Madison, WI, USA). Specifically, incubation with Proteinase K (Promega, Madison, WI, USA) and Incubation Buffer (Promega, Madison, WI, USA) was performed for one hour at 56 °C. DNA concentration was determined using the Qubit Fluorometer (Life Technologies, Carlsbad, CA, Single portions were then snap-frozen in liquid nitrogen for subsequent DNA extraction.

Sample Processing
Each tumor region was used for genomic DNA extraction using a modification of the Maxwell 16 FFPE Tissue LEV DNA Purification Kit (Promega, Madison, WI, USA). Specifically, incubation with Proteinase K (Promega, Madison, WI, USA) and Incubation Buffer (Promega, Madison, WI, USA) was performed for one hour at 56 • C. DNA concentration was determined using the Qubit Fluorometer (Life Technologies, Carlsbad, CA, USA) and the quality was assessed using the Agilent 2200 Tapestation (Agilent Technologies, Santa Clara, CA, USA) system.

Multi-Region Whole-Exome Sequence (WES) Analysis
For each tumor portion, an exome library was prepared from 50 ng DNA using Nextera Exome Kit (Illumina, San Diego, CA, USA) following the manufacturer's instructions. Each NGS run included 9 pooled libraries loaded into one NextSeq High Output cartridge (300 Cycles; Illumina). Paired-end sequencing was performed on a NextSeq 500 system (Illumina) with 151 bp sequencing.

Bioinformatic Analysis and Data Interpretation
Quality control was performed on fastq files using FastQC (v011.9) and FastqScreen (v0.14.1) [16]. Paired-end reads were aligned to GRCh37/hg19 using BWA-MEM [17] aligner algorithm with default parameters. Mapped reads were sorted with SMA and MarkDuplicates tools (http://broadinstitute.github.io/picard, accessed on 16 March 2020). Target coverage of the exome sequencing of all analyzed samples averaged 70.4X. Aligned reads were processed using GATK [18] to remove low mapping quality reads (MPQ ≥ 20) and realigned in the genomic regions around potential indels. Base quality scores were recalibrated for the BAM files using GATK.

Whole-Exome Sequencing
A total of 27 intra-tumor regions, derived from three different human primary GBs, were subjected to whole exome sequencing (WES). Data analysis showed in coding regions an average of more than 33, 31 and 29 million of mapping reads, in GB01, GB02 and GB03, respectively. Total mean coverage was 70.36X (SD = 2.39). Further information on the run is shown in Table S1. Nucleic acid sequence data are available from The European Genome-phenome Archive (EGA): https://ega-archive.org, accessed on 8 March 2021) (submission ID: EGAD00001007063).

Copy Number Variation Analysis
Copy number variation (CNV) analysis was carried out as a first analysis to provide an overview of the chromosomal structural abnormalities and genomic instability, both in the whole tumor and in the single intra-tumor portions. Through CNApp [23] we calculated the chromosomal alterations of the entire tumor ( Figure 2A).
As shown in Figure 2A the alteration (gain or loss) of the short (p) and long (q) arm of each chromosome is displayed as a percentage variable from 0 to 100%. This percentage represents the number of tumor portions out of the total (nine portions) that are affected by the specific alteration ( Figure 2A). In addition, with the CNApp, we also clustered the 27 intra-tumor regions ( Figure 2B), to identify the differences and similarities in the copy number variations among regions of the same tumor and different tumors. The CNV analysis revealed a considerable aneuploidy with whole-arm and whole-chromosome alterations (gains and losses), demonstrating chromosomal instability across our samples (Figure 2A,B). The most evident chromosome alterations present in all three GB samples are chromosome 9 and 10 deletion and chromosome seven amplification. The long arm of chromosome 10 has major losses in all 27 tumor regions (Figure 2A,B). The short arm of chromosome 10 is deeply deleted in all regions of GB01 and GB02 while it has a smaller deletion in the GB03 portions (Figure 2A,B). 9p is also deleted in all 27 tumor regions, whereas 9q has deletions in only few regions (eight/nine in GB01, five/nine in GB02, and three/nine in GB03) (Figure 2A Figure 2B). Other common chromosomal alterations observed in the three different tumors are 13q deletion (all regions of GB03 and region 09 of GB01) and 19q amplification (seven/nine and six/nine regions of GB01 and GB02, respectively) ( Figure 2B). GB03 shows loss of the 19p arm, whereas in some regions of GB01 and GB02 are amplified (Figure 2A,B). GB01 has a loss of 20q while in GB02 shows a gain of the entire chromosome (Figure 2A,B). GB01 displays increased numbers of altered DNA copies, especially in chromosomes 3, 8, 15, 16, and 22 (Figure 2A,B). Notably, all regions of GB01 share common alterations in chromosomes 15 and 22, with a deletion of the long arm (Figure 2A,B). Three out of nine regions of GB01 (GB01-04/05/06) have, on chromosomes 3 and 17, common gene copy number alterations distributed throughout the chromosomes ( Figure 2B). On chromosomes 8 and 16, one/nine regions of GB01 (GB01-09) has exclusive gene copy number deletions located on both the long and short arms ( Figure 2B).
Genome-wide DNA copy number alteration profiles have previously been correlated with the four molecular subtypes of GB: Classical, Mesenchymal, Neural, and Proneural, according to the study initiated by Phillips et al. [26] and the subsequent Verhaak classification [27]. This classification is based on fourteen amplifications and seven homozygous or hemizygous deletion events, both broad and focal, on chromosomes 4 (4q12); 7 (7p11.2, 7q21.2, 7q31.2 and 7q34); 9 (9p21.3); 10 (10q23); 13 (13q14) and 17 (17q11.2). Considering the aforementioned chromosome alterations, we used the CNApp classifier [23] to calculate the global score (GCS) that the system assigns during re-segmentation by classifying and weighting CNVs according to their length and width. Next, the GCS was used to correlate the 27 intra-tumor regions with the 480 GBs from the TCGA-GB data collection (https://www.cancer.gov/tcga, accessed on 4 December 2020), to determine which GB molecular subclass they belonged to or resembled the most. The correlation analysis is shown in Figure 2C, which consists of a clustering dendrogram of our 27 samples based on their GCSs (upper section) and a table with the correlation coefficients related to the four molecular subclasses and CGS values (lower section). Based on the correlation coefficients, the classifier assigns a molecular subclass (third-to-last column, Classifier prediction group, Figure 2C). In cases where the system was unable to assign a sample to any of the four molecular subclasses, we manually associated the sample to the class with the highest correlation coefficient (second-to-last column, Most related molecular subclass for NC predicted regions, Figure 2C). The analysis revealed a highly variable CGS (−1.3 to 3.8) both among the three tumors and among regions of the same tumor, further confirming the existence of distinct CNV profiles (last column, Figure 2C). Fourteen out of twentyseven regions were correlated by the classifier to one of the four GB molecular subclasses ( Figure 2C). Eleven out of thirteen regions, not automatically correlated by the system, were manually associated with the molecular class for which they had the highest correlation coefficient ( Figure 2C). Two regions of GB02 (01 and 02) were not correlated by the classifier or manually associated with a specific molecular subtype because they had no correlation with any of the four molecular subclasses (correlation coefficient 0) ( Figure 2C). GB01 was mostly assigned to the Classical subtype (five/nine regions), while three/nine regions of GB02 are associated with the Neural subclass. GB03 was mostly assigned to the Mesenchymal subclass (four/nine) ( Figure 2D). Nevertheless, for all three samples, within the same tumor, we found different molecular subtypes ( Figure 2D).  in the classifier as "not classified" (NC) and correlation was performed using the global score (GCS) that CNApp assigns during resegmentation, by which it classifies and weights CNVs based on their length and width. The values in the table are the correlation coefficients that each region has with each of the four molecular subclasses and with the NC class. The "Classifier prediction group" column reports the molecular subtype that the system has correlated and predicted. The "Most related molecular subclass for NC predicted regions" column reports the best correlation found with one of the four molecular subclasses when the tumor region was associated with NC by the classifier. (D) Pie chart of the molecular composition of each tumor. Each tumor is composed of several molecular subclasses, each associated with an intra-tumoral region.

Mutational Landscape
Initially we analyzed the genomic alterations by selecting all mutations in the coding region (known to alter protein function), as well as all splicing mutations, excluding synonymous mutations that had no predicted impact on splicing. We discovered a total of 97,732, 50,562 and 35,498 variants for tumor GB01, GB02 and GB03, respectively.
The variant classification plot ( Figure 3A), shows the number and classification of the variants in each tumor region, providing percentages of missense, nonsense, splice-site, and frame shift mutations (both insertion and deletion) out of the total number of mutations.  (D) Oncoplot of the distribution of mutations found in our samples in the most frequently mutated genes in GB. Each column represents one sample (9 regions per tumor) and each row a different gene. Colored squares show mutated genes, while empty (gray) squares show no mutated genes. The different types of mutations are colored according to the type of variant: orange, splice site mutation; blue, frameshift deletion; green, missense mutation; red, nonsense mutation; and black, multi-hit mutation. Genes annotated as "multi-hit" have more than one type of mutation in the same region. The bar graph on the right shows the total number and percentage of mutated regions for each gene, out of the total 27 regions, colored according to mutation type. The upper graph shows the total number and type (different color) of mutations for each tumor region.
GB01 has fewer missense mutations overall than GB02 and GB03, counting more frameshift deletions and nonsense mutations. In Figure 3B, mutations are classified by type (single nucleotide variants-SNVs, deletion, and insertion), represented as a percentage of the total in the whole tumor. Again, the mutation type classification of GB01 differs from the other two tumors, having fewer SNVs and more deletions than GB02 and GB03, which instead maintain a more similar profile to each other ( Figure 3B). Specifically, we observed 96,640 SNVs, 914 insertions and 375 deletions in GB01; 48,050 SNVs, 992 insertions and 1621 deletions in GB02; and 33,864 SNVs, 233 insertions and 1446 deletions in GB03 ( Figure 3B). Within SNVs we find a similar percentage in the substitution of T > C, C > T, and C > G alleles among the three samples ( Figure 3C). GB02 has a higher percentage of T > G, GB03 of T > A and GB01 of C > A compared to the other tumor samples ( Figure 3C).

Somatic Mutations
Somatic mutations were called using Mutect2 [19]. We started by focusing on the genes known to be altered in GB IDH1-WT, referring to two databases: my cancer genome-MCG, mycancergenome.org and American Brain Tumor Association (ABTA), abta.org. We found 30 mutated genes in the three GB samples from a total of 264 different variants. The oncoplot in Figure 3D represents each of the 30 genes: the type(s) of gene mutation present in each of the 27 intra-tumor regions (colored square, middle section of the graph), the percentage of tumor regions out of the total (27) reporting the specific gene mutation(s) (bar graph on the right) and the total number of variants, considering all 30 genes, for each specific tumor region (bar graph at the top). Some genes (APC, ATM, ATR, and FAT1) are mutated in all 27 regions of the three GB samples ( Figure 3D). In particular, ATM, ATR, and FAT1 show multi-hit mutations (>2 types of mutations on the same gene, black squares, Figure 3D). EGFR carries mutations in all tumor regions of GB02 and GB03 while only three regions of GB01 have EGFR mutations ( Figure 3D). ERBB3 has mutations in all tumor regions of GB01 and GB03, whereas GB02 has no ERBB3 mutations in any of its regions ( Figure 3D). All regions of GB01 and GB02 share the mutated BRIP1 gene, whereas all regions of GB03 demonstrate no mutations in BRIP1 ( Figure 3D). GB01 has mutations in PTEN and IK in all nine regions, as well as GB03 has all regions with mutations in MSH6, PIK3CA and TP53 ( Figure 3D). There are also several genes that carry mutations in only one or two of the nine regions within the same tumor, such as ARID2, BAD, BCL2, BRAF, CBL, CDKN2B, FBXW7, KMT2D, MDM2, MET, MYC, PIK3CA, PTEN, SF3B1, SOX2, TERT and TP53 ( Figure 3D).

Rare Variants
In order to differentiate our somatic mutations as common or rare variants in the population, we set a minor allele frequency (MAF) threshold of less than 0.01 to define rare variants. We discovered a total of 17,068, 4255, and 3039 rare variants for GB01, GB02, and GB03 tumors, respectively.
The variant classification plot in Figure 4A shows the number and classification of variants in each tumor region, providing percentages of missense, nonsense, splice-site, and frame shift mutations (both insertion and deletion) out of the total number of mutations. population, we set a minor allele frequency (MAF) threshold of less than 0.01 to define rare variants. We discovered a total of 17,068, 4255, and 3039 rare variants for GB01, GB02, and GB03 tumors, respectively.
The variant classification plot in Figure 4A shows the number and classification of variants in each tumor region, providing percentages of missense, nonsense, splice-site, and frame shift mutations (both insertion and deletion) out of the total number of mutations.  To assess the similarities and differences between each of the intra-tumoral regions, as well as between different tumors, we performed hierarchical clustering on the rare mutations ( Figure 4C). The hierarchical analysis generates three different clustering profiles, one per tumor, based on the dissimilarity of the 9 intra-tumoral portions due to the presence/absence of rare tumor mutations. The length of the vertical lines in each dendogram is proportional to the differences in rare mutations among individual tumor regions or among groups of regions. Tumor GB01 generates a dendrogram different from the other two tumors. GB01 is divided into two main groups, GB01-06/07/09 on one side and GB01-01-05/08 on the other side. In the second group GB01-08 differs from the other 4 regions (01-05) which are more similar to each other. In contrast, the intra-tumoral regions of GB02 appear to have profound dissimilarities with each other, and the cluster shows a progressive subdivision into gradually smaller subgroups with increasingly similar regions. Furthermore, in GB03, all nine intra-tumoral portions show great dissimilarity from each other. With the exception of regions GB03-02 and G03-07, which appear to have more similarity, the other regions show more independent profiles.

Phylogenetic and Clonal Evolution
To further select potential disease-related genes, an additional filtering step was added to rare variants. Rare mutations in genes already described and annotated in COSMIC in GB tumors were considered for these two analyses.
We used phylogenetic tree inference (PTI) [31] to infer the rooted phylogenetic tree between different regions of the same tumor. Once the mutations are identified for all samples and the number of shared mutations are defined, the system forms the trunk of the tree. Then, through the unique mutations of individual regions and those shared by two or more tumor regions, the PTI finds the optimal branch division until all tumor regions reach the leaf nodes ( Figure 6A and Table S2).
The trees are built on the presence/absence of mutations; leaves from the same branch share some mutations but differ in others. Thus, each leaf (tumor region) will be unique due to the presence of a specific combination of mutations. By following the early (mutations shared between regions, present in the trunk/branches) or late (mutation at the level of the leaves) occurrence, it will be possible to reconstruct the mutational history of these tumors. Each trunk/branch/leaf shown in Figure 6A has a length dictated by the number of mutations present. The distinction between intra-tumoral regions is dictated only by the mutations in the genes that appear on the phylogenetic tree structure in Figure 6A. Among these mutated genes, we also found driver genes (in bold in Figure 6A): KMT2C, ATM, ATRX, and TSC2 in GB01; LZTR1, TGFBR2, and POLE in GB02; and PDGFRA in GB03. Some driver gene mutations are early mutations because they are present in branches shared by more than one tumor region and thus arose earlier than others exclusive to the regions. Among these, we find ATRX and TP53 in GB01 and, LZTR1 and TGFBR2 in GB02.
To gain insight into the evolutionary status of these three GB cases and to reconstruct multi-sample tumor phylogenies and decompose tumor subclones, we used LICHeE (Lineage Inference for Cancer Heterogeneity and Evolution). This system uses the variant allele frequencies of rare somatic mutations to reconstruct multi-sample cell lineage trees and infer the subclonal composition of the samples. The evolution trees and subclonal composition of GB tumors are shown in Figure 6B. The complete list of representative mutations in each clonal subpopulation is in Table S3.
For all three GBs, each tumor region is characterized by at least one clonal subpopulation exclusive to that region (colored rectangles within the square representing the tree leaf). In some regions (02, 04, 06 and 09 in GB01 and 04 in GB02) we note that two clonal subpopulations coexist, one within the other. As an example, in the case of GB01-06 ( Figure 6B

Phylogenetic and Clonal Evolution
To further select potential disease-related genes, an additional filtering step was added to rare variants. Rare mutations in genes already described and annotated in COSMIC in GB tumors were considered for these two analyses.
We used phylogenetic tree inference (PTI) [31] to infer the rooted phylogenetic tree between different regions of the same tumor. Once the mutations are identified for all samples and the number of shared mutations are defined, the system forms the trunk of the tree. Then, through the unique mutations of individual regions and those shared by two or more tumor regions, the PTI finds the optimal branch division until all tumor regions reach the leaf nodes ( Figure 6A and Table S2). The trees are built on the presence/absence of mutations; leaves from the same branch share some mutations but differ in others. Thus, each leaf (tumor region) will be unique due to the presence of a specific combination of mutations. By following the early (mutations shared between regions, present in the trunk/branches) or late (mutation at the level of the leaves) occurrence, it will be possible to reconstruct the mutational history of these tumors. Each trunk/branch/leaf shown in Figure 6A has a length dictated by the number of mutations present. The distinction between intra-tumoral regions is dictated only by the mutations in the genes that appear on the phylogenetic tree structure in Figure 6A. Among these mutated genes, we also found driver genes (in bold in Figure 6A): KMT2C, ATM, ATRX, and TSC2 in GB01; LZTR1, TGFBR2, and POLE in GB02; and PDGFRA in Figure 6. Tumor phylogeny and clonal evolution. (A) Phylogenetic trees. For each sample, a rooted tree was created whose leaf nodes are tumor regions. The length of the branches is equal to the number of mutations. Genes annotated on the tree are those with which the system has distinguished (through the absence or presence of specific mutations) and divided the different leaf branches representing the 9 tumor portions. Genes in bold are driver genes. (B) Clonal evolution lineage tree and sample composition. The lineage tree was built based on the constraint network using Lineage Inference for Cancer Heterogeneity and Evolution (LICHeE). Each node (circle) represents a subpopulation of GB cells. All nodes arose from a single hypothetical clone called Germline (GL), representing the genetic architecture of normal tissue from the same patient. Numbers within the circles indicate the number of nucleotide variants shared by the cluster; numbers outside the circles show the average variant allele frequency (VAF) of the variants in each cluster. Squares represent each individual tumor region, with colored rectangles indicating the cell fraction represented by the clonal cluster.

Discussion
GB is the most frequent brain tumor in adults and is characterized by an invariably fatal prognosis [32]. With optimal treatment, the median survival of patients is just over 1 year and only 18% survive two years [33]. GB is considered one of the most feared of all human diseases, both because of the lack of a cure and the associated loss of cognitive function as part of the disease process [27]. Currently, there are few biomarkers of favorable prognosis and consequently few therapies that strongly influence disease outcome. This is primarily due to the fact that the extreme heterogeneity of these tumors makes therapies increasingly challenging.
Along with inter-tumor heterogeneity, intra-tumor heterogeneity represents a pivotal area of investigation of GB. Several genomic studies have highlighted this heterogeneity by demonstrating a difference in somatic alterations, expression subtypes, and epigenetic modifications between different GBs and within the same tumor [13,27,[34][35][36][37][38][39][40]. Other studies have attempted to go even further by specifically examining regional heterogeneity between multiple sectors of the same primary GB, sometimes also coupled to relapsed tumor, through whole exome and transcriptome sequencing [2,3,[41][42][43]. Notably, multiregional exome sequencing studies of primary GB have been conducted on up to four tumor portions [2,41,42].
In this study, we aimed to further improve the analysis by going into even more detail in tumor heterogeneity. Therefore, we increased the number of tumor portions to be analyzed by choosing nine regions and analyzing them by whole exome sequencing. We also sought to ensure cohort homogeneity by selecting only GB IDH-WT without 1p/19q codeletion and patients of the same sex (women), with similar age and clinical course. Through a multiregional study of tumors from post-surgical samples, we aim to increase our understanding of tumor evolution and highlight the importance of tumor sampling from spatially distinct regions to avoid misinterpretation of genomic data from a single sample collection. Despite the well-known difficulties in collecting multiple samples of a glioblastoma in surgical procedures, the goal of deepening the understanding of tumor heterogeneity is also useful in devising potential new solutions to address the complexities of tumor heterogeneity in the face of the reality of therapeutic decisions based on limited access to tumor tissue.
All three of our samples are characterized by an extremely high level of DNA copy number alteration, as shown in Figure 2A,B. Several copy number variations, reported within and among the three tumors, commonly characterize primary IDH-WT GB, such as: deletion of chromosome 9 and 10 and amplification of chromosome 7 [44,45]. However, there are intra-tumoral differences in gain and loss at loci 9q, 7p, and 7q. Other observed alterations are deletion of 13q and amplification of 19q, both previously described in GB [46], both with intra-tumoral differences. In many other chromosomal loci (15q, 19p, 20q, 22q, and entire chromosomes 3, 8, 16, and 17), alterations, both deletions and amplifications, are present and are reflected quite unevenly across tumor regions.
A molecular classification of GB performed by The Cancer Genome Atlas (TCGA) identified four molecular subgroups with presumed prognostic significance [27]. The subgroups of GB described by TCGA, namely classical, neural, proneural, and mesenchymal, were identified and classified based on transcriptional profiles and supplemented with mutational and DNA copy number alteration profiles [27]. We relied on this classification system to define the molecular subtype of our tumor regions by correlating changes in DNA copy number. As displayed in Figure 2C,D, for all three samples, within the same tumor, more than one molecular subtype was correlated and identified. This confirms the findings of Sottoriva and colleagues [3], who through gene expression analysis, have observed the coexistence of different molecular subtypes within single GB tumors. From a mutational perspective, we first looked for somatic mutations in genes known to be mutated in GB. Although some genes (APC, ATM, ATR, FAT1, EGFR, ERBB3, BRIP1, PTEN, IK, MSH6, PIK3CA and TP53) were found to be mutated in all regions of the same tumor, the type of alteration is not always the same in all nine regions. In addition, there are also several genes that carry mutations in only one or two of the nine regions of the same tumor (ARID2, BAD, BCL2, BRAF, CBL, CDKN2B, FBXW7, KMT2D, MDM2, MET, MYC, PIK3CA, PTEN, SF3B1, SOX2, TERT and TP53). All of these variants were then ranked for their predicted pathogenicity, and a total of 59 variants in 18 genes were reported as pathogenic or probably pathogenic. None of these mutations are shared by all regions of the three tumors. Of these 59 variants found, 86% are mutations that affect only a single region in the nine within the tumor, meaning that the mutation is present in only 11% of the cells of the entire tumor. Twenty-six variants are known COSMIC mutations, and among them, two variants have already been described in GB [47,48]: EGFR Arg108Lys present in all regions of GB02 and KMT2D Arg2830Ter in region 07 of GB01.
We then further performed several analyses selecting the rare somatic mutations found in our samples. We first performed an unsupervised analysis to group the nine regions of each tumor according to the mutations they carried. In all three tumors, the tumor regions appear to be very dissimilar to each other, except for GB01, where we note four regions (01-05) that are more closely related to each other. The regions of tumors GB02 and GB03 cluster in distinct ways; in GB02, the cluster shows a gradual division into smaller and smaller groups with increasingly similar regions to each other; in GB03, all nine intra-tumoral portions show great dissimilarity to each other. In each tumor, however, the regions show highly independent profiles.
Next, we investigated how rare mutations might affect the 10 canonical cancer molecular pathways: cell cycle, Hippo, MYC, NOTCH, NRF2, PI-3-Kinase/AKT, RTK-RAS, TGFβ signaling, p53, and β-catenin/WNT. Within each of these pathways, we looked in detail at the mutated genes and alterations that occur in each tumor region. As reported in Figure 5B, there are both suppressor genes and tumor oncogenes mutated in all nine/nine parts of the tumor, such as EP300 and NCOR2 (NOTCH), PTEN (PI3K), EGFR (RTK-RAS) and FAT1 (Hippo). Other genes belonging to these oncogenic pathways were found mutated in all three tumors in at least one of their regions: DCHS1 and TAOK2 (Hippo), MET, JAK2, DAB2IP and RASGRF1 (RTK-RAS), MTOR and INPP4B (PI3K), NCOR1 and MAML3 (NOTCH), TGFBR2 (TGF-β) and WNT7a (WNT). Most importantly, many tumor suppressor genes or oncogenes are mutated in only one of the nine tumor portions.
From the list of mutated genes and pathways that emerges from our analyses, it is important to underline how our results confirm data already reported in the literature. In fact, it is well known as genetic alterations in GB typically deregulate pathways involving the tumor suppressors p53 (87%), RB (78%), and receptor-tyrosine kinase (RTK)/RAS/PI3K (88%) [49]. The PI3K/Akt/mTOR signaling pathway after EGFR activation is one of the most significant signaling pathways in cancer cells, and as clinical research has revealed, mutations in both EGFR and PTEN would lead to continuous activation of the signaling pathway, thus contributing to tumorigenesis and resistance to therapy [50]. INPP4B has been identified as a tumor suppressor often associated with PTEN in several cancers, and loss of INPP4B protein has also been correlated with reduced patient survival [51]. Mutations in MET, as well as the dysregulation of other regulators of crosstalk with MET signaling pathways in glioma, have also been previously identified and associated with proliferation, survival, migration, invasion, angiogenesis, stem cell characteristics, therapeutic resistance, and glioblastoma recurrence [52]. JAK2, in the RTK-RAS pathway has also garnered significant interest in the past as a key driver of cancer cell survival, proliferation and invasion in GB [53].
The Notch pathway has already been described in glioblastoma tumorigenesis [54] and correlated with GB development and proliferation as well as with its prognosis [55]. Specifically, EP300 and NCOR2 expression levels were correlated with survival and therapeutic response in GB [56]. In cancer, EP300 has been previously reported to target both mutations and structural alterations [57], which could be an important contributor to malignant transformation [58][59][60]. Growing evidence has demonstrated the role of Hippo signaling in cancer biology, and its alteration has been closely linked to tumorigenesis in many different cancers, also playing a key role in glioma development, including cell proliferation, apoptosis, and invasion [61,62]. FAT1 mutations in GB have already been reported, and studies strongly suggest that members of the FAT gene family are important players in cancer development [63][64][65].
As a concluding analysis, we aimed to reconstruct phylogenetic and clonal evolution trees ( Figure 6A,B), either by accounting for the presence/absence of a mutation in different tumor regions and by decomposing GB into tumor subclones. We performed this analysis using rare mutations of genes already described and annotated in COSMIC in GB tumors.
Although we find early mutations, including those of driver genes (ATRX, TP53, LZTR1 and TGFBR2), present in branches shared by more than one tumor region, the results of this analysis demonstrate a high degree of intra-tumoral heterogeneity in the three patients, with private leaf node mutations dominating the mutational landscape. By analyzing the allele frequencies of rare somatic mutation variants to reconstruct multi-sample cell lineage trees and infer the subclonal composition of the samples, we found that each tumor region is characterized by at least one clonal subpopulation exclusive to that region. Indeed, some regions show two clonal subpopulations, one within the other, while others are characterized by the co-presence of two distinct subclones.
In this work, we conducted most of the mutational investigations focusing on altered genes already described in GB; however, the majority of the variations found have not yet been annotated. The large number of mutations identified in our samples may therefore also include mutations that have not been functionally characterized in the current literature. Further future studies are certainly needed to investigate their functional roles.
The high level of heterogeneity of these three GBs is shown both at the mutational level, in somatic mutations, even rare ones, in DNA copies and also in the chronology of mutational evolution and in the subclonal composition of individual tumor regions. As never reported before, by increasing the regional subdivision of the tumor, up to 9 portions, the intra-tumor heterogeneity remains highly represented. All this suggests that most mutations are transient mutations, and that each intra-tumoral region is characterized by its own alterations. This seems in agreement with what was previously discussed by Sottoriva [3]. In particular, we can speculate that there are multiple clones with different fitness co-existing within GB.
The study and understanding of these patterns of heterogeneity could then be used to stratify individual patients and select an appropriate therapeutic strategy. To do this, however, multiple sampling is strictly necessary, as it is highly unlikely that a single biopsy could represent the complete set of mutations present underestimating the landscape of all alterations present. There are also hypotheses for which characterization of GB intra-tumor heterogeneity using single cell sequencing technology would be necessary [66]. Heterogeneity among tumor cells does not arise only as a result of molecular and genetic changes but is also influenced by different microenvironments within the tumor and reversible changes in cellular properties [67,68]. GB is composed of an interactive network of neoplastic and non-neoplastic cells, characterized by cell-cell interactions and connections of cellular compartments, creating a complex microenvironment that constantly gives signals, activating cell migration and ultimately developing permissive niches that promote tumor cell survival and proliferation [69,70] and influence therapeutic response [70][71][72]. Correlating key histological features of GB, different tumor microenvironments, and intracranial tumor locations to mutational patterns would allow for a greater and more complete understanding of tumor heterogeneity.
As Spiteri [43] also suggests, not only primary tumor should be studied, but residual disease in relapsed GB cases is equally important to be investigated to understand how treatment-resistant disease develops and to further investigate tumor evolution. It would also be important and useful to develop improvements during tumor surgery such as harvesting the subventricular zone and infiltrative margin, which is currently not possible due to the lack of reliable tumor cell markers [43].

Conclusions
In conclusion, our study is the first to use multiregional WES in such a large number of portions (nine tumor regions per tumor) to study the tumor heterogeneity. Although the sample size of the current study is relatively small (three primary IDH-WT GBs), our data provided additional insights into intra-tumoral heterogeneity, painting a still incomplete picture that warrants further investigation. Despite therapeutic opportunities for GB remain limited, any continued effort to study and understand tumor heterogeneity may help find a cure for this devastating disease.

Supplementary Materials:
The following are available online at https://www.mdpi.com/article/ 10.3390/cancers13092044/s1. Table S1: Sequencing run information for each tumor region. Mapped reads, total number and percentage of reads mapped inside target regions; Coverage, mean of the coverage depth; Quality, mean mapping quality of the mapped reads; CG percentage, GC percentage in the mapped reads, Table S2: Binary matrix output from PTI analysis for GB01, GB02 and GB03. The first two columns record unique gene and mutation IDs. The middle columns record whether a single tumor portion (leaf) has a certain mutation or two/more tumor regions (trunk or branches) share it. 1 = mutation, 0 = no mutation. The last column notes whether the genes involved in the mutation are genes drivers, Table S3: Representative mutations in each subclonal population. Sample, tumor; Region, tumor region; Size, number of mutations in that specific subpopulation; VAF, variant allele frequency of representative mutations; Gene, mutated genes; Mutation, mutation detail; AA change, amino acid change.

Conflicts of Interest:
The authors declare no conflict of interest.