Identification of a Prognostic Microenvironment-Related Gene Signature in Glioblastoma Patients Treated with Carmustine Wafers

Simple Summary Carmustine wafer (CW) implantation into the resection cavity of patients operated for glioblastoma (GBM) was approved as an adjuvant treatment before the Stupp Protocol. Although contrasting clinical results limited its use, our retrospective study on 116 GBM treated with CW showed a significant benefit in terms of OS in a subgroup of patients. Since GBM growth, progression, and drug resistance are supported by the surrounding environment, and since the tumor microenvironment (TME) is the source of druggable targets, we hypothesized that the TME of patients who benefited from CW could have different characteristics compared to patients who did not show any advantage. Exploiting a human in vitro model of glioma microenvironment and a transcriptomic approach, we found a different gene signature suggesting the importance of developing in vitro models that mimic the properties of human cancers and that can help to study individual patient characteristics at the cellular and molecular level. Abstract Despite the state-of-the-art treatment, patients diagnosed with glioblastoma (GBM) have a median overall survival (OS) of 14 months. The insertion of carmustine wafers (CWs) into the resection cavity as adjuvant treatment represents a promising option, although its use has been limited due to contrasting clinical results. Our retrospective evaluation of CW efficacy showed a significant improvement in terms of OS in a subgroup of patients. Given the crucial role of the tumor microenvironment (TME) in GBM progression and response to therapy, we hypothesized that the TME of patients who benefited from CW could have different properties compared to that of patients who did not show any advantage. Using an in vitro model of the glioma microenvironment, represented by glioma-associated-stem cells (GASC), we performed a transcriptomic analysis of GASC isolated from tumors of patients responsive and not responsive to CW to identify differentially expressed genes. We found different transcriptomic profiles, and we identified four genes, specifically down-regulated in GASC isolated from long-term survivors, correlated with clinical data deposited in the TCGA–GBM dataset. Our results highlight that studying the in vitro properties of patient-specific glioma microenvironments can help to identify molecular determinants potentially prognostic for patients treated with CW.


Introduction
The current standard therapy for glioblastoma (GBM), represented by maximal surgical resection combined with chemo-and radiotherapy, offers only a palliative treatment including their motility, in vitro [27]. We have also already demonstrated that GASCs interact with tumor cells (both Glioma Stem Cells (GSC) and immortalized tumor cell lines) through the release of exosomes [27,28]. For these reasons, GASCs represent bona fide Glioma Stromal Cells residing in the tumor microenvironment. Furthermore, we have reported that the transcriptomic profiling of GASC allows for the identification of prognostic gene signatures in low-grade and high-grade gliomas, suggesting their clinical relevance [29][30][31].
In the present work, we employed a transcriptomic approach to describe the gene expression profile of GASC isolated from tumors explanted before CW implantation. We chose to analyze and compare GASC derived from long-term (GASC-LS, OS ≥ 30 months and mean PFS of 24 ± 9 months) and short-term (GASC-SS, OS ≤ 14 months and mean PFS of 8 ± 2 months) survivors. We found 374 differentially expressed transcripts, and by focusing the analysis on a panel represented by 78 protein-coding genes, we observed that GASC-LS were characterized by a statistically significant down-regulation of four clinically relevant genes. These findings suggest that the TME could provide important biomarkers allowing to identify the group of patients for which a CW-based therapy is recommended over traditional approaches.

Surgical Procedure
Tumor tissues were collected from patients (age ≥ 18 years) affected by newly diagnosed GBMs arising de novo. All of the patients underwent surgery at the Neurosurgical Department of the Academic Hospital Udine (ASUFC). The criteria for CW implantation have been previously described [19]. Briefly, all of the patients did not undergo any previous surgery nor any adjuvant therapy and exhibited a closed surgical cavity. Multifocal lesions and/or lesions extending across the corpus callosum were excluded. The patients had at least 12 months of follow-up and were followed with neurological assessment and post-operative MRI every 4 months.
CWs were implanted in all patients after surgical tumor removal and the intraoperative pathological confirmation of GBM. They were not implanted when tumor removal required the creation of a large opening of a ventricle and/or the basal cistern. The research was approved by the local ethics committee (Parere 196/2014/Em). Written informed consent was obtained from all of the patients. Clinical investigations were conducted according to the principles expressed in the Declaration of Helsinki.

Isolation and Culture of Glioma Associated Stem Cells (GASC)
Glioma-associated stem cells (GASC) were isolated from surgical samples of a removed tumor before CW implantation and maintained in vitro, applying, with minor modifications, a protocol optimized for culturing multipotent adult stem cells from normal and neoplastic human tissues [26]. Briefly, the GBM fragments were first mechanically disaggregated with scalpels and then enzymatically dissociated in a 0.025% Collagenase type II solution (Worthington, Columbus, OH, USA) in Joklik-modified Eagle's Medium (Sigma-Aldrich, Saint Louis, MO, USA), for 5 min at 37 • C. Collagenase activity was stopped by adding 10% Fetal Bovine Serum in Joklik-modified Eagle's Medium (Sigma-Aldrich, Saint Louis, MO, USA). The cell suspension was centrifuged at 500× g for 10 min after less than 40 µm in diameter. Then, 2.0 × 10 6 freshly-isolated human cells were plated onto 100-mm diameter human fibronectin (Sigma-Aldrich, Saint Louis, MO, USA)-coated dishes (BD Falcon, Franklin Lakes, NJ, USA) in an expansion medium that was composed as follows: 60% low glucose DMEM (Invitrogen, Waltham, MA, USA), 40% MCDB-201, 1 mg/mL linoleic acid-BSA, 10 −9 M dexamethasone, 10 −4 M ascorbic acid-2 phosphate, 1× insulin-transferrin-sodium selenite (all from Sigma-Aldrich), 2% fetal bovine serum (StemCell Technologies, Cambridge, UK), 10 ng/mL of human PDGF-BB, 10 ng/mL of human EGF (both from Peprotech EC, London, UK). The clinical characteristics of all of the

GASC RNA Extraction, Library Preparation and Sequencing
We analyzed the GASC cell lines derived from patients who showed an OS ≥ 30 months (GASC-LS; n = 3) and an OS ≤ 14 months (GASC-SS; n = 2) evaluated after CW implantation at the time of surgery. The total RNA was extracted using the RNAeasy mini kit (Qiagen GmbH, Hilden, Germany) according to the manufacturer's instructions. The Universal Plus mRNA-Seq kit (Tecan Genomics, Redwood City, CA, USA) was used for library preparation following the manufacturer's instructions (library type: fr-secondstrand). The RNA samples were quantified and quality tested by Agilent 2100 Bioanalyzer RNA assay (Agilent Technologies, Santa Clara, CA, USA) or by Caliper LabChip GX (PerkinElmer, Waltham, MA, USA). The final libraries were checked with both Qubit 2.0 Fluorometer (Invitrogen, Carlsbad, CA, USA) and Agilent Bioanalyzer DNA assay or by the Caliper LabChip GX (PerkinElmer, Waltham, MA, USA). The libraries were then prepared for sequencing and sequenced using the single-end 150 bp mode on a NovaSeq 6000 (Illumina, San Diego, CA, USA). Quality control for the raw sequencing reads was performed using the FastQC (v0.11.9) [32] and MultiQC (v1.09) software [33]. Quality, adapters, and contamination filtering were performed using the Trimmomatic command-line tool [34]. The processed reads were aligned to the NCBI GRCh38 human reference using STAR (v2.7.1a) [35]. Transcripts assembly and the number of reads per gene were determined using Stringtie (v1.3.6) [36]. Differentially expressed genes were identified using the DESeq2 (v1.26.0) R/Bioconductor package [37], applying the Wald test; we adjusted for multiple hypothesis testing by employing the Benjamini-Hochberg correction, considering as statistically significant the results having abs(log2FC) ≥ 1, FDR < 0.05. Extended gene annotations (including HGNC gene symbol, description, and transcript type) were obtained using the biomaRt (v2.42.0) R/Bioconductor package [38].

Functional Enrichment Analysis
Significantly up-and down-regulated genes identified in the differential gene expression analysis were used to query the MSigDB database, investigating statistically relevant biological associations. Specifically, we examined the C4, C6, CP, and GeneOntology: Bio-logicalProcess, H, and IMMUNESIGDB ontologies: C4 gathers cancer-related signatures originated by the data mining of the large microarray data, C6 collects signatures related to pathways deregulated in cancer, CP includes data from canonical pathways (i.e., KEGG, BIOCARTA, Reactome, PID, and WikiPathways), and IMMUNESIGDB contains gene signatures derived from chemical and genetic perturbations of the immune system. The top-20 significantly-enriched (FDR q-value < 0.05) genesets were retrieved.

RT2 Profiler™ PCR Array
Seventy-Eight genes were selected to create a customized RT2 Profiler PCR Array (CLAH43115; Qiagen Gmbh, Germany). cDNA was synthesized using the RT2 First Strand Kit (Qiagen) following the manufacturer's instructions. According to the manufacturer's protocol, real-time PCR was performed using the RT2 Profiler PCR Arrays in combination with the RT2 SYBR Green/ROX PCR Master Mix (Qiagen). A 102-µL cDNA synthesis reaction volume was mixed with 2 × RT2 SYBR Green Mastermix and RNase-free water to obtain a total volume of 2700 µL. Subsequently, 25 µL of the PCR component mix was dispensed into each of the 96-well PCR arrays. The cycling program, performed on a Roche LightCycler 480, comprises three steps: 95 • C for 10 min for 1 cycle, 45 cycles at 95 • C for 15 s, and 60 • C for 60 s. Data analyses were performed using the manufacturer's software (https: //geneglobe.qiagen.com/us/analyze/ (accessed on 20 January 2022)), which calculates the fold change/regulation of the investigated genes using the delta Ct (∆Ct) method. Briefly, ∆Ct was calculated between each gene and an average of the housekeeping genes (ACTB, B2M, GAPDH, HPRT1, and RPLP0). Then, ∆∆Ct was extrapolated as the difference between ∆Ct of genes in the tested group (GASC-LS) and ∆Ct of the same genes in the control group (GASC-SS). Finally, fold change was calculated using 2 (−∆∆Ct) formula.

Differential Gene Expression Analysis in the TCGA-GBM Dataset
The significantly differentially expressed genes identified with the RT2 Profiler array were tested by comparing their expression levels in the TCGA-GBM dataset (n = 163) with respect to the TCGA-GTEx matched normal samples (n = 207). The data were obtained from the GEPIA web server (http://gepia.cancer-pku.cn/index.html, last accessed on 1 February 2022) and summarized as boxplots; |Log2FC| Cutoff: 1; p-value Cutoff: 0.01.

Survival Analysis
The genes differentially expressed according to the RT2 profiler array underwent a survival analysis to define their prognostic value in the TCGA-GBM dataset on a single gene basis in terms of overall survival (OS), using the GEPIA2 web server (http://http: //gepia2.cancer-pku.cn/, last accessed on 1 February 2022). Afterward, we defined a four-genes signature and tested it in the same way. Median or tertile cutpoints were used to stratify patients.

Transcriptomic Characterization of GASC through RNA-Seq
The gene expression profile of GASC isolated from GBM of patients who underwent CW implantation during surgery before the standard Stupp protocol was evaluated by transcriptomic analysis. We compared three GASC cell lines derived from patients with OS ≥ 30 months, ranging from 30 to 43 months, and median PFS of 33 months, ranging from 30 to 41 months, (GASC-LS) and two GASC cell lines derived from patients with OS ≤ 14 months, ranging from 12 to 16 months, and median PFS of 6 months, ranging from 4 to 8 months (GASC-SS). Supplementary Table S2 displays the 374 differentially expressed transcripts (protein-coding, pseudogenes, and lncRNA), 221 up-and 153 downregulated in GASC-LS compared to GASC-SS. The significantly up-and down-regulated genes identified in the differential gene expression analysis were used separately to query the MSigDB database, investigating statistically relevant biological associations. Specifically, we examined five major collections of gene sets present in the database: GeneOntology: BiologicalProcess, cancer-related signatures, pathways deregulated in cancer, canonical pathways (i.e., KEGG, BIOCARTA, Reactome, PID, and WikiPathways), and gene signatures deriving from chemical and genetic perturbations of the immune system. The top 20 significantly enriched (FDR q-value < 0.05) genesets for up (Table 1) and down-regulated genes ( Table 2) in GASC-LS were reported. The most represented enriched terms found between genes up-regulated in GASC-LS were associated with neurogenesis, cell and neuron differentiation, central nervous system development, and secretion, while the down-regulated genes were associated with cell-cell signaling, defense response, and the regulation of the immune system. Interestingly, both up-and down-regulated genes were associated with two common biological processes, namely regulation of ion transport and regulation of transport.

A Customized RT-PCR Array Identified Four Genes Specifically Modulated in GASC-LS
To give insight into the differences between GASC-LS and GASC-SS transcriptomic profiles, we manually selected, among the 374 differentially expressed genes obtained by RNA-seq, transcripts with abs(logFC) ≥ 2 described as "protein-coding", and we identified 78 genes. To evaluate their expression by real-time PCR, we created a customized RT2 profiler PCR array panel, and we validated their expression on the same cell lines used for RNAseq and extended the analysis to three other GASC-LS (OS ≥ 30 months and mean PFS of 24 ± 9 months) and three other GASC-SS (OS ≤ 14 months and mean PFS of 8 ± 2 months). Figure 1 shows the volcano plot representing the changes in gene expression, derived from the comparison of six GASC-LS with five GASC-SS by plotting the log2 of the fold changes Cancers 2022, 14, 3413 6 of 14 (greater than 2, red dots, or less than 2, green dots) on the x-axis versus their statistical significance (p ≤ 0.05). For most of the genes in the array, we did not find any statistically significant difference. However, four genes resulted in significantly down-regulated (green dots in the upper left side of the graph) in GASC-LS: ALPL (Alkaline Phosphatase), GPR68 (G protein-coupled receptor 68), NETO1 (neuropilin and tolloid-like 1), and VGF (or VGF nerve growth factor inducible), (Table 3). Table 1. The Top 20 significantly enriched MSigDB pathways associated with the specific lists of genes up-regulated in GASC-LS. The table reports the overall number (#) of genes involved in the annotated pathways, the number (#) of genes provided as input, and the associated p-value.  their statistical significance (p ≤ 0.05). For most of the genes in the array, we did not find any statistically significant difference. However, four genes resulted in significantly down-regulated (green dots in the upper left side of the graph) in GASC-LS: ALPL (Alkaline Phosphatase), GPR68 (G protein-coupled receptor 68), NETO1 (neuropilin and tolloid-like 1), and VGF (or VGF nerve growth factor inducible), (Table 3).

Clinical Correlation of Genes Down-Regulated in GASC-LS
To investigate if the four genes significantly down-regulated in GASC-LS compared to GASC-SS, according to the RT2 profiler array, could have clinical relevance, we evaluated their expression in 163 GBM patients included in the Cancer Genome Atlas GBM dataset (TCGA-GBM). To this aim, we queried the GEPIA web server comparing the TCGA RNA-Seq data of GBM patients to that of 207 TCGA/GTEx matched normal samples. As shown in Figure 2, ALPL, GPR68, NETO1, and VGF were all up-regulated in normal brain tissue, suggesting their specificity for non-tumoral cells-microenvironmental, as GASCs were described. Afterward, we performed a survival analysis with these genes to establish their prognostic value in the TCGA-GBM dataset on a single gene basis in terms of overall survival (OS, defined as the time between cancer diagnosis and death for any cause), using the GEPIA2 web server. We observed that patients with low expression of ALPL, GPR68, and VGF show a significantly longer OS. A similar, although not statistically significant, the trend was reported for NETO1. Moreover, these four genes, when taken together as a whole signature, retained their prognostic value (Figure 3). When we repeated the same type of analysis in terms of PFS, we found that the low expression of these genes is correlated with a longer PFS, although the correlation is statistically significant only for GPR68 (Supplementary Figure S1). Altogether, these results suggest the identification of four genes that should be specifically expressed in cells of the GBM microenvironment. Interestingly, they were found to be down-regulated in patients that, after CW treatment, have shown an OS ≥30 months, thus suggesting the presence of a less tumor-supporting environment.

Discussion
Outside the standard of care for GBM treatment, carmustine wafer insertion into the resection cavity has been approved as adjuvant therapy, giving a bridge between surgery and the initiation of the Stupp protocol [7,8]. The benefit of CW implantation is debated in the literature due to reported contrasting results, side effects, its costs, and the exclusion of treated patients from other trials [39,40]. However, a retrospective study performed in the neurosurgery unit of the Academic Hospital of Udine (ASUFC) found a category of patients that really benefited from CW, in terms of OS (≥36 months), with negligible side effects [19]. In the present work, we took advantage of case studies treated with CW and of an in vitro model of the GBM microenvironment, represented by the GASC cell lines, to give more insights into which determinants could influence the response to CW at the cellular level. We initially performed an explorative transcriptomic analysis of GASC-LS and GASC-SS to find out which differences characterized the TME of patients who showed a long-term OS (≥30 months) with respect to those who did not show any improvement of OS (≤14 months), after receiving CW, at the time of surgery. We found that 374 genes were differentially expressed in GASC-LS compared to GASC-SS in a statistically significant manner. The functional enrichment analysis revealed that the most represented enriched terms found in up-regulated genes in GASC-LS were associated with neurogenesis, cell and neuron differentiation, and central nervous system development. These findings possibly indicate that the microenvironment of GBM responsive to CW has an asset more differentiated and more committed to administrating central nervous system basic functions rather than supporting tumor growth. Moreover, the down-regulated genes were associated with cell-cell signaling and the regulation of the immune system, describing a poor communicative TME, thus suggesting a possible interference in the cross-talk with the tumor, which is thought to support its growth and recurrence. Interestingly, both up-and down-regulated genes in GASC-LS were associated with two common pathways: the regulation of ion transport and regulation of transport, represented mainly by genes coding for ion channels.
The transport of ions across the cell membrane is a fundamental process for maintaining normal cellular functions and activity (cell cycle, cell death, cell volume regulation, and proliferation) [41,42]. Moreover, it is well accepted that cancers of the nervous system cross-talk within the local tumor microenvironment. Communication between cancer cells and neurones is driven by synapses through the release of neurotransmitters and voltagegated mechanisms to regulate cancer cell growth [43,44]. Furthermore, glioma cells can electrically integrate into neural circuits through neuro-glioma synapses [43]. Ion channels are also involved in pathways modulating tumor vascularisation and tumor-immune cell interactions [45]. It is likely that the ion channels' activity and their dysregulations underlie several known hallmarks of glioma, such as proliferative capacity, apoptotic escape, and invasion, and they are now regarded as new therapeutic targets [46][47][48][49]. The finding that most of the genes differentially regulated in GASC-LS and GASC-SS belong to the "ion channels" category suggests that they play an important role in our model, worthy of further investigation. To analyze the expression profile of GASC-LS and GASC-SS in more detail, we selected 78 genes with an abs(logFc) ≥ 2, described as "protein-coding", and we performed a gene expression assay using a customized RT2-Profiler PCR array. Most of these genes are not differentially expressed between GASC-LS and GASC-SS, except for ALPL, GPR68, NETO1, and VGF, which were significantly down-regulated in GASC-LS. The role of these four genes in the GBM microenvironment is unknown, although their expression has been described in other types of tumors. ALPL is reported as a marker of embryonic stem cells, and it is highly expressed in induced-Pluripotent Stem Cells (iPSCs) obtained from human fibroblasts [50]. Moreover, a role of ALPL was described in malignant leukemia, in which changes in its levels can be used to identify rare populations of highly refractory malignant cells [51].
GPR68 is a GPCRs (G protein-coupled receptor)'s family member. GPCRs modulate signal transduction pathways and cellular processes that are critical for physiological functions [52] and for the initiation and progression of tumors [53,54]. GPR68 is a protonsensing G-protein-coupled receptor regulated by extracellular acidity and involved in the regulation of a variety of cellular functions. Acidosis is considered a hallmark of the tumor microenvironment. Multiple factors, such as hypoxia, inflammation, and glycolytic cell metabolism, contribute to creating an acidic TME [55][56][57], and it has been shown that acidosis modulates cell proliferation, apoptosis tumor progression, metastasis, anti-tumor immunity, and angiogenesis [58,59]. The expression of GPR68 is highly up-regulated in numerous types of cancer, including prostate, colon and pancreatic tumors, melanoma, myelodysplastic syndrome, and medulloblastoma [60]. Emerging evidence has revealed that GPR68 may play crucial roles in the biology of pancreatic ductal adenocarcinoma (PDAC). Its activation induces Cancer-Associated Fibroblast (CAFs) to acquire an enhanced tumor-promoting phenotype and promotes PDAC cell proliferation [61]. NETO1 is a gene coding for a recently described protein involved in the modulation of glutamate receptor activity [62,63]. In particular, Neto1 is an auxiliary subunit modulating the gating properties of Kainate receptors (KARs), a subfamily of ionotropic glutamate receptors that mediate excitatory transmission and regulate neurotransmitter release at the preand post-synaptic level [64,65]. Although Neto1 was not previously described in glioma, KARs have been implicated in epilepsy and other neuropsychiatric conditions [66,67]; moreover, the high expression of Neto1 was associated with metastatic ovarian carcinomas, in which it increases the migration and invasion capability of cancer cells, modulating actin cytoskeletal dynamics [68]. For these reasons, investigating their implication and modulation on highly invasive tumors such as GBM would be of great interest. Finally, VGF encodes a neuropeptide precursor expressed in several types of neurons in the central and peripheral nervous system [69][70][71][72]. The role of VGF in tumors is poorly described. In lung cancer, it seems to be associated with the acquisition of resistance to EGFR inhibitors [73], and in breast cancer, it is involved in the epithelial-to-mesenchymal transition [74]. A recent study identified VGF as a key player in the bidirectional influence between Glioma stem cells (GSC) and their differentiated progeny (DGCs, differentiated glioma cells) [75]. In the model described by Wang et al., VGF contributes to glioma progression by maintaining selfrenewal and proliferation of both GSC and DGCs. Further investigation will clarify if the different expression of VGF, in our model, correlates with a different level of communication intervening into the GBM's TME of patients with long-term and short-term survival, causing an opposed tumor supporting action. Since the present study is only descriptive and performed on a restricted number of cases, we attempted to improve our results by evaluating the expression of ALPL, GPR68, NETO1, and VGF in the Cancer Genome Atlas GBM dataset (TCGA-GBM). We found that these four genes are predominantly expressed in healthy brains rather than in GBM. Moreover, a survival analysis highlighted that the low expression of each single gene was associated with patients showing a significantly longer OS. This result was also confirmed when ALPL, GPR68, NETO1, and VGF were evaluated together. Moreover, the low expression of these four genes showed a correlation with a longer PFS, although in a not statistically significant manner, thus confirming their possible clinical relevance. In summary, here, we described the gene expression profile of an in vitro model of GBM's TME, represented by the GASC cell lines. We found that GASC-LS displayed a profile not directly correlated with response to therapy but suggestive of a less tumor-promoting action. Of course, future studies will be required to better explain these results and to clarify the mechanistic functions of ALPL, GPR68, NETO1, and VGF in the GBM TME and how they can mediate the success of CW application in the tumor resection cavity.

Conclusions
The efficacy of carmustine wafers in the treatment of GBM is a still debated question, and the identification of patients who may really benefit from them is needed. Here, we showed that studying an in vitro model of patient-derived TME using a transcriptomic approach highlighted differences in terms of gene expression between patients selected for the response to CW. Deepening the knowledge of these differences will help to explain how the TME molecular landscape is connected with the success of CW application.

Informed Consent Statement:
Written informed consent has been obtained from the patient(s) to publish this paper.

Data Availability Statement:
The RNAseq data are available as GEO accession GSE199407 Go to https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE199407. Enter token evqzukemxzgfhap into the box.