Abstract
Gliomas are the most common primary brain cancers. In recent years, IDH mutation and 1p/19q codeletion have been suggested as biomarkers for the diagnosis, treatment and prognosis of gliomas. However, these biomarkers are only effective for a part of glioma patients and thus more biomarkers are still emergently needed. Recently, an electrochemical communication between normal neurons and glioma cells by neuro-glioma synapse has been reported. Moreover, it was discovered that breast-to-brain metastasis tumor cells have pseudo synapses with neurons and these synapses were indicated to promote tumor progression and metastasis. Based on the above observations, we first curated a panel of 66 SA genes and then proposed a metric, SA score, to quantify the synapseness for each sample of 12 glioma gene expression datasets from TCGA, CGGA, and GEO. Strikingly, SA score showed excellent predictive ability for the prognosis, diagnosis, and grading of gliomas. Moreover, being compared with the two established biomarkers, IDH mutation and 1p/19q codeletion, SA score was demonstrated independent and better predictive performance. In conclusion, this study revealed that SA genes contribute to glioma formation and development, and proposed a quantitative method, SA score, as an efficient biomarker for monitoring gliomas.
Introduction
Brain and other nervous system cancers are estimated to take up 1.4% of new cancers but 2.9% of cancer deaths in 20191. Gliomas are the most frequent cancers of these cancers, including astrocytoma (including glioblastoma), oligodendroglioma, ependymoma, oligoastrocytoma (mixed glioma), malignant glioma, not otherwise specified (NOS), and a few rare histologies2. World Health Organization (WHO) classified gliomas into grades I to IV, and introduced biomarkers of IDH mutation and 1p/19q codeletion in 2016 edition3,4. Glioblastoma (WHO grade IV) accounts for about half of gliomas, with a median survival of less than two years2,5. Gliomas with lower grade have a diverse prognosis, either progressing to be as poor as glioblastoma or living more than 10 years after effective treatment6.
Over the years, with the fast improvement of omics and big data technology, RNA sequencing has been developing towards lower cost and higher throughput, producing a large amount of biological and medical data, which provides great convenience for life science research7,8. Impelled by advantage of big data analysis, numerous biomarkers have been found in the diagnosis and prognosis of gliomas9,10. In addition, single-cell sequencing showed more detailed landscapes of gliomas11-13. Gene set enrichment analysis (GSEA) provides a facility to extract effective information from a large number of RNA expression data14. Moreover, single sample GSEA (ssGSEA) can calculate without group information and give every sample an enrichment score 15. The Biomarkers such as IDH mutation and 1p/19q codeletion provided helps for monitoring the development and prognosis of gliomas but are only effective for a part of patients16. Therefore, given the enormous severity of gliomas, more biomarkers are emergently needed.
It is recently reported that neuron and glioma have electrochemical communication through AMPA receptor-dependent synapses between presynaptic neurons and postsynaptic glioma cells17,18. These observations suggest that the neural synaptic electrochemical connections promote glioma progression. Simultaneously, an appearance of glutamatergic ‘pseudo-tripartite’ synapses between breast-to-brain metastasis tumor cells and neurons was observed19. Based on these anatomical and cytological findings, we hypothesized that the synapse assembly (SA) genes can be used as a biomarker for glioma prognosis. To confirm this hypothesis, here we first curated a list of genes involved in SA and then performed ssGSEA analysis for glioma gene expression datasets from the Cancer Genome Atlas (TCGA), the Chinese Glioma Genome Atlas (CGGA), and the Gene Expression Omnibus (GEO). Interestingly, the tumor synapses may be involved in immune escape and multiple molecular pathways. Strikingly, SA was found to be an independent and effective biomarker for gliomas.
Results
The panel of synapse assembly genes serves as a novel biomarker for gliomas
In order to investigate whether SA can be biomarker for glioma patients, we first curated a list of SA genes. As a result, 66 genes were collected (Supplementary File 1). Then we performed ssGSEA to TCGA and CGGA datasets for survival analysis. The results showed glioma patients with higher SA scores have longer overall survival time (Fig 1). Cox regression analysis also shows the same results (Table 1). Moreover, patients with higher WHO grade have significantly lower SA scores (Fig. 2), which agrees with the survival analysis. In addition, it is deserved to be mentioned that SA score shows an ability to distinguish between glioma and normal brain tissue, which reveals a potential diagnostic application of SA (Fig. 2e-g, Supplementary Fig. 1; AUCs of ROC analysis for GSE4290, GSE16011, and GSE50161 datasets are 0.91, 0.99, and 0.93, respectively).
Comparison of SA score with established biomarkers
IDH mutation and 1p/19q codeletion are two established biomarkers for gliomas. Both biomarkers provided great helps for monitoring glioma development but both are effective on only a part of patients. Therefore, it is interesting to explore whether SA score is an independent biomarker and whether SA score is better than the established biomarkers or not. For doing so, we first analyzed the relationship of SA scores with IDH mutation and 1p/19q codeletion status. We found that IDH-mut gliomas were associated with significantly higher SA scores than IDH-wt ones (Fig. 3a-d). And 1p/19q codeletion gliomas represent a higher SA scores than non-codeletion ones (Fig. 3e-h). Moreover, after moving the effects of the two established biomarkers using multivariate Cox regression model, we revealed that SA score is an independent biomarker for predicting prolonged overall survival in gliomas (Table 1). In addition, the grading ability of SA score is also independent to IDH mutation and 1p/19q codeletion (Supplementary Fig. 2). Finally, we compared the predictive performance of SA score, IDH mutation and 1p/19q codeletion status (Table 2). In most instances, SA score outperforms IDH mutation and 1p/19q codeletion.
Discussion
Given the recently revealed roles of neuro-glioma synapse in glioma development, here we curated a panel of 66 synapse assembly (SA) genes and proposed the SA score as a biomarker for the prognosis, grading, and diagnosis of gliomas. The SA score was validated by over 3000 samples of 12 datasets from TCGA, CGGA, and GEO.
To confirm whether the SA score does represent the process of synapse assembly or not, we performed Spearman’s correlation analysis between SA score and expressions of synapse receptor genes 17 (Supplementary Fig. 3a). The result showed that a number of genes, especially AMPA (α-amino-3-hydroxy-5-methyl-4-isoxazole propionic acid) receptor genes (not contain in SA genes), which were focused in recent studies17,18, have a positive correlation with SA scores. Then, we applied pathway analysis to genes with significant correlation with SA scores (Methods). As a result, the positively correlated genes were significantly enriched in synaptic related pathways (Supplementary Fig. 3b, c) and cell adhesion molecules (CAMs, hsa04514) in neural system (Supplementary Fig. 4), suggesting that the calculated SA score can represent the process of synapse assembly.
In spite of its ability as glioma biomarker for the identified SA gene panel, it should be especially noted that the result seems opposite to existed knowledge. That is, it was reported that neuro-glioma synapse could promote tumor progression and metastasis17-19, which thus can infer that SA should result in a poorer prognosis, however, we revealed it is associated with a better but not poorer prognosis. Molecular processes may play different roles in various cells, organs and diseases. For example, as an important discovery in glioma research, IDH mutation is identified as one of the early events of glioma, and the epigenetic changes caused by IDH mutation are considered as a main tumor driver 20. Nevertheless, clinical studies have found that IDH mutation can lead to a longer survival time 21. Similarly, immunotherapy, which has been widely used, was criticized for producing serious side effects 22,23. These instances suggest that the SA gene panel could also have multiple aspects.
To better understanding the molecular processes of the SA genes, we sought to pathway analysis. Besides the proliferative effect suggested in previous studies, we revealed that neuro-glioma synapses may also show a participation in the immune response, cyclic adenosine monophosphate (cAMP, hsa04024) and nuclear factor kappa-B (NF-κB, hsa04064) pathway (Supplementary Fig. 3b, c). The SA genes are related to the activation of cAMP pathway (Supplementary Fig. 3c, 5), which has been proved to be involved in the synaptic plasticity process, including AMPA receptor trafficking, which may be a reason why glioma cells form synapses and produce AMPA receptors24. Moreover, cAMP can also activate phosphoinositide 3-kinases (PI3K)/protein kinase B (PKB, also known as Akt), thus promoting cell growth, survival and proliferation, which may be one of the explanations for the proliferation promotion effect of the newly discovered neuro-glioma synapse 25. Since the synaptic structure was also observed in the breast-to-brain cancer cells, we hypothesized that it was not the cell origin but the environment that determined this morphological change in tumor. It is worth mentioning that the cAMP pathway can be activated by neurotransmitters and Ca2+, which is extensively distributed in synaptic cleft and easily accessible for glioma26-28.
The NF-κB signaling pathway and CAMs of immune system were also found to be vastly enriched with genes negatively associated with SA (Supplementary Fig 3c, 4, 6). Previous studies on NF-κB are mainly focused on immune cells, which are involved in cell activation and survival, as well as the production of cytokines29,30. In tumors, the stimulation of NF-κB signal leads to up-regulation of anti-apoptotic genes that allow tumors to tolerate inflammatory environment. Meanwhile, NF-κB induces cytokines and CAMs, which give rise to the recruitment of immune cells to the tumor site31. As can be seen, the relationship between NF-κB-caused inflammation and cancer is complex, which can either eliminate cancer cells or lead to immune escape. Additionally, stimulated cAMP could suppress NF-κB activation 32. This finding, together with the relationship between NF-κB and CAMs, is in good agreement with our discoveries. We also demonstrated that SA is associated with IDH mutation status, which is a momentous survival risk factor. In a previous study, significantly higher tumor infiltrating lymphocytes (TILs) infiltration and programmed death ligand 1 (PD-L1) expression, which are targets of immune checkpoint inhibitor, were found to be associated with IDH-wt33. However, the others have found that IDH-mut tumors show reductions of various immune cells, including microglia, macrophages and polymorphonuclear leukocytes, which are consistent with the decrease of chemokines in tumors34.
As it was shown in recent studies17-19, the neuro-glioma synapses are involved in cancer proliferation and metastasis, which are beneficial to tumor progression. But reversely, we found that higher SA signifies lower glioma grade and longer overall survival. Analogously, IDH-mut and 1p/19q codeletion are typically biomarkers that promote glioma progression but benefit prognosis. Existing studies have focused on mechanisms that promote glioma, but the reasons for better prognosis are generally reported by clinical studies, such as better chemoradiotherapy sensitivity35. We conjectured that SA, IDH mutation and 1p/19q codeletion shared a part of the mechanism that resulted in the observed phenomenon. Apoptosis caused by suppression of NF-κB signaling pathway may also be a motivation. The causations in SA, mutation, pathway and immune response are remained to be explored.
In summary, we found that SA is associated with activation of the cAMP signaling pathway and down-regulation of NF-κB pathway and immune-related CAMs. Although the mechanism is unclear, we revealed that SA genes are involved in glioma process and the proposed SA score is an independent and potentially better biomarker for glioma overall survival, and shows a predictive capacity in different grade gliomas and normal brain tissue, which could be useful in the prognosis, grading and diagnosis of gliomas.
Methods
Gene expression datasets and analysis
We obtained sample and clinical information from TCGA (https://portal.gdc.cancer.gov/). RNAseq data of TCGA were downloaded from FireBrowse (http://firebrowse.org/), which were processed and normalized uniformly. Histology type, WHO grade, IDH mutation status and 1p/19q codeletion status were obtained from the study by Ceccarelli et al36. CGGA (http://www.cgga.org.cn/) provides tumor gene expression data for thousands of glioma patients (including one microarray and two RNAseq batches), as well as corresponding clinical data. The calculation and presentation of the results will be conducted separately due to different platforms and batches. In addition, glioma microarray gene expression profiling data (GSE4290, GSE16011, GSE50161, GSE52009, GSE54004, GSE61374, GSE107850) were available at GEO (https://www.ncbi.nlm.nih.gov/gds/). Gene expression data were structured with gene symbols as row names, sample ids as column names, duplicate gene symbols were averaged using their median value.
Single sample gene set enrichment analysis
The human synapse assembly gene set was extracted from the Gene Ontology (GO) by term ‘positive regulation of synapse assembly’ (GO:0051965). The ssGSEA scores were calculated by python (v3.6.8) package gseapy (v0.9.13), which is a python wrapper for GSEA and ssGSEA.
GO and KEGG pathway analysis
We performed Spearman’s correlation analysis between SA ssGSEA score and gene expression in the TCGA low grade glioma (LGG) and CGGA (3 subdatasets) gliomas datasets. Top 5000 (ordered by rho) significantly positively/negatively correlated genes were selected respectively and the results of four datasets were intersected. R package clusterProfiler37 (v3.10.1) was used for GO and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis. P values were adjusted by Benjamini and Hochberg methods and background genes were set as genes which exist in all four datasets. KEGG pathway graphs were rendered by R package pathview38 (v1.22.3).
Statistical analysis
Kaplan-Meier (K-M) curves and Cox proportional hazards regression were performed by R package survival (v2.44-1.1) and survminer (v0.4.6). Log rank test was used to calculate the difference between two K-M curves. Significance of difference between two groups of continuous variables was analyzed by two-side Wilcoxon rank sum test. Receiver operating characteristic (ROC) curve and area under ROC curve (AUROC) were processed by R package pROC39 (v1.15.3). All statistical significances above were calculated by R (v3.5.2). Spearman’s correlation analysis was applied to evaluate the correlation using python package scipy (v1.2.1). P values < 0.05 were considered significant.
Data Availability
The datasets supporting the conclusions of this article are available in TCGA repository, https://portal.gdc.cancer.gov/, the FireBrowse repository, http://firebrowse.org/, the CGGA repository, http://www.cgga.org.cn/ and GEO repository, https://www.ncbi.nlm.nih.gov/gds/. WHO grade, IDH mutation status and 1p/19q codeletion status of TCGA are available in the supplementary file of the published article by Ceccarelli et al, https://doi.org/10.1016/j.cell.2015.12.028. The 66 SA genes are available in the supplementary file 1.
https://portal.gdc.cancer.gov/
Funding
This work was supported by the grants from the Natural Science Foundation of China (81670462, 81970440, and 81921001 to QC).
Author contributions
QC conceived the project. XJ performed the analysis and conducted the experiments. XJ, HZ and QC wrote the manuscript. All authors read and approved the final manuscript.
Supplementary information
Supplementary File 1: List of the 66 SA genes.