Glioblastomas within the Subventricular Zone Are Region-Specific Enriched for Mesenchymal Transition Markers: An Intratumoral Gene Expression Analysis

Simple Summary Involvement of the subventricular zone (SVZ) in glioblastoma is associated with poor prognosis and is associated with specific tumor-biological characteristics. In this study, we demonstrate that patient-derived glioblastoma samples from within the SVZ region show increased (epithelial-)mesenchymal transition and angiogenesis/hypoxia signaling as compared to glioblastoma samples from the same patient from outside the SVZ. These results suggest that intratumoral alterations in oncogenic signaling could be mediated by the SVZ microenvironment. Our findings offer rationale for specific targeting of the SVZ in the development of glioblastoma therapy. Abstract Background: Involvement of the subventricular zone (SVZ) in glioblastoma is associated with poor prognosis and is associated with specific tumor-biological characteristics. The SVZ microenvironment can influence gene expression in glioblastoma cells in preclinical models. We aimed to investigate whether the SVZ microenvironment has any influence on intratumoral gene expression patterns in glioblastoma patients. Methods: The publicly available Ivy Glioblastoma database contains clinical, radiological and whole exome sequencing data from multiple regions from resected glioblastomas. SVZ involvement of the various tissue samples was evaluated on MRI scans. In tumors that contacted the SVZ, we performed gene expression analyses and gene set enrichment analyses to compare gene (set) expression in tumor regions within the SVZ to tumor regions outside the SVZ. We also compared these samples to glioblastomas that did not contact the SVZ. Results: Within glioblastomas that contacted the SVZ, tissue samples within the SVZ showed enrichment of gene sets involved in (epithelial-)mesenchymal transition, NF-κB and STAT3 signaling, angiogenesis and hypoxia, compared to the samples outside of the SVZ region from the same tumors (p < 0.05, FDR < 0.25). Comparison of glioblastoma samples within the SVZ region to samples from tumors that did not contact the SVZ yielded similar results. In contrast, we observed no differences when comparing the samples outside of the SVZ from SVZ-contacting glioblastomas with samples from glioblastomas that did not contact the SVZ at all. Conclusion: Glioblastoma samples in the SVZ region are enriched for increased (epithelial-)mesenchymal transition and angiogenesis/hypoxia signaling, possibly mediated by the SVZ microenvironment.


Introduction
Glioblastoma is the most common primary malignant brain tumor. The prognosis is poor, with a median survival of about 15 months, despite intensive treatment including maximal safe resection and temozolomide-based chemoradiation [1].
The subventricular zone (SVZ), lining the walls of the lateral ventricles of the brain, has been of increasing interest in recent glioblastoma research. One key aspect of the SVZ is the neural stem cell niche it harbors. These astrocyte-like neural stem cells are hypothesized as the cells of origin of glioblastoma [2][3][4]. SVZ contact is associated with poor prognosis in glioblastoma patients. The underlying mechanisms of this adverse prognostic effect are possibly mediated by intrinsic tumor-biological characteristics. However, studies that have focused on exploring tumor-biological differences between glioblastomas with and without SVZ contact have shown inconclusive results [5][6][7][8][9][10][11][12]. In all previous gene expression studies on this subject, bulk tumor data were analyzed. Intratumoral heterogeneity, a known characteristic in glioblastoma, was not taken into account [13,14]. For example, it has been shown that multiple subtypes of the Verhaak classifier for glioblastoma (i.e., proneural, neural, classical and mesenchymal) can be found in the same tumor [14]. It is proposed that the tumor microenvironment contributes to this intratumoral heterogeneity [15].
Hence, our study focused on exploring oncogenic signaling in patient glioblastoma samples within the SVZ microenvironment, to explore a possible influence of the SVZ microenvironment on glioblastoma gene expression. We made use of glioblastoma tissue samples from the publicly available Ivy Glioblastoma Atlas, which contains clinical, radiological and whole exome sequencing data from multiple regions from en bloc resected glioblastomas [16]. In the tumors that contacted the SVZ, we performed gene expression analyses and gene set enrichment analyses to compare gene (set) expression in tumor regions within the SVZ to tumor regions outside the SVZ, within the same tumors. We also compared these samples to glioblastomas that made no contact with the SVZ.
We hypothesized that gene expression in glioblastoma samples from within the SVZ region show unfavorable oncogenic signaling characteristics, as a possible effect of the SVZ microenvironment.

Patient Cohort and Baseline Characteristics
Patients were selected from the Ivy Glioblastoma Atlas, an online accessible database that contains clinical, genomic and histologic information on patients with glioblastoma (http://glioblastoma.alleninstitute.org/ accessed on 15 March 2019) [16]. Adult patients with a primary IDH1-wildtype glioblastoma were included. Subsequently, clinical features including gender, age, initial KPS, extent of resection, adjuvant radiotherapy and/or chemotherapy and survival were collected. Additionally, molecular characteristics including EGFR amplification, MGMT methylation, PTEN status and molecular subtype (classical, mesenchymal, proneural and neural) [13] were collected. The baseline characteristics were compared between patients with SVZ-contacting glioblastoma and patients with glioblastoma without SVZ contact. Survival of patients with SVZ-contacting glioblastoma and patients with glioblastoma without SVZ involvement was compared with Kaplan-Meier curves and a log rank test.

MRI Analysis and Tissue Block Selection
SVZ involvement in glioblastomas was inspected on pre-operative T1-weighted postcontrast MRI images. SVZ contact was defined as contact of contrast-enhancing tumor with a lateral ventricle. The presence of SVZ contact was evaluated independently by two investigators (D.J.Z.D. and S.B.), with blinding to the clinical and genomic data. In case of disagreement, a third investigator (T.J.S.) gave the conclusive statement.
Along with the MRI scans, macroscopic images of en bloc resected tumors were available. The spatial orientation of the resected tumor was provided (i.e., anterior, posterior, lateral and medial side). The tumor tissue had been divided subsequently into multiple tissue blocks, for in situ hybridization and RNA expression analysis, as described in the original paper [16]. In this way, tissue blocks could be matched to their anatomical location on the MRI images, and tissue blocks from the SVZ region were selected.
These selection processes resulted in three groups of samples for further analysis and comparison: (a) glioblastoma samples within the SVZ (withinSVZ-samples), (b) samples outside of the SVZ from the same SVZ-contacting glioblastomas (outsideSVZ-samples), and (c) samples from tumors without any SVZ contact (noSVZcontact-samples).

RNA Sequencing Data Selection
Sampling of tissue blocks for RNA sequencing is described in the original paper [16] and was based on anatomical structural features that are commonly seen by neuropathologists in glioblastoma tissue sections stained with hematoxylin and eosin (H&E). The major structural regions were Leading Edge, Cellular Tumor and Infiltrating Tumor. Leading Edge was defined as the outermost boundary of the tumor, where the ratio of tumor to normal cells is about 1-3/100. Cellular Tumor constitutes the major part of core, where the ratio of tumor to normal cells is about 100/1 to 500/1. Infiltrating Tumor is the intermediate zone between the Leading Edge and Cellular Tumor, where the ratio of tumor cells to normal cells is about 10-20/100 [16]. In our study, only RNA sequencing data sampled by Cellular Tumor were included, as this yields the most tumor-specific data.

RNA Sequencing Data Analysis
Detailed information on tissue processing, RNA isolation, sequencing, quality control and alignment was described previously [16]. Subsequent RNA sequencing analysis was performed with aligned files in bam format with R (version 3.5.2). Count matrices were generated with the GenomicAlignments package (version 1.18.1). Genes with low read counts were dropped and TMM normalization for library size was performed to eliminate composition biases between libraries (edgeR package, version 3.24.3). We performed differential expression analyses by fitting a linear mixed model, with the location of the sample as fixed effect and patient ID as random effect, with the dream (differential expression analysis for repeated measures) analysis from the variancePartition package (version 1.12.3). We used mixed models, as this analysis method allows for correction of gene expression correlations in tumor samples from the same patient. Genes with a p-value adjusted for multiple testing with Benjamini-Hochberg's False Discovery Rate (FDR) below 0.05 were considered to be significantly differentially expressed. A heatmap was constructed to visualize the most differentially expressed genes across tumor locations (gplots package).
Subsequently, we performed Quantitative Set Analysis for Gene Expression (qusage package) in order to identify differential enrichment of gene sets between the groups. The ggen function of the qusage package (version 2.16.1) allowed us to incorporate the linear mixed model with the location of the tissue sample as fixed effect and patient ID as random effect. This analysis was performed with the hallmark gene sets [17] from the Molecular Signatures Database. Enriched gene sets with a p < 0.05 and a false discovery rate (FDR) <0.25 were considered significant. The results were visualized with pathway enrichment plots (qusage package).

Results
The Ivy Glioblastoma Atlas contains a total of 41 glioblastoma patients, of whom 34 had IDH1 wildtype primary glioblastoma ( Figure 1). Due to absence of spatial orientation data of the tumor or absence of RNA sequencing data, eight patients were excluded ( Figure 1). As a result, a total of 26 patients were included in our analysis ( Figure 2).    Patients with SVZ-contacting glioblastoma did not differ significantly from patients with glioblastoma without SVZ contact with respect to age, initial KPS, extent of resection, treatment, PTEN status, MGMT methylation status, EGFR amplification status and EG-FRvIII mutation status (Table 1). Median overall survival was shorter in patients with SVZ-contacting glioblastoma (442 days vs. 544 days in patients with glioblastoma without SVZ contact), but this observation did not reach a statistically significant difference. (Table  1, Supplementary Figure S1). Finally, molecular tumor subtype(s) were comparable in SVZ-contacting glioblastoma and glioblastoma without SVZ contact. Often, more than one subtype was found in the same tumor (Supplementary Table S1).  Patients with SVZ-contacting glioblastoma did not differ significantly from patients with glioblastoma without SVZ contact with respect to age, initial KPS, extent of resection, treatment, PTEN status, MGMT methylation status, EGFR amplification status and EGFRvIII mutation status (Table 1). Median overall survival was shorter in patients with SVZ-contacting glioblastoma (442 days vs. 544 days in patients with glioblastoma without SVZ contact), but this observation did not reach a statistically significant difference. (Table 1, Supplementary Figure S1). Finally, molecular tumor subtype(s) were comparable in SVZ-contacting glioblastoma and glioblastoma without SVZ contact. Often, more than one subtype was found in the same tumor (Supplementary Table S1).

Differential Enrichment of Oncogenic Gene Sets across Tumor Regions: Intratumoral Analysis
First, we analyzed RNA expression in all tumors with SVZ contact. Within the same tumors, we compared the withinSVZ-samples to the outsideSVZ-samples. In this analysis, 45 tumor samples (nine withinSVZ-samples and 36 outsideSVZ-samples) were included, from sixteen patients (Figure 1).
After adjustment for multiple testing, no single gene was significantly differentially expressed between tumor samples from within and outside of the SVZ (Supplementary Figure S2). Gene set expression analysis, however, showed differential expression of 24 gene sets (p < 0.05, FDR < 0.25) (Supplementary Table S2, Figure 3). In the with-inSVZ-samples, we observed the most increased enrichment of gene sets involved in (epithelial-)mesenchymal transition, angiogenesis, and NF-κB and JAK/STAT3 signaling compared to the outsideSVZ-samples from the same tumors (Supplementary Table S2).

Differential Enrichment of Oncogenic Gene Sets in Glioblastomas in the SVZ Region: Intertumoral Analysis
In this analysis, 37 tumor samples (nine withinSVZ-samples and 28 noSVZcontactsamples, Figure 2) from 17 patients were included.
The only significantly differentially expressed gene between the withinSVZ-samples and the noSVZcontact-samples was DCAF4, which was downregulated in the withinSVZgroup (logFC = −1.50, adjusted p = 0.002). The separate expression values (log counts per million) of DCAF4 per sample are shown in Supplementary Figure S3. The results for the fifty genes with the lowest (unadjusted) p-values are shown in a heatmap (Supplementary Figure S4).
Gene set enrichment analysis showed differential expression of 27 gene sets (p < 0.05, FDR < 0.25, Supplementary Table S3, Figure 4). Again, we found relatively increased enrichment of gene sets involved in (epithelial-)mesenchymal transition, angiogenesis, and NF-κB and JAK/STAT3 signaling in the withinSVZ-samples.

No Difference in Oncogenic Signaling between Glioblastoma Samples from Outside the SVZ and Tumors without SVZ Contact: Intertumoral Analysis
In this analysis, 64 tumor samples (36 outsideSVZ-samples and 28 noSVZcontactsamples) were included, from 24 patients (Figure 1).
No significantly differentially expressed genes between the outsideSVZ-samples and noSVZcontact-samples were found (adjusted p-values > 0.50 for all genes). The results for the 50 genes with the lowest adjusted p-values are shown in a heatmap (Supplementary Figure S5). No differential enrichment of gene sets (with an FDR less than 0.25) were found.

Discussion
With intra-and intertumoral gene expression analyses, we show that gene sets associated with mesenchymal transition are relatively enriched in glioblastoma tissue in the SVZ region, as compared to tumor tissue outside the SVZ region. To our knowledge, this is the first study in which gene expression in glioblastomas contacting the SVZ is regionally explored to this end so far. Our findings suggest that the SVZ microenvironment could influence oncogenic signaling in glioblastoma.
The differential enrichment of gene sets was more pronounced in the comparison of withinSVZ-samples versus noSVZcontact-samples, than in the comparison of different samples from the SVZ-contacting tumors (withinSVZ versus outsideSVZ), which could reflect a greater degree of similarity in the samples in the latter comparison. The comparison of outsideSVZ-samples from the SVZ-contacting tumors with noSVZcontact-samples from glioblastomas without SVZ contact showed no significantly different gene or pathway expression. This is in line with the hypothesis that the SVZ microenvironment could influence region-specific gene expression in glioblastoma.
Few studies have focused on the interaction of the SVZ microenvironment and glioblastoma (stem) cells (GSCs) [18][19][20]. This is of particular interest, as the SVZ niche is believed to serve as a GSC reservoir which contributes to resistance to therapy. It is proposed that the microenvironment in the SVZ closely interacts with GSCs in order to establish this protective niche [18,19,21]. One study showed that GSCs in the SVZ appear to have an enhanced mesenchymal signature compared with their counterparts from the tumor [19]. These mesenchymal features, including a higher expression level of N-cadherin and vimentin, were shown to be upregulated by SVZ-released CXCL12. Moreover, inhibition of the CXCL12/CXCR4 signaling axis with AMD3100 (a CXCL12/CXCR4 antagonist) weakened the tumor's mesenchymal signature in the SVZ and increased the tumor's sensitivity to radiotherapy [19]. This correlation of mesenchymal activation in glioblastoma and resistance to radiotherapy (and chemotherapy) has been reported in other studies as well [22]. Another study found a similar predominance of the mesenchymal subtype in glioblastoma samples from the SVZ as compared to SVZ distant samples from the same tumor [20]. Moreover, in this study, isolated GSCs from the SVZ and GSCs from other tumor mass of the same glioblastoma showed different patterns of response to therapies [20]. In addition, neural-precursor cells can secrete factors such as PTN and ROCK that were shown to drive glioblastoma invasion to the SVZ. In our analysis, for instance, ROCK1 is included in, among others, the HALLMARK_APOPTOSIS pathway, which is enriched in withinSVZ tissue samples [23].
In this light, our results provide further evidence to this SVZ region-specific mesenchymal transition in glioblastoma. However, a deeper insight into the complex interaction between glioblastoma cells (and GSCs in particular) and the SVZ micro-environment is needed to unravel the mechanism behind this mesenchymal shift.
In addition to the (epithelial-)mesenchymal transition signaling pathway, several other gene sets linked to aggressive growth of glioblastoma were relatively upregulated in the withinSVZ-samples group, including TNF-α-mediated NF-κB activation, IL-6 induced STAT3 activation, TGF-β signaling, p53 signaling, KRas signaling, genes upregulated by reactive oxygen species, angiogenesis, hypoxia, coagulation, complement activation and inflammatory response. Activation of and interaction between several of these pathways have been linked to induction of (epithelial-)mesenchymal transition in glioblastoma as well [22,[24][25][26][27][28][29][30]. For example, the mesenchymal subtype of glioblastoma is characterized by increased levels of NF-κB signaling components and moreover, glioma sphere cultures of the proneural subtype can transform to a mesenchymal state under TNF-α-mediated NF-κB activation [22]. Additionally, it has been demonstrated that STAT3 activation increases when initial proneural tumors experience a mesenchymal shift upon recurrence [31,32]. Furthermore, STAT3 was shown to be highly expressed in glioblastoma stem cells (GSCs), which are tumor cells with self-renewing properties that contribute to tumor initiation and therapeutic resistance [33]. Finally, hypoxia and reactive oxygen species could also potentially induce mesenchymal transition in glioblastoma [27,34].
The only significantly differentially expressed gene (DCAF4) was observed in the with-inSVZ-samples vs. noSVZcontact-samples analysis. DCAF4 (DDB1 and CUL4-associated factor 4) encodes a WD (Trp-Asp) repeat-containing protein that interacts with the CUL4-DDB1 E3 ligase macromolecular complex. The CUL4-DDB1 ubiquitin ligase is involved in cell proliferation, survival, DNA repair and genomic integrity [35]. CUL4A demonstrated the potential to promote (epithelial-)mesenchymal transition through the activation of ZEB1 in breast cancer cells [36]. One study showed that overexpression of DCAF4L2, a paralog of DCAF4, could support (epithelial-)mesenchymal transition in colorectal cancer cells through activation of the NF-κB pathway [37], and genetic variations in DCAF4 have been previously associated with leukocyte telomere length, keloid formation and lung cancer risk [38][39][40]. However, as we found a relative downregulation of DCAF4 in nearly all withinSVZ-samples, the relevance of this observation compared to increased expression of gene sets involved in (epithelial-) mesenchymal transition in tumor samples from the SVZ region remains unclear.
Our study has some technical limitations that should be considered. We relied on the availability of data from the comprehensive Ivy Glioblastoma Atlas. RNA sequencing data were not available for all tissue regions and spatial orientation of the tumor tissues was not clear in some cases. As a consequence, only a limited number of patients could be included in our analyses. Especially for the samples within the SVZ region, the number of available tumor samples was limited (nine samples from seven patients). Given the small number of patients, we were unable to study our findings in relevant subgroups, such as different glioblastoma subtypes or varying subventricular subregions.
Although the gene set expression results of our separate subanalyses are in line with each other, gene set analyses in relatively small-sized cohorts should be regarded as exploratory, and must be validated in other studies prior to drawing any conclusions. Future studies could elucidate the region-specific molecular profiles in SVZ-regions with macroscopic tumor involvement as well as normal-appearing SVZ, through prospective collection of relevant tissue and detailed analysis with single-cell RNA sequencing.

Conclusions
Glioblastoma tissues within the SVZ region show increased expression of gene sets involved in (epithelial-)mesenchymal transition, angiogenesis, and NF-κB and JAK/STAT3 signaling, compared to tumor tissues from non-SVZ-regions, either from the same SVZcontacting tumors or from glioblastomas that do not contact the SVZ. Our results suggest that the SVZ microenvironment could influence regional gene expression in glioblastoma, to induce (epithelial-)mesenchymal transition and possibly a more invasive glioblastoma phenotype.

Data Availability Statement:
Publicly available datasets were analyzed in this study. This data can be found here: https://glioblastoma.alleninstitute.org/.

Conflicts of Interest:
The authors have no competing interests to declare.