Characterizing the Genomic Profile in High-Grade Gliomas: From Tumor Core to Peritumoral Brain Zone, Passing through Glioma-Derived Tumorspheres

Simple Summary The genomic landscape of the stem cell compartment and peritumoral brain zone of glioblastoma is still incomplete. The key role of the stem component in tumor maintenance and progression, as well as in drug-resistance spreading, has already been demonstrated. In recent years, the importance of the marginal area of the neoplasm has been considered since this is where the tumor recurrences appear in up to 90% of cases. In this study, we carried out a 360° genomic profile analysis of the different glioblastoma components in order to understand how they are characterized and how they work. Studying their genomic constitution is the starting point necessary to finally develop target treatments associable with the standard ones, as a new hope for glioblastoma patients. Abstract Glioblastoma is an extremely heterogeneous disease. Treatment failure and tumor recurrence primarily reflect the presence in the tumor core (TC) of the glioma stem cells (GSCs), and secondly the contribution, still to be defined, of the peritumoral brain zone (PBZ). Using the array-CGH platform, we deepened the genomic knowledge about the different components of GBM and we identified new specific biomarkers useful for new therapies. We firstly investigated the genomic profile of 20 TCs of GBM; then, for 14 cases and 7 cases, respectively, we compared these genomic profiles with those of the related GSC cultures and PBZ biopsies. The analysis on 20 TCs confirmed the intertumoral heterogeneity and a high percentage of copy number alterations (CNAs) in GBM canonical pathways. Comparing the genomic profiles of 14 TC-GSC pairs, we evidenced a robust similarity among the two samples of each patient. The shared imbalanced genes are related to the development and progression of cancer and in metabolic pathways, as shown by bioinformatic analysis using DAVID. Finally, the comparison between 7 TC-PBZ pairs leads to the identification of PBZ-unique alterations that require further investigation.


Introduction
Glioblastoma (GBM) is the most common and malignant primary brain tumor, characterized by rapid progression, invasion, high genomic instability, intense angiogenesis and resistance to therapies [1]. Despite aggressive standard treatments, the prognosis remains extremely poor with a mean survival of 20.9 months [2]. The failure of current therapies is mainly due to the striking inter-and intratumoral heterogeneity of the disease [3][4][5], supported by the presence, within the tumor mass, of cells with stem-like properties, called glioma stem cells (GSCs) [6,7]. To complicate this picture, there is also the presence, even if there is still a role to be defined, of an area named the peritumoral brain zone (PBZ), at the margin of the tumor central core [8].
Although the original identification of GSCs dates back more than a decade, the purification and characterization of GSCs remain challenging. Since they play important roles in mediating therapeutic resistance through supporting radio-resistance, chemo-resistance, angiogenesis, invasion and recurrence, we need a deeper understanding of how to selectively target and ablate these tumor-initiating and -propagating populations [9][10][11][12][13]. The most compelling reason to study glioma biology with GSCs is the fact that they have been shown to be very tumorigenic in vivo and form diffuse and invasive tumors that are highly resistant to conventional treatments, indicative of actual patient disease in clinic [14,15].
The GBM-PBZ is a region that radiologically and macroscopically resembles normal brain tissue, but with a particular cellular content, which consists of infiltrating tumor cells, reactive astrocytes, inflammatory cells and other stromal cells [8]. Furthermore, since in 90% of cases tumor recurrence occurs at the margin of the surgical cavity, even after a complete tumor resection and chemo-radiotherapies, the glioma microenvironment seems to be a critical regulator of tumor progression [16]. Therefore, a better understanding of this area is crucial to unravel the mechanisms underlying the GBM relapse and to develop new therapeutic approaches [8,17].
At the DNA level, GBMs are usually characterized by high levels of genomic instability with high rates of copy number alterations (CNAs), easily identifiable by array comparative genomic hybridization (array-CGH) [9,[18][19][20]. Thanks to these studies, frequently amplified genes, such as EGFR, MET, PDGFRA, MDM2, PIK3CA, CDK4 and CDK6, and deleted genes such as CDKN2A/B, PTEN and RB1, have been highlighted [21,22]. However, the constant improvement in genetic characterization of GBMs is still failing to be translated to clinical practice, suggesting that other discovery paradigms should be examined.
Considering the importance of CNA data, in this work we performed an in depth study using array-CGH in order to outline the genomic profiles of 20 tumors and 7 peritumoral biopsies; in addition, we compared the genomic profiles of 14 tumor biopsies with their derived tumorspheres to identify new specific biomarkers useful for new therapies.

Study Population
The patient population of this study consists of 20 adults of both sexes diagnosed with high-grade gliomas. The study was approved by the ethic committee "Comitato Etico Monza e Brianza" (study number: 0031436-GLIODRUG-V, approved on 3 January 2020). Patients undergoing a craniotomy for a high-grade glioma were enrolled between January 2020 and September 2021 by the Neurosurgery Unit of the San Gerardo Hospital (Monza Brianza, Italy) after informed consent was signed. The criteria for histological analyses were based on the recommendations of the 2021 WHO classification of CNS tumors [23] and samples from patients without a confirmed high-grade glioma were excluded from the study. Demographic and clinical data are reported in Table 1. Table 1. Clinicopathological characteristics and immunomolecular phenotype of glioma patients (GP) enrolled in the study. M: male; F: female; GBM: glioblastoma IDH-wildtype; ASG4: astrocytoma IDH-mutant, grade 4; ASG3: astrocytoma IDH-mutant, grade 3; wt: wild-type; mut: mutated; +: <30% of positive cells; -: no staining; met: methylated; follow-up period expressed in months; DOD: patient dead of disease; R: relapsed patient, follow up ongoing; PF: progression free patient, follow up ongoing; n.a.: no information available.

Biopsy Collection
After surgical removal of the tumor with perilesional dissection [24], the surgical cavity was checked with an intraoperative ultrasound (BK Medical, Herlev, Denmark) for tumor remnants. All patients received a resection of at least 95% of contrast enhancing tumor. The majority of the tumor was sent for formal histological diagnosis. When the surgical bed was considered tumor free, a sample of the surgical cavity was taken. The non-neoplastic tumor margin (PBZ) was collected from what was considered far from functionally eloquent areas under neuronavigation (BrainLab, Munich, Germany) and ultrasound guidance (BK Medical, Denmark). Finally, in only 10 out of 20 patients, it was possible to safely collect PBZ samples. When PBZ was collected, a part of the specimen was sent for formal histology in order to check for the absence of tumor remnants. Biopsy samples were placed in cooled PBS (Euroclone S.p.A., Milan, Italy) containing 10% antibiotics (streptomycin/penicillin, Euroclone S.p.A., Milan, Italy) for the isolation of GSCs [25]; small pieces were also stored for DNA extraction.

Immuno-Molecular Analysis
Immuno-molecular analysis was performed according to routine diagnostic procedures. Immunohistochemical staining was performed on FFPE (4% formalin) sections of 1-µm thickness, according to the manufacturer's protocols, using the automated instrument Dako Omnis (Agilent Technologies, Santa Clara, CA, USA). All antibodies were purchased from Dako. For p53, only nuclear staining was considered positive. MGMT promoter methylation and IDH mutation analysis were performed as previously reported [26].

GSCs Isolation from Tumor Tissues
GSCs were isolated directly from the TC samples after surgery. Briefly, samples obtained from tumor cores were washed with PBS and placed in a Petri dish. Then, they were disaggregated mechanically and enzymatically with a 1X trypsin-EDTA solution (Euroclone S.p.A., Milan, Italy). The digested tissue was passed through a cell strainer (70 µm) and finally subjected to lysis of red blood cells. Lastly, single cell suspension was seeded in a complete neural stem cell (NSC) culture medium (see below) at a density of 40,000 cells/cm 2 .

Primary GSC Cultures Conditions
GSCs were cultured in a selective medium for NSC, composed of DMEM F-12 and Neurobasal 1:1, B-27 supplement without vitamin A (Life Technologies Italia, Milan, Italy), 2 mM L-glutamine (Euroclone S.p.A., Milan, Italy), 10 ng/mL recombinant human bFGF and 20 ng/mL recombinant human EGF (Miltenyi Biotec, Bergish Gladbach, Germany), 20 UI/mL penicillin and 20 µg/mL streptomycin (Euroclone S.p.A., Milan, Italy). After isolation, the medium was replaced every 3 days to remove stroma and red blood cells residues, catabolic products and to supply fresh nutrients. Debris and adherent death cells generally were eliminated after a couple of passages. The isolated cells propagate in culture as free-floating spheres defined as tumorspheres [6], which appeared in 15-20 days of culture after isolation. When tumorspheres reached an average size of 100 µm in diameter, the culture was ready to be passed and expanded. At each passage (P), tumorspheres were mechanically dissociated using a sterilized p200 pipette set at 180-200 µL and pipetting up and down 100-150 times to achieve a single-cell suspension.

Established Glioma Stem Cell Lines
Two established glioma stem cell lines were used as a positive control of stemness. G166 and G179 cell lines were kindly provided by Professor A. Smith of the Wellcome Trust Medical Research Council Stem Cell Institute, University of Cambridge, Cambridge (UK). These cell lines were extensively characterized by Pollard and Baronchelli [9,14]. The established stem cell lines were cultured as the GSC primary cultures; a G166 line grows as floating spheres, otherwise G179 grows in semi-adhesion.

Clonal Assay
Mechanically dissociated tumorspheres were seeded into 96-well plates at a density of 10 cells per mL in culture medium. Colony formation was scored 7-10 days after initial seeding. The self-renewal efficiency or the percentage of cells that formed spheres was determined by the following formula: n = Y X · 100; where Y is the number of wells in which one tumorsphere is developed from a single cell and X is the number of wells in which a single cell was present [27]. Wells containing either none or more than one cell were excluded from the analysis.

Differentiation Assay
Mechanically dissociated tumorspheres were seeded at a density of 1 × 10 5 cells/well into 6-well plates, with a coverslip on the bottom of each well, and into culture medium permissive for differentiation without EGF and with 5% FBS (Euroclone S.p.A., Milan, Italy). After 4 days of culture, the medium was replaced with fresh medium with 5% FBS, without growth factors. Under these conditions, the detection of the three neural lineages was evidenced at 7 days after plating by immunofluorescence.

DNA Extraction and Purification
DNA was extracted from GSC primary culture pellets (between P4 and P6), from tumor and peritumor biopsies and from patients' blood (used as reference) using the automatic extractor iPrep TM (Thermo Fisher Scientific, Waltham, MA, USA) and using kits supplied with the instrument: iPrep tissue, for DNA extraction from cell pellet, TC and PBZ; and iPrep whole blood, for DNA extraction from peripheral patients' blood. Then, DNA was purified using Genomic DNA clean & concentration kit (Zymo Research, Irvine, CA, USA) based on several washings and elution on column in order to obtain DNA ultra-pure. The concentration and the purity of the extracted DNA were determined by measuring the absorbance (A260/280) of the sample with NanoDrop ND-1000 Spectrophotometer (Thermo Fisher Scientific, Waltham, MA, USA). In some cases, the extracted DNA from primary culture was not sufficient to perform the analysis and it was amplified using the GenomePlex Whole Genome Amplification (WGA) Kit (Sigma-Aldrich, St. Louis, MI, USA), according to the manufacturer's instructions. Amplified DNA was tested for purity and concentration as above.

Array-CGH
Array-CGH analysis was performed using 60-mer oligonucleotide probe technology (SurePrint G3 Human CGH 8 × 60 K, Agilent Technologies, Santa Clara, CA, USA), according to the manufacturer's instructions. Agilent Feature Extraction was exploited to generate raw data, which were further analyzed using Cytogenomics 5.1 with the ADAM-2 algorithm (Agilent Technologies, Santa Clara, CA, USA). A minimum of three consecutive probes/regions was considered as a filter. The threshold for genomic deletion is x = −1; the threshold for genomic gain is x = +0.58. The estimated percentage of mosaicism was calculated using the formula reported in [9]. Notably, in a mosaic scenario, the threshold is between −1 and 0 for deletions and between 0 and +0.58 for duplications. Amplifications and homozygous deletions are considered with threshold >+2 and <−1, respectively.

Bioinformatics Analysis
The Database for Annotation, Visualization and Integrated Discovery (DAVID), v 6.8 https://david.ncifcrf.gov/summary.jsp/ (accessed on 7 June 2021) [28,29], was used to analyze the lists of genes included in CNAs shared in at least 3 samples. The chart function was used to identify pathways in which genes in gain and in loss are involved. The clustering function was used to cluster the pathways found in groups with their own enrichment score.

Statistical Analysis
A chi-square test was used to compare data relating to patient-derived primary GSC cultures and those obtained from control GSC cell lines. Pearson correlation was used to compare genetic alteration profiles in matched pairs of TCs and patient derived GSCs and TCs and related PBZ samples. The statistical significance for each pair of correlations was calculated by consulting the Table of Critical Values for Pearson's R. The level of significance for a two-tailed test was set to α = 0.05. The EASE Score, a modified Fisher's Exact Test, is used by the DAVID bioinformatics program. p < 0.05 was considered statistically significant.
To maximize the ability to correlate our results with clinical features, as the patients were enrolled from January 2020 to September 2021, they were further classified into 3 categories based on the actual follow-up, if available: dead of disease (DOD); progression free (PF), in case of alive patients lacking recurrence; and alive with recurrence (R).

Isolation, Expansion and Characterization of GSC Primary Cultures
We were able to generate and expand 15 primary cell cultures with an isolation efficiency of 75%. About 15-20 days after the isolation, we observed the formation of tumorspheres in 15 out of 20 patient-derived cultures ( Figure 1a). However, as the sole formation of tumorspheres does not provide, per se, a demonstration of the presence of GSCs, we evaluated the self-renewal efficiency and multipotency. We performed clonal assay experiments on our GSC primary cultures and their average self-renewal efficiency (78%) was not statistically different from those of two established GSC lines (83%) (chisquare test) (Figure 1b) [14]. The stem nature of our GSC primary cultures was further demonstrated by evaluating the expression of specific markers of stemness (CD133 and nestin). In addition, we also assessed their ability to differentiate in neural lineages by evaluating GFAP, β III tubulin and MBP expression. All the investigated markers were expressed in our GSC primary cultures similarly to the control cell lines (Figure 1c,d).

Genomic Profiles of Tumor Biopsies Confirmed Canonical Alterations of GBM and Inter-Tumor Heterogeneity
We performed a genomic characterization of 20 tumor core biopsies (TCs) by array-CGH. A total of 78 copy number alterations (CNAs) disrupting the three canonical pathways involved in gliomagenesis, p53, Rb and PIK3KC, were found in our samples ( Table 2). The mean mutation burden of canonical CNAs was 3.9 per sample, ranging from zero (TC6 and TC15) to even 7-8 (TC13, TC18 and TC27). The most altered pathway, with 35 alterations (~45% of the total canonical CNAs), was the Rb pathway: 17 out of 20 tumor biopsies (85%) had imbalances in this pathway. The 9p21.3 locus (CDKN2A/B) was lost in 12 tumors; 11 biopsies had gains in 7q21.2 (CDK6); 12q14.1 locus (CDK4) was affected in 6 samples (five gains and one loss); 4 samples led to a loss in the 13q14.2 locus (RB1); finally, locus 12p13.32 (CCND2) was altered in two samples. The second pathway with the greatest number of alterations (34,~44% of the total canonical CNAs) was the PIK3KC pathway: 17 out of 20 tumors (85%) were found to be altered in its genes. The EGFR locus (7p11.2) had gains in 14 biopsies; the PTEN locus (10q23.31) was affected in 14 samples (13 losses and 1 gain); 3 samples had a gain in the PDGFRA locus (4q12), while 2 tumors lost in the NF1 locus (17q11.2); one tumor had a gain in the 5q13.1 locus (PIK3R1), but none had alterations in the 3q26.32 locus (PIK3CA). The p53 pathway was the least affected canonical pathway, but it was also the one with the fewest genes considered:~11% of the total CNAs in only 9 out of 20 tumors (45%). The MDM2 locus (12q15) had 4 imbalances (3 gain and 1 loss); the TP53 locus (17p13.1) has 2 losses and 1 gain, while MDM4 locus (1q32.1) had 1 loss and 1 gain. Definitely, the most altered canonical loci in our cohort were 7p11.2 (EGFR), 10q23.31 (PTEN), 9p21.3 (CDKN2A/B) and 7q21.2 (CDK6). assay experiments on our GSC primary cultures and their average self-renewal efficiency (78%) was not statistically different from those of two established GSC lines (83%) (chisquare test) (Figure 1b) [14]. The stem nature of our GSC primary cultures was further demonstrated by evaluating the expression of specific markers of stemness (CD133 and nestin). In addition, we also assessed their ability to differentiate in neural lineages by evaluating GFAP, β III tubulin and MBP expression. All the investigated markers were expressed in our GSC primary cultures similarly to the control cell lines (Figure 1c,d).  Moreover, considering the total burden of genomic alterations, not just in the canonical pathways, the total CNAs' load was 529, with a median of 23.5 and range 7-87.
Patients with a heavier burden of CNAs (≥23) relapsed (GP11 and GP18) or died of disease (GP7, GP10, GP12, GP23, GP24 and GP27) within a year from the diagnosis. However, patients with a lighter burden of CNAs (<23) relapsed (GP9) or died of disease (GP6, GP13 and GP14) one year after diagnosis, pointing out that other factors contribute to worsening of prognosis, in addition to CNAs (Tables 1 and 2).

Genomic Profiles of Tumor Biopsies and Matched Derived Tumorspheres Showed a Good Correlation
In this study, we performed a genomic analysis of 15 patient-derived GSCs in order to further demonstrate their stemness properties, i.e., the ability to contain the genomic aberrations typical of GBM and perpetuate them indefinitely. In one case (GSC15) no alterations were reported, probably due to a low enrichment of the stem components in the culture, so it was not included in further analysis (Table 3). Table 2. Summary of copy number alterations affecting the three main canonical pathways altered in high-grade gliomas. The analysis was conducted on patient tumor core biopsies (TCs). The canonical CNAs load (Canonical CNA) and the total burden of imbalances (Total CNA) per sample were also reported.  We compared the genomic profiles of 14 GSCs and their relative TCs. Shared imbalances were evaluated for their dimension (base pair length) and their representativeness (mosaicism percentage) ( Figure S1, Table S1). Furthermore, in order to estimate the similarity between the respective pairs of genomic profiles, we calculated the Pearson correlation metric. We divided the patients into three groups, based on the strength of correlation [30]: very strong correlation (Pearson value R ≥ 0.80), moderately strong correlation (0.40 ≤ R < 0.80); and low correlation (R < 0.40) ( Table 4). Table 4. Report of Pearson correlation values (R) for TC-GSC-genomic profiles comparison. GPs are ranked into three groups: ones with a very strong correlation between the pair of genomic profiles (R ≥ 0.80); ones with a moderately strong correlation (0.40 ≤ R < 0.80); ones with a low correlation (R < 0.40). The statistical significance for each pair of correlations was calculated consulting the Seven patients belong to the first group, with Pearson values from 0.80 to 0.94. Three out of 14 patients are part of the second group, with Pearson values from 0.55 to 0.79. In these two groups, shared CNAs concerned mainly the loss of genomic material; numerical imbalances were less represented in TCs, demonstrated by the presence of mosaicism in most cases, while in GSCs the alterations were more homogeneous, thanks to a probable clonal selection. Notably, the number of alterations increased in GSC cultures, with respect to the matched TCs, except GP10, GP12, GP17 and GP18 (Tables 2 and 3). Patient GP18 relapsed within one year from surgery, they had a Pearson value score of 0.81. Patients GP7, GP10, GP12, GP13, GP23 and GP24 died of disease within a year after diagnosis and their Pearson values were, respectively, 0.94, 0.9, 0.87, 0.93, 0.92 and 0.55. All p values of our correlation metric resulted <0.05, showing a significant and strong correlation between TC and GSC of each patient of these two groups.
Four out of 14 patients were part of the low correlation group, with Pearson values from 0.07 to 0.16. Furthermore, for this group shared CNAs concerned mainly the loss of genomic material, and the imbalances were more represented in GSC cultures than in biopsies. Patient GP6 was the only member of this group to die of disease within a year from the diagnosis, they had a Pearson value score of 0.07. The p value of our correlation metric resulted <0.05 only in patient GP20.

Comparing the Genetic Profiles of Matched TCs and GSCs of Different Patients to Find New Targets for Therapies
We used DAVID (Database for Annotation, Visualization and Integrated Discovery) in order to investigate the pathways of genes involved in "common" aberrations detected in matched TCs and GSCs and shared among at least three different patients. We presented the top ten pathways recognized by the KEGG pathway database and the top ten GO functions (biological process: BP, cellular component: CC, and molecular function: MF) with statistically significant results. Gained genes were mainly involved in biological processes such as mRNA splicing via spliceosome, leukocyte migration and endocytosis.
Cytological composition analysis showed that most parts of the genes were significantly involved in the composition of membrane, extracellular region and extracellular space. The molecular functions were mainly concentrated in receptor binding, carbohydrate binding and lipid binding. The KEGG-pathway showed that these genes were mainly involved in the pathways in cancer, the Rap1 signaling pathway, Huntington's disease and the cAMP signaling pathway (Table 5). Lost genes were mainly involved in biological processes such as oxidation-reduction and the lipid catabolic process. Most parts of these genes were significantly involved in the cytological composition of cytosol, nucleoplasm and mitochondrion. Their molecular functions were mainly involved in oxygen binding and cytoskeletal protein binding. The KEGG-pathway showed that these genes were mainly involved in metabolic pathways, biosynthesis of antibiotics, the WNT signaling pathway and carbon metabolism ( Table 6).    Additionally, we evaluated the gene enrichments. Regarding gained genes, 157 genes were enriched in KEGG pathways, 39 in MF, 33 in BP and none in CC. Then, we observed that 42 lost genes were enriched in BP, 24 in MF, 9 in CC and none in KEGG pathways. Subsequently, we generated Venn diagrams for pathways/GO functions, in order to show the shared genes, both in gain and in loss, by the four datasets. The diagram with the gained genes shows that 12, 4 and 3 genes were shared between two datasets: BP and KEGG, BP and MF, KEGG and MF, respectively (Figure 2a). Similarly, the diagram with the lost genes shows that 15 genes were shared between BP and MF, while only one gene was shared in CC and MF (Figure 2b). Furthermore, both diagrams show no gene shared by all datasets simultaneously.

Revelations from the Comparison between Genomic Profiles of Tumor Core Biopsies (TCs) and Peritumoral Brain Zone (PBZs) Samples
In 10 patients it was also possible to collect the non-neoplastic peritumor margin; therefore, we have also extended the genomic analysis on these biopsies. However, three PBZ samples (PBZ6, PBZ9 and PBZ20) did not report any CNA, so, they were not included in further analysis. We compared the genomic profiles of 7 TCs and the matched PBZs and we calculated the Pearson correlation metrics, dividing the patients into three groups, based on the strength of the correlation [31]: very strong correlation (Pearson value R ≥ 0.80), moderately strong correlation (0.40 ≤ R < 0.80); and low correlation (R < 0.40) ( Table  7).

Revelations from the Comparison between Genomic Profiles of Tumor Core Biopsies (TCs) and Peritumoral Brain Zone (PBZs) Samples
In 10 patients it was also possible to collect the non-neoplastic peritumor margin; therefore, we have also extended the genomic analysis on these biopsies. However, three PBZ samples (PBZ6, PBZ9 and PBZ20) did not report any CNA, so, they were not included in further analysis. We compared the genomic profiles of 7 TCs and the matched PBZs and we calculated the Pearson correlation metrics, dividing the patients into three groups, based on the strength of the correlation [31]: very strong correlation (Pearson value R ≥ 0.80), moderately strong correlation (0.40 ≤ R < 0.80); and low correlation (R < 0.40) ( Table 7). Table 7. Report of Pearson correlation values (R) for TC-PBZ-genomic profiles comparison. GPs are ranked into three groups: ones with a very strong correlation between the pair of genomic profiles (R ≥ 0.80); ones with a moderately strong correlation (0.40 ≤ R < 0.80); ones with a low correlation (R < 0.40). The statistical significance for each pair of correlations was calculated consulting the GP11 and GP27 had a great correlation, with Pearson values, respectively, of 0.94 and 0.91, both were statistically significant (p < 0.05). Patient GP11 relapsed within one year from surgery and GP27 died of disease. GBM canonical alterations already present in TCs were confirmed in the related PBZs, confirming a very good overlapping of the genomic profiles ( Figure S3, Table S2).
Two out of 7 patients, GP10 and GP17, belong to the second group, with Pearson values, respectively, of 0.62 and 0.42, both are statistically significant (p < 0.05). In this group, shared CNAs concerned both loss and gain of genomic material. Numerical imbalances were less represented in PBZs, demonstrated by the presence of mosaicism, confirming the infiltration of cancer cells in this area. Interestingly, patient GP10 died of disease within a year after diagnosis and is the patient with the highest number of shared CNAs between TC and PBZ ( Figure S3, Table S2).
Three patients are part of the low correlation group. GP24 died of disease within one year from surgery and had a Pearson value of 0,37. Notably, patient GP15 and GP22 had a negative correlation, with a Pearson value score, respectively, of −0,87 and −0,52 (p < 0.05); a negative correlation attests that the two genomic profiles not only are not overlapping, but they are exclusive and unique ( Figure S3, Table S2).
However, it is interesting to pay attention to some imbalances exclusive of PBZs and not shared with other samples of the same patient (Table S3). In this case, gains seem more numerous than losses. Noteworthy are two gained regions: 11p11.2, evidenced in PBZ of GP9 and GP22; and 16p13.3, evidenced in GP15 and GP27 (Table S3). Both regions include genes of interest for glioblastoma, such as EXT2, SSTR5, SSTR5-AS1, C1QTNF8 and CACNA1H. Furthermore, two regions in gain (1p34.2 and 2q14.2) were evidenced specifically only in two PBZs of two patients (GP10 and GP27, respectively).
Finally, in three patients (GP10, GP22 and GP24), we evidenced several CNAs shared between PBZ with their matched GSCs and TCs. Curiously, in some cases the alterations were identified only in the PBZ and in the matched GSCs, suggesting that the tumor cells infiltrated in the PBZ may be the GSCs (Table S4).

Discussion
GBM is the most common and fatal primary brain tumor with a median survival rate of only 15 months after the first diagnosis. The current standard of care for GBM, as proposed by Stupp in 2017 [2], consists of maximal safe surgical resection, radiotherapy with concomitant temozolomide, followed by adjuvant temozolomide and tumor-treating fields. Unfortunately, despite the treatment, about 70% of these tumors recur with de novo or acquired resistance, which leads to a low five-year survival rate [32][33][34]. The reasons for this high failure rate are different and closely interconnected. First of all, similarly to what happens in benign tumors such as meningiomas, surgery cannot be considered curative, as the presence of cortical and subcortical functional areas, combined with the widespread infiltrative pattern of tumor growth, makes it difficult to perform a complete surgical resection enclosing a wide area of peritumoral brain [35]. As a consequence, a variable number of invading tumor cells are invariably left behind in the peritumoral brain zone (PBZ), which is most often the site of recurrence [16,36]. Secondly, the huge variability of this type of tumor plays a paramount role: interpatient, intratumoral, functional and molecular heterogeneity, widely described in literature, and which is most likely supported and nourished by GSC subpopulations, makes this tumor very resistant to adjuvant treatments that should destroy the neoplastic cells left behind by the surgeon. In light of this, it is urgent to find more effective therapies that target specifically GCS subpopulations, and simultaneously, to develop in vitro models able to reliably recapitulate the original tumor and that can be used for preliminary rapid and cost-effective testing [37,38]. Establishing primary GSC primary cultures would provide a valuable and accurate model of the human tumor, would give insight into the origins of tumor heterogeneity, and finally, would direct towards the most effective therapies for the patient [39]. In our study, we have shown that a simple and not expensive protocol is valid to establish primary cultures enriched with GSCs in at least 70% of cases. Indeed, the characteristics of stemness of our primary cultures, such as the formation of perpetuating tumorspheres in different passages and the ability to differentiate in neural lineages, have been confirmed. One of the ways to understand if the isolated cells faithfully represent the GSC subpopulations is to evaluate the similarity of the genomic profiles with their tumor cores (TCs). Moreover, we have already shown that the tumorspheres had specific genomic profiles, which can be used as a specific tracer of these subpopulations [9].
In the last years, several studies focusing on genomic, transcriptomic and methylomic analyses of GBM were completed [9,20,[40][41][42][43]. In addition, Lemée J. et al. also characterized the peritumoral brain zone of GBM, examining its genomic, transcriptomic and proteomic profiles, but it does not include the GSC component [8]. In this work, for the first time to our knowledge, we obtained a complete view of the genomic profiles of GBM, studying, wherever possible, the tumor core, the relative GSCs and also including the PBZ, all derived from the same patient.
First of all, we obtained the genomic profile of 20 TCs, highlighting that our cohort almost completely recapitulates the information reported in literature. In fact, the most frequent copy number alterations were 7p11.2 amplifications (EGFR locus), 9p21.3 deletions (CDKN2A locus) and 10q23.31 deletions (PTEN locus), that are typical features of primary GBM, essential for gliomagenesis [44]. However, tumor-specific profiles in terms of CNAs distribution were also observed, in accordance with the well-known heterogeneity of GBM.
At the same time, we were able to isolate and expand GSC primary cultures from 15 TCs and we compared their genomic profiles with the matched TC. We showed that 14 out of 15 GSC primary cultures very faithfully reflected the genomic profile of their original TC, as evidenced by the Pearson correlation and by the high number of shared CNAs. We observed that losses are better maintained than gains, and that all the alterations identified in TCs become more homogeneous and represented in GSC cultures, thanks to in vitro clonal selection. These shared CNAs are enriched by genes mainly involved in the pathways of cancer, as evidenced by GO and KEGG enrichment analysis, reinforcing the hypothesis that they are important for the progression and maintenance of cancer cells. We evidenced that a very large number of genes involved in gains and shared between TCs and GSCs belong to the GO term cellular component, particularly extracellular space, extracellular region and membrane, indicating a strong involvement in intercellular signaling. Interestingly, several gained genes shared between at least two datasets were reported in the literature as overexpressed in high-grade gliomas, associated with poor prognosis, or in any case, involved in canonical pathways of GBM (EIF2AK1, FTL, SRC, GNG8, GNG11, GNG7, FZD9, GNB2, TGM2) [45][46][47][48][49]. Conversely, only one lost gene, AKR1C2, shared between at least two datasets, was reported as downregulated and associated with high-risk-score glioma in the REMBRANDT data set [50]. These genes could represent interesting markers for prognosis and for new therapies in the future (Table 8). Table 8. Summary of relevant genes shared at least between two datasets of the enrichment analysis with DAVID. The following is reported: the cytogenetic localization, the function (GeneCards: The Human Gene Database, https: //www.genecards.org/, accessed on 20 September 2021) [51], the association with GBM and possible new targeted therapy (DrugBank database; https://go.drugbank.com/, accessed on 20 September 2021) [52].

Gene
Cytogenetic Band Function In GBM Target Drug The protein encoded by this gene acts at the level of translation initiation to downregulate protein synthesis in response to stress.
It is over expressed in glioma tissue rather than in normal brain tissue. Its expression is not variable among different GBM subtype [45]. Its expression is enriched in high-grade glioma (HGG) and together with IDH1/2 wildtype it is significantly associated with an unfavorable prognosis of glioma patients [46].
This proto-oncogene is part of the EGFR/SRC/ERK pathway. It may play a role in the regulation of embryonic development and cell growth.
Its role in YTHDF2 phosphorylation leads to an increase proliferation, invasiveness and tumorigenesis in GBM [53].
Dasinatib (phase 2 clinical trials) GNG7 19p13.3 A member of the guanine nucleotide-binding proteins (G proteins) gamma family that is involved as a modulator or transducer in various transmembrane signaling systems.
It has been associated with GBM due to their genomic localization close to WNT5A and WNT10B genes whose pathway is known to be altered in GBM [47].
This gene encodes heterotrimeric guanine nucleotide-binding proteins (G proteins), which integrate signals between receptors and effector proteins, and are composed of an alpha, a beta and a gamma subunit.
It was observed that GBM patients with a high GNB2 level are associated with higher pathological grade, wild-type IDH and unmethylated MGMT promoter [48].

TGM2 20q11.23
Transglutaminases are enzymes that catalyze the crosslinking of proteins by epsilon-gamma glutamic lysine isopeptide bonds.
Its high expression in glioma tissues has been reported [49]. It is expressed in the brain and its levels are significantly higher in malignant astrocytomas than in low-grade astrocytomas [54].

AKR1C2
10p15.1 This gene encodes a member of the aldo/keto reductase superfamily.
It was reported as downregulated and associated with high-risk-score glioma [50].
In this study, we were able to draw a complete picture of genomic profiles of GBM in 10 patients, adding information from the PBZ. There are still very few data on this area in the literature, and they are generally based on gene expression or imaging analysis [55,56]. The larger sharing of genomic profiles would testify to a greater infiltration of tumor cells in healthy areas of the brain, which, after surgical removal of the tumor core, would remain in the site, with a high risk of recurrence formation. Our data confirmed a previous study, which reported that tumor cell infiltration was found in one third of PBZs, despite radiological and macroscopic analysis revealing normal brain tissues [8]. However, histopathological examination did not find tumor cell infiltration in our series of PBZ, even if the overlap between PBZ-TC genomic profiles would prove otherwise for more than 50% of cases. Interestingly, we evidenced two exclusive PBZ-CNAs that could identify a specific signature in this specific area, because it was not present in other samples, nor in other patients. The first includes the SCMH1 gene, associated with the Polycomb group (PcG) multiprotein complexes, required to maintain the transcriptionally repressive state of some genes [57]. The second region includes the GLI2 gene, which promotes cell proliferation and migration in glioma [58,59]. Other potentially interesting regions are two gains, in 11p11.2 and 16p13.3, shared by two PBZs and evidenced also in other TCs and/or GSCs from other patients. Several genes included in these regions could hide an important role in the progression of glioma, as reported by the literature. For example, the expression level of the EXT2 gene is increased in glioblastoma [60]; C1QTNF8 promotes temozolomide resistance [61]; CACNA1H promotes GBM cell proliferation and migration [62]. Finally, we noticed that the CNAs in gains are more represented in PBZ samples, compared to the losses. We have already highlighted this phenomenon in other types of tumor [63], associating it with a possible mechanism of endoreduplication-polyploidization [64]. In support of this, the karyotype of GSCs of GP22 is tetraploid (data not shown), reinforcing the importance of this mechanism in glioma tumorigenesis.

Conclusions
In this study we were able to draw a complete picture of genomic profiles of GBM, from TCs to PBZs, including GSCs. Although the numbers are too few to draw conclusions, it is curious to note how some patients have died (GP24 and GP10) and their TCs, GSCs and PBZs have very similar genomic profiles. We have confirmed some strong points of glioma tumorigenesis, but we have also contributed to the identification of new genes and new regions potentially useful for future therapeutic strategies. The high heterogeneity of GBM requires the enrollment of a high number of patients in order to achieve statistically significant results and to understand the complexity of the disease. The limitations of our study derive from the small number of enrolled patients and the limited follow up. Nevertheless, the study of genomic profiles is a good starting point for a comparative analysis among the different components of GBM, and this is a strong point of our work. In the future we will move to high throughput approaches, such as single-cell analysis, in order to have a more complete overview of the disease.

Institutional Review Board Statement:
The study was conducted according to the guidelines of the Declaration of Helsinki and approved by the Ethics Committee of "Comitato Etico della provincia Monza e Brianza" (study GLIODRUG-V, approved on 3 January 2020).

Informed Consent Statement:
Informed consent was obtained from all subjects involved in the study.

Data Availability Statement:
The data used to support the findings of this study are available from the corresponding author upon request.