Deciphering of Adult Glioma Vulnerabilities through Expression Pattern Analysis of GABA, Glutamate and Calcium Neurotransmitter Genes

Adult infiltrating gliomas are highly aggressive tumors of the central nervous system with a dismal prognosis despite intensive multimodal therapy (chemotherapy and/or radiotherapy). In this study, we studied the expression, methylation and interacting miRNA profiles of GABA-, glutamate- and calcium-related genes in 661 adult infiltrating gliomas available through the TCGA database. Neurotransmitter-based unsupervised clustering identified three established glioma molecular subgroups that parallel major World Health Organization glioma subclasses (IDH-wildtype astrocytomas, IDH-mutant astrocytomas, IDH-mutant oligodendroglioma). In addition, this analysis also defined a novel, neurotransmitter-related glioma subgroup (NT-1), mostly comprised of IDH-mutated gliomas and characterized by the overexpression of neurotransmitter-related genes. Lower expression of neurotransmission-related genes was correlated with increased aggressivity in hypomethylated IDH-wildtype tumors. There were also significant differences in the composition of the tumor inflammatory microenvironment between neurotransmission-based tumor categories, with lower estimated pools of M2-phenotype macrophages in NT-1 gliomas. This multi-omics analysis of the neurotransmission expression landscape of TCGA gliomas—which highlights the existence of neurotransmission-based glioma categories with different expression, epigenetic and inflammatory profiles—supports the existence of operational neurotransmitter signaling pathways in adult gliomas. These findings could shed new light on potential vulnerabilities to exploit in future glioma-targeting drug therapies.


Introduction
Infiltrating gliomas are the most common malignant tumors of the central nervous system in adults. They represent nearly 45-50% of malignant primary brain neoplasms [1] and are associated with relatively short survivals [2]. The infiltrative nature, intratumoral heterogeneity [3] and cellular signaling complexities of these aggressive tumors make them a major challenge to overcome in terms of therapy and personalized patient management [4,5]. Pathological and adaptative interactions with the surrounding microenvironment (non-neoplastic glia, neurons and immune cells) further contribute to tumor aggressivity and progression [6]. Current World Health Organization diagnostic classification schemes for adult infiltrating gliomas are based on the presence of molecular alterations such as isocitrate dehydrogenase 1/2 mutations (IDH1/2) and 1p/19q codeletion [1]. While IDH-wildtype gliomas usually present with high-grade histology and correlate to short survival [7,8], IDH-mutated, 1p/19q-non-codeleted astrocytomas and IDH-mutated 1p/19q-codeleted oligodendrogliomas are usually first diagnosed as low-grade gliomas and associated with longer survival [7,9]. IDH1/2 genes code for isocitrate dehydrogenase cytoplasmic enzymes that play a crucial role in cellular energy metabolism 7 . One of their major functions is the catalysis of citrate oxidative decarboxylation to alpha-ketoglutarate (α-KG), which is a well-known intermediate of the TCA cycle [10]. Further, α-KG is involved in neurotransmission by serving as a precursor of two important human brain neurotransmitters: glutamate and gamma-aminobutyric acid (GABA) [11]. In glioma, mutations in IDH1/2 lead neoconversion of α-KG into the oncometabolite D-2-hydroxyglutarate (D-2-HG) [12]. Consequently, the IDH status of gliomas impacts not only their energetic metabolism but also their integration into surrounding neural circuits due to the potential dysregulation in neurotransmitter metabolism, thereby affecting tumor progression [13].
Aside from their leading role in neuronal synaptic transmission, GABA and glutamate can regulate many other brain biological processes. Several lines of evidence point to an important modulatory role of GABA and glutamate systems on neuroinflammatory processes [14] in the mature brain. This connection to signaling mechanisms of immunity makes them important therapeutic targets to reverse the detrimental effects of chronic neuroinflammation, as it is pivotal to the pathophysiology of numerous diseases of the central nervous system. In the developing brain, both molecules modulate neural precursor cell proliferation, differentiation and neuron migration through a reciprocal relationship with neurotransmitter-sensitive immune cells such as microglia [15,16]. Signaling mechanisms involved in these processes are thought to be operational in glioma [17].
Starting with the postulate that neurotransmission signaling activity is of importance to glioma biology, we initially interrogated the transcriptome of 661 gliomas available through the TCGA database for GABA-, glutamate-and calcium-related gene expression patterns. Four glioma clusters were generated by unsupervised clustering and compared to established glioma molecular diagnostic subgroups, as per the World Health Organization glioma classification. We subsequently used theses analyses as a starting point in deciphering ties to epigenetic (DNA methylation and microRNA) and immune mechanisms to find novel vulnerabilities.

Sample Extraction
Data from the low-grade glioma [18] (LGG) and glioblastoma multiforme [19] (GBM) projects were extracted from the TCGA (The Cancer Genome Atlas, https://portal.gdc. cancer.gov/, accessed on 1 September 2020). These datasets include clinical and epidemiological data, gene and miRNA expression, methylation quantification, copy number variation and simple nucleotide variation. Biomolecular information (TERT promoter, ATRX, MGMT promoter status) extracted from the Ceccarelli study [20] was also added to our analysis.
Two RNA-seq datasets were extracted from the Chinese Glioma Genome Atlas [21][22][23] (CGGA, http://cgga.org.cn/, accessed on 1 September 2020) and were merged into one dataset to simplify the analysis. Clinical data (sex, age or overall survival) associated with these samples were also extracted.

Gene Expression Normalization
Expression normalization was performed to enable comparisons between glioma samples. TCGA raw-count expression data were normalized using the Variance Stabilizing Transformation function available in the DESeq2 [25] v1.26.0 R package. Low-expression genes whose maximum did not pass 10 counts were excluded.

Unsupervised Clustering
Unsupervised clustering was performed to reveal hidden patterns of neurotransmissionrelated gene expression. Entropy values were calculated, varying the cluster number in order to select the optimal number of clusters that minimized the heterogeneity of IDH mutation and 1p19q codeletion within a cluster. Both unsupervised clustering and heatmaps were generated using the Complex Heatmap [26] v2.4.2 R package, based on the Ward method and Spearman correlation as distance.

Differential Gene Expression Analysis
The DESeq2 [25] R package v1.26.0 was used to retrieve differentially expressed genes between clusters. Bonferroni-corrected p-values below 0.05 were considered significant. Volcano plots of differential gene expression were generated using the EnhancedVolcano v1.8.0 R package.

miRNAs Interactome Profiling
MiRNAs interacting with the neurotransmission-related gene set were extracted from the RNA Interactome database [27]. Using the TCGA miRNA dataset, differential gene expression analysis was performed to filter differentially expressed miRNAs between identified glioma clusters (DESeq2 Bonferroni-adjusted p < 0.001 and absolute log2 fold change superior to 1). MiRNAs with an average expression lower than 1 RPKM were not included in the analysis.

Snakemake Pipepeline Creation
The study pipeline was built using snakemake [28] workflow manager v5.32.0. The software and packages used in the pipeline were downloaded from the bioconda channel via the package manager conda, and steps with high computational cost were executed using Compute Canada structures. This pipeline is available in GitHub at this address: https://github.com/hoang31/gaba_glutamatate_TCGA_profiling.git (accessed on 1 September 2020).

Unsupervised Clustering Based on GABA, Glutamate and Calcium Gene Expression Distinguishes Four Clusters with Distinct Neurotransmission Profiles
A total of 421 neurotransmission-related genes were extracted from the KEGG metabolic database, more specifically from four metabolic pathways: GABAergic neuron (kegg id: 04727), glutamatergic neuron (kegg id: 04724), glutamate metabolism (kegg id: 00250) and calcium signaling and endocrine (kegg ids: 04020 and 04961) pathways. Of these, 351 minimally expressed genes were kept (Supplementary Data Table S1), and there was limited overlap between the different gene sets (Supplementary Data Figure S1). Unsupervised clustering of 661 TCGA glioma samples was then performed in order to evaluate the expression pattern of these genes (Figure 1). IDH and 1p/19q codeletion-based entropy analysis showed an optimal number of clusters equal to four (Supplementary Data Figure S2). We thus generated four clusters, which were renamed NT-1 (n = 168), NT-2 (n = 188), NT-3 (n = 81) and NT-4 (n = 224) for neurotransmission-related clusters (Supplementary Data Table S2). This clustering analysis identified four main clusters that differ greatly in their neurotransmission-related gene expression profiles. Figure 1. GABA, glutamate and calcium pathway-related gene expression signatures of 661 TCGA gliomas samples. Four clusters were generated by hierarchical unsupervised clustering using Pearson correlation (glioma samples in column) and Minkowski distance (genes in rows). Ward.D2 clustering method was selected for this clustering. Colors in rows represent different KEGG metabolic pathways of the included neurotransmission genes.

Clinical and Epidemiological Characterization of Neurotransmission-Based Glioma Clusters
We evaluated the clinical and histopathological characteristics of the four glioma subgroups identified based on their expression patterns of glutamate, GABA and calcium genes. Comparison analysis (Table 1) did not show significant age differences between the NT-1, NT-3 and NT-4 gliomas (average of 42.69 yo, p = 0.08 and 0.18). On the other hand, NT-2 gliomas affected significantly older patients (average of 58.12 yo, p = 1.06 × 10 −21 ; p = 4.29 × 10 −12 ; p = 9.95 × 10 −29 ). Gliomas were more frequent in males than females, with no significant gender distribution differences between the different clusters. Survival analysis showed similar survival for NT-1 and NT-3 (log-rank test p = 0.44, Figure 3A). However, survival was shorter for NT-2 gliomas (log-rank test, p = 7.10 × 10 −32 , p = 1.06 × 10 −20 and 1.09 × 10 −33 compared to NT-1, NT-3 and NT-4, respectively). The same pattern was observed with regard to the Karnofsky's performance scores associated with the different glioma NT clusters; NT-2 glioma patients were associated with significantly lower performance scores (Fisher's exact test, p = 5.00 × 10 −4 Figure 3B).  Red, blue, green and purple represent the various NT-1-4 glioma clusters, respectively. Cluster-to-cluster significance was calculated using the log-rank test. (B) Karnofsky performance score distribution among NT-1-4 gliomas. These scores reflect the patient's ability to perform ordinary tasks and range from 100 (patient without disabilities) to 0 (patient death).
We performed a univariate and multivariate cox regression analysis to validate the prognosis associated with each glioma cluster ( Table 2). As expected from the literature, univariate regression suggested that higher age at diagnosis, clustering within the NT-2 cluster or higher grade had significant negative impact on patient survival (beta = 0.066 with p = 9.31 × 10 −42 ; beta = 2.216 with p = 6.62 × 10 −28 ; beta = 2.976 with p = 6.38 × 10 −43 for NT-2 cluster, G4 glioma and age at diagnosis, respectively). The multivariate cox regression model confirmed that clustering within the NT-2 glioma cluster is an independent factor impacting patient survival (beta = 0.901; p = 1.93 × 10 −4 ) regardless of gender, age or grade. As expected, age at diagnosis and grade were also found to be independent, significant prognostic factors for gliomas. NT-2 gliomas, which are mostly comprised of IDH-wt gliomas, bear the worst prognosis out of the neurotransmission-related clusters. NT-1 gliomas had similar prognosis compared to NT-3 and NT-4 gliomas.

Lower Expression of Neurotransmission Genes Correlates with Increased Aggressivity in the NT-1, NT-2, NT-3 and NT-4 Gliomas
Following the identification of four glioma clusters with distinct GABA, glutamate and calcium-related gene expression patterns, we searched for specific genes expressed by cancer cells or their microenvironment that may impact gliomagenesis by performing differential gene cluster-to-cluster expression analysis. We retrieved 232 unique neurotransmissionassociated genes that were significantly expressed (Bonferroni-adjusted p < 0.001 and absolute value of log2foldchange greater than 1). In this analysis, the highest number of significantly differentially expressed genes (43 genes) belonged to the NT-1 cluster comparison ( Figure 4A).
Cluster-to-cluster expression profile analyses for calcium endocrine, calcium signaling, GABA synapse, glutamate metabolism and glutamate synapse metabolic pathways are presented in Figure 4B and Supplementary Data Table S4. Amongst all samples (glioma and healthy patients), we observed that neurotransmission-related genes were the most highly expressed in healthy samples when compared to NT glioma clusters. (18/23, 82/132, 43/55, 6/15 and 49/65 genes for calcium endocrine, calcium signaling, GABA synapse, glutamate metabolism and glutamate synapse metabolic pathways, respectively). With regard to the comparisons between NT glioma clusters, we found that NT-1 gliomas were associated with the largest number of overexpressed genes (15/23, 70/132, 41/55, 8/15 and 43/65 genes for calcium endocrine, calcium signaling, GABA synapse, glutamate metabolism and glutamate synapse metabolic pathways, respectively). Conversely, NT-2 gliomas were associated with the largest number of underexpressed genes ascribed to these pathways (14/23, 60/132, 36/55, 7/15 and 48/65 for calcium endocrine, calcium signaling, GABA synapse, glutamate metabolism and glutamate synapse metabolic pathways, respectively). These observations were more significant for GABA synapse, glutamate metabolism and glutamate synapse pathways. NT-3 and NT-4 cluster expression profiles showed intermediate levels of expression for all pathways. We also identified three genes that were significantly differentially expressed amongst all six cluster-to-cluster analyses: the cholinergic receptor muscarinic 1 (CHRM1), the cholinergic receptor muscarinic 3 (CHRM3) and the glutamate ionotropic receptor NMDA type subunit 1 (GRIN1) genes. These three genes were significantly overexpressed in healthy samples and NT-1 gliomas ( Figure 4C). Overall, we found that NT-1 gliomas are characterized by the overexpression of neurotransmission-related genes when compared to other glioma clusters.

Correlation between DNA Hypermethylation and Gene Expression Is Preserved in NT-1 Gliomas
We then sought to evaluate the role of DNA methylation, which is an important epigenetic mechanism involved in gliomagenesis, in modulating neurotransmission gene expression in NT gliomas. Using the TCGA methylation quantification dataset, we first examined the DNA methylation levels of GABA-, glutamate-and calcium-related genes used for NT-1-4 glioma unsupervised clustering. For this, we extracted the methylation beta values of the differentially expressed neurotransmission-related genes between NT-1-4 gliomas (Wilcoxon rank sum test Bonferroni-corrected p-value less than 0.001 and absolute log2 fold change superior to 1). Methylation beta values equal to 0 and equal to 1 reflect DNA hypomethylation and hypermethylation, respectively. We found that NT-2 gliomas were associated with lower average beta values than NT-1, NT-3 and NT-4 gliomas ( Figure 5A, Wilcoxon rank sum test adjusted by Bonferroni correction; p = 4.9 × 10 −10 ; p = 5.8 × 10 −14 ; p = 1.1 × 10 −5 , respectively). The average beta values of NT-2 were also significantly higher than in healthy samples albeit with lower statistical significance because of the number of healthy samples in the analysis (n = 2) (Wilcoxon rank sum test p = 1.2 × 10 −9 ). We further investigated DNA methylation levels for individual genes ascribed to calcium endocrine, calcium signaling, GABA synapse, glutamate metabolism and glutamate synapse signaling pathways by generating methylation profiles ( Figure 5B). Again, NT-2 gliomas and healthy samples were associated with overall DNA hypomethylation. Interestingly, a significant negative correlation between DNA methylation and expression levels was only maintained for NT-1 gliomas (r = −0.53; p = 3.02 × 10 −13 ; Figure 4D). However, DNA methylation levels in NT-2, NT-3 and NT-4 gliomas correlated weakly with gene expression levels (r = −0.23; p = 1.21 × 10 −5 ; Figure 5C).

NT-1 and NT-2 Gliomas Are Regulated by More Complex Epigenetic Mechanisms Involving Differential Expression of microRNA
MicroRNAs (miRNAs) are small, single-stranded, non-coding RNA molecules (21-25 nucleotides in length) that play an important role in tumorigenesis through RNA silencing and post-transcriptional regulation of gene expression [31]. We explored their potential regulatory roles on GABA-, glutamate-and calcium-associated gene expression in gliomas by extracting from the RNA Interactome database [27] the miRNAs that interacted directly with our 351 neurotransmission-related genes. We identified 73 differentially expressed miRNAs (DE-miRNAs) in this analysis. The largest counts of highest-expressed DE-miRNAs were found in NT-1 and NT-2 gliomas (26, 25, 14 and 8 DE-miRNAs for NT-1, NT-2, NT-3 and NT-4, respectively; Figure 6A) and the largest counts of lowest-expressed DE-miRNAs were found in NT-2 and NT-3 gliomas (31, 25, 9 and 8 for NT-2, NT-3, NT-1 and NT-4, respectively).  Healthy samples (n = 2) were also added to the analysis. (C) Correlation between average neurotransmission-related gene DNA hypermethylation and expression in NT-1-4 gliomas (n = 661). Each point represents a glioma sample, and the color is specifically related to the NT glioma cluster (*** p < 0.001).  A total of 29 DE-miRNAs were found to be associated with higher expression compared to the average expression of all DE-miRNAs (average of 4.3 DESeq2 normalized counts; Supplementary Data Table S5). Amongst these 29, the top 4 were the miRNAs hsa-mir-100, hsa-mir-183, hsa-mir-128-2 and hsa-mir-23a ( Figure 6B and Supplementary Data Table S5). The hsa-mir-100 and hsa-mir-23a mirRNAs were significantly overexpressed in NT-2 gliomas when compared to NT-1, NT-3, and NT-4 clusters (Bonferroni-adjusted p < 0.05 for all cluster pairwise comparison). The hsa-mir-128-2 was overexpressed in NT-1 (Bonferroniadjusted p < 10 × 10 −5 for all cluster pairwise comparison). Expression of hsa-mir-183 was lower in NT-2 gliomas (Bonferroni-adjusted p < 10 × 10 −5 for all cluster pairwise comparison). Overall, neurotransmitter-related miRNAs appeared to be more deregulated in NT-2 gliomas when compared to other glioma subgroups.

Neurotransmission-Related Gene Expression Correlates with the Immune Response Signaling Pathways in NT-1-4 Glioma Clusters
To look into the regulation of specific cellular signaling pathways by neurotransmitterrelated genes in glioma, we performed a correlation analysis of the TCGA glioma dataset between previously-identified NT1-4 discriminatory genes (CHRM1, CHRM3 and GRIN1) and non-neurotransmission-related genes. We retrieved 7758 and 7346 genes with negative and positive correlation, respectively (Bonferroni-adjusted p < 0.05, Supplementary Data Table S6).
Analyses of negatively-correlated genes revealed enrichment for genes pertaining to 15 specific cellular pathways such as cell activation, immune effector process, cell population proliferation, response to biotic stimulus, primary metabolic process, symbiotic process, movement of cell or subcellular component, response to external stimulus, etc. (Supplementary Data Figure S3). Amongst these pathways, 10 out of 15, such as cell activation, cell population proliferation and movement of cell and/or subcellular component process groups, were mainly composed of immune cellular processes such as positive regulation of leukocyte activation, T cell activation and leukocyte proliferation (Supplementary Data Figure S4).
Positively-correlated genes were significantly associated with signaling pathways such as system process, macromolecule localization, establishment of localization, regulation of biological quality, cellular component organization or biogenesis and cell communication (Supplementary data Figure S5). In summary, neurotransmission gene expression in gliomas correlates with various cellular pathways related to immunity.

Immune Cell Characterization Reveals Different Tumor Immune Microenvironment Composition
The identification of a high number of immune processes with the CHRM1, CHRM3 and GRIN-1 neurotransmission-related gene signature may suggest that NT-related gliomas possess distinct tumor immune microenvironments. We first performed an ESTIMATE R tumor purity calculation on NT-1-4 gliomas, as this tool evaluates immune and stromal gene expression signatures. We found that the NT-1 gliomas were associated with significantly higher tumor purity scores when compared to the NT-2 (Wilcoxon rank sum test p = 3.01 × 10 −45 ) and NT-4 (Wilcoxon rank sum test p = 5.2 × 10 −26 , Figure 7A) gliomas, reflecting lower levels of immune cells in the NT-1 gliomas. There were no significant differences observed between NT-1 and NT-3 gliomas.  To further substantiate this observation, we also inferred the immune cell type composition of each glioma cluster using CIBERSORT [32] and CIBERSORTx [33] tools ( Figure 7B). The CIBERSORTx analysis showed similar absolute immune cell quantification between the NT-1 and NT-4 gliomas (Wilcoxon rank sum test p = 0.28). NT-2 gliomas were associated with a higher number of immune cells (Wilcoxon rank sum test p = 6.42 × 10 −32 , p = 1.94 × 10 −27 and p = 1.70 × 10 −26 compared to NT-1, NT-3 and NT4, respectively). When analyzing the CIBERSORT results, NT-2, NT-3 and NT-4 gliomas were significantly associated with a higher fraction of M2-phenotype macrophages when compared to NT-1 gliomas (Wilcoxon rank sum test p = 3.90 × 10 −30 , p = 2.01 × 10 −8 and p = 1.49 × 10 −17 , respectively). In addition, NT-1 gliomas had a higher fraction of plasma B cells (Wilcoxon rank sum test p = 1.73 × 10 −42 , p = 2.42 × 10 −9 and p = 2.33 × 10 −35 , respectively) when compared to the three other groups. NT-1 gliomas also had a higher fraction of monocytes when compared to NT-2 and NT-3 gliomas (Wilcoxon rank sum test p = 7.45 × 10 −10 , p = 1.23 × 10 −5 and p = 8.36 × 10 −1 , respectively). Other immune cell types were identified to be significantly enriched in NT-1 gliomas and are described in the Supplementary Data Table S7.
Principal component analysis (PCA) on the CIBERSORTx immune data showed that the first and second components explained 82.3% of the total variability (75.1% and 7.2% for the first and second component, respectively, Figure 7C). The first principal component mainly consists of the M2 macrophage variable (93.87%). Monocytes were the main contributors to the second principal component (59.89%), with macrophages M0 (14.06%), mast cells resting (12.15%) and B cells plasma (8.01%) also contributing.

Neurotransmission-Related Transcriptomic Profiling on the Chinese Glioma Genome Atlas Cohort
We tested the reproducibility of our findings using 889 glioma samples from the Chinese Glioma Genome Atlas cohort (CGGA) with the neurotransmission-related gene set (351 genes) used in the analysis of the TCGA cohort. We performed unsupervised clustering and generated four different clusters, which were identified as NT-1-like (n = 189), NT-2-like (n = 259), NT-3-like (n = 166) and NT-4-like (n = 275) glioma clusters ( Figure 8A). Similar to the TCGA cohort, the NT-2-like cluster was mainly composed of IDH-wt gliomas (222/259 or 85.71%) whereas IDH-mutated gliomas were predominant in the NT-3-and NT-4-like clusters (122/166, 73.49% and 226/275, 82.18% for the NT-3-and NT-4-like clusters, respectively; Supplementary Data Table S8). The NT-1-like glioma cluster was composed of a mixture of IDH-mutated (112/189 or 59.26%) and IDH-wt gliomas (77/189 or 40.74%). The 1p/19q chromosomal co-deletion did not predominate in any cluster.  The four clusters were generated by hierarchical unsupervised clustering using Pearson correlation (glioma samples in column) and Minkowski distance (genes in rows). Ward.D2 clustering method was selected for this clustering. (B) Kaplan-Meier survival curves associated with the CGGA NT-1-4 like glioma clusters. Red, blue, green and purple represent the various CGGA NT-1-4-like glioma clusters, respectively. Cluster-to-cluster significance was calculated using the log-rank test. (C) Immune cell type inference for NT-1-4 like gliomas using CIBERSORT. CIBERSORT infers the immune cell fraction relative to the total immune cell content.
In terms of immune cell type composition, NT-2-like gliomas were associated with a higher cell fraction of M2-phenotype macrophages, similar to the TCGA cohort analysis findings. (Wilcoxon rank sum test p = 4.76 × 10 −22 , p = 1.03 × 10 −29 and p = 2.09 × 10 −17 compared to the NT-1-like, NT-3-like and NT4-like, respectively; Figure 8C). In general terms and similar to the TCGA analyses, neurotransmission-based unsupervised clustering of the CGGA glioma expression dataset recapitulated four subgroups with similar survival and inflammatory microenvironment characteristics.

Discussion
Current glioma therapies lead to limited improvement in median overall survival in patients with high-grade infiltrating gliomas [34]. Beyond targeting tumor cells, there is currently a shift of focus towards understanding components of the tumor microenvironment, such as surrounding neural cells (neurons and glia), hematopoietic cells (monocyte/macrophage/microglia and T cells) or blood vessels, in an attempt to overcome redundant compensatory mechanisms [35]. In this large-scale multi-omics analysis, IDH-wildtype and IDH-mutated infiltrating gliomas were studied through the lens of neurotransmissionrelated (GABA, glutamate and calcium) gene expression patterns with the aim of unraveling specific vulnerabilities and cellular pathways.
Neurotransmission-based unsupervised clustering enabled the proper classification of the majority of infiltrating gliomas into current WHO tumor categories (IDH-wt gliomas, IDH-mutated, 1p/19q oligodendrogliomas and IDH-mutated astrocytomas), suggesting that neurotransmission-related pathways are differentially regulated in tumor cells and/or their microenvironment according to tumor subtype, and also reaffirming the importance of IDH1/2 mutations and of 1p/19q-codeletion for glioma stratification. Interestingly, this strategy also identified a novel NT-1 glioma subgroup mostly comprised of IDH-mutated gliomas, which also included "normal-like" IDH-wt gliomas associated with a longer survival [1,29].
The NT-1 subgroup primarily distinguishes itself by its overexpression of GABA, glutamate and calcium genes. Conversely, the NT-2 cluster, which is principally comprised of more-aggressive IDH-wildtype glioblastomas, is associated with a lower expression of neurotransmission-related genes. IDH-mutated gliomas with (NT-3) or without (NT-4) the 1p/19q codeletion show intermediate levels of expression. This supports the importance of the IDH mutation in regulating neurotransmission-related gene programs in glioma [36,37].
The role of the IDH mutation as a trigger of the CpG island methylator phenotype in IDH-mutated gliomas via the production of 2-hydroxyglutarate (2-HG) oncometabolite is well-established [38,39]. It is interesting to note that-while the average methylation of neurotransmission genes in NT-1 gliomas is high, as expected from a group mainly composed of IDH-mutated tumors-these tumors are distinct from other IDH-mutated gliomas by the partial preservation of their capacity to silence neurotransmission genes through methylation. We can speculate that epigenetic DNA methylation events that follow IDH mutation in early gliomagenesis target neurotransmission genes randomly, thereby accounting for heterogeneous neurotransmission-related profiles and selective vulnerabilities within IDH-mutated tumors. Neurotransmitter-related genes amenable to epigenetic reprogramming may impact the therapeutic efficacy of experimental anti-cancer DNA demethylating drugs [40].
Correlations between the CHRM1, CHRM3 and GRIN-1 NT-1-4 glioma intersecting gene signatures and various immune signaling pathways suggested that NT gliomas may be endowed with distinct tumor immune microenvironments. Further, 2-HG is an important mediator of tumor immunity in IDH-mutated gliomas. It acts as a suppressor of antitumor T-cell activity and also impedes macrophage recruitment in gliomas by altering tryptophan metabolism [59,60]. Lower pools of M2-phenotype macrophages were detected in NT-1 gliomas when compared to NT-2-4 gliomas. Considering the anti-tumor and immunosuppressive functions of this type of macrophage in gliomas, this observation could be in keeping with a less suppressive and less tumor-supportive inflammatory microenvironment in NT-1 gliomas, partly explained by their IDH status but also promoted by operational neurotransmitter signaling pathways modulated by released GABA or glutamate in cancer or microenvironment cells [61].
We identified novel, neurotransmission-based glioma subgroups with their peculiarities in terms of expression, epigenetics (methylation and miRNAs) and inflammatory microenvironment. These findings may be of clinical relevance should neurotransmitter pathways impacting tumor aggressiveness be actionable or reactivable in gliomas, whether it is through epigenetic-based strategies for IDH-mutated tumors or any other strategy for IDH-wt gliomas.
The adult gliomas included in this study are infiltrative by nature. It is thus expected that brain biopsy samples will include non-neoplastic brain cells. We do expect that some of the differentially expressed NT genes identified in this study belong to the non-neoplastic tumor microenvironment, where they can still impact tumor aggressiveness. Further bioinformatics studies on microdissected tumor tissue and single-cell RNA sequencing data will help to specify more precisely the location of neurotransmission targets in tumor cells versus non-neoplastic tumor cells of the microenvironment. Experimental studies targeting neurotransmitter signaling elements and relevant miRNAs in vitro will also further our knowledge regarding the impact of these pathways on glioma aggressiveness.

Conclusions
This multi-omics analysis revealed the existence of neurotransmission-based glioma categories with significant differences in regard to neurotransmission-related gene expression, methylation, and miRNA profiles in adult gliomas. It also revealed alterations in the nature of the tumor inflammatory microenvironment between NT glioma subgroups. Deciphering operational neurotransmitter signaling pathways and underpinning mechanisms that may represent actionable targets is a promising personalized treatment avenue to explore for glioma patients as a complement to current radiotherapy and chemotherapy treatments in an attempt to improve clinical outcomes.
Supplementary Materials: The following are available online at https://www.mdpi.com/article/10.3 390/jpm12040633/s1. Supplementary Data Figure S1: Upset plot describing the 351 neurotransmissionrelated genes used in our analyses. These gene sets were taken from the KEGG metabolic database, Supplementary data Table S1: List of 351 neurotransmitter-related genes used for unsupervised clustering (>10 counts), Supplementary Data Figure S2: Average entropy of the TCGA glioma expression dataset unsupervised clustering based on the IDH and 1p/19q codeletion status, according to the number of clusters generated. The optimal number of clusters is equal to 4 (red dot), Supplementary Data Table S2: Distribution of the NT-1-4 glioma clusters among the 661 TCGA glioma samples, Supplementary Data Figure S3: Gene ontology enrichment analysis on correlated genes between previously-identified NT1-4 intersecting genes (CHRM1, CHRM3 and GRIN1) and non-neurotransmission-related genes., Supplementary Data Table S3: Distribution of biomolecular and epidemiological variables among NT-1-4 glioma clusters, Supplementary Data Figure S4: Significant Gene Ontology group terms associated with cell activation, cell proliferation and movement of cell or subcellular component gene ontology level 3 gene ontology group terms. Colors represent the correlation type (blue and red for negative and positive correlations, respectively). Supplementary Data Table S4: Distribution of overexpressed and underexpressed neurotransmission-related genes in NT-1-4 gliomas, Supplementary Data Figure S5. Significant Gene Ontology group terms associated with system process, macromolecule localization, establishment of localization, regulation of biological quality, cellular component organization or biogenesis and cell communication gene level 3 gene ontology group terms. Colors represent the correlation type (blue and red for negative and positive correlations, respectively), Supplementary Data Table S5: Average expression levels of 73 differentially expressed miRNAs in NT-1-4 gliomas. Supplementary Data Table S6: Positive and negative correlations between CHRM1, CHRM3 and GRIN-1 genes and nonrelated neurotransmission genes. Supplementary Data Table S7: Comparison of immune cell type enrichment among NT-1-4 glioma clusters generated by the CIBERSORT analysis tool. Supplementary Data Table S8

Institutional Review Board Statement: Not applicable.
Informed Consent Statement: Not applicable.

Data Availability Statement:
The data used in the study are available in the TCGA public databases. The datasets supporting the conclusions of this article are included within the article (and its additional files).