Identification of differentially expressed genes and fusion genes associated with malignant progression of spinal cord gliomas by transcriptome analysis

Glioma, the most common histological subtype of primary spinal cord tumors, is considered as a rare central nervous system neoplasm. In this study, 9 glioma samples (4 of grade II and 5 of grade IV with H3K27M positive) were analyzed to examine the molecular mechanisms underlying the malignant progression of gliomas, transcriptome sequencing. Differentially expressed genes (DEGs) in grade IV vs. grade II were analyzed by using the Limma package in R. Enrichment analysis was performed for the individual DEGs through VennPlex software and the Database for Annotation. Gene mutations and fusions were analyzed using the Genome Analysis Toolkit and STAR-Fusion. A total of 416 DEGs were identified in grade IV vs. grade II. Functional analysis of the DEGs showed that GALR1 and GRM5 of neuroactive ligand-receptor interactions signaling pathways may be relaed to malignant progression of gliomas. Further systematic transcriptional profiling identified 11 in-frame/frameshift gene fusions in the tumors. Notably, one novel gene fusions, GATSL2-GTF2I was detected in all of the grade II samples. In summary, the molecular alterations observed in glioma progression may improve the characterization of different human spinal cord glioma grades. The transcriptome analysis of intramedullary spinal cord glioma will provide a new candidate gene list for further mechanism research.

Intramedullary spinal cord gliomas account for about 2 to 4% of central nervous system (CNS) tumors [1][2][3] . Of these spinal gliomas, astrocytomas is the most common subtype. Spinal cord gliomas is classified from grade I (least aggressive) to grade IV (most aggressive) according to conventional histopathological criteria based on the World Health Organization (WHO) grading system 4 . This grading system was dramatically revised in 2016,and the the genetic profiles of tumors 5 was added into it. As one of the deadliest human cancers, grade IV glioma, also termed as glioblastoma (GBM), has a poor median survival time of 12 to 18 months 6,7 . Although the majority of GBMs are primary tumors, around 20% of gliomas still progress from grade II or III 4 . Therefore, the identification of early clinical biomarkers for glioma diagnosis and prognosis is urgently required in clinical practice.
Molecular variations in IDH1, IDH2, H3F3A, HIST1H3B, and BRAF have been frequently identified in different grades of brain gliomas. However, there are very few researches concerning the molecular indicators of spinal gliomas. High-throughout transcriptome sequencing is capable of capturing the transcriptomic landscape of tumors, including protein-coding and non-coding gene expression, the identification of non-coding RNAs (ncRNAs) and/or fusion genes, and the determination of gene mutations and/or alternative splicing 8,9 . The detection of molecular variants across the each grades of spinal gliomas at the transcriptional level will provide further evidence of the mechanisms in the development and progression of spinal glioma.
In this study, transcriptome sequencing of 9 spinal cord glioma samples were performed. The samples included grade II (control) and grade IV tissues. DEGs and gene fusions were identified in different grades of spinal gliomas samples including tissues of grade II (control) and grade IV.he samples included grade II (control) and grade Identification of DEGs. To gain insight into the molecular pathogenesis of spinal gliomas among the Chinese population, we searched for genetic alterations in 9 spinal gliomas by transcriptome analysis.
RNA-Seq analysis included 264 million reads. An average of 29 million short-sequence reads was obtained for each individual sample and more than 81.05% of total reads were identified as clean reads (Supplementary Table S1).
The total clean reads were then mapped to the human whole-genome assembly. The expression level of each gene was normalized and evaluated by the FPKM approach. By comparing the expression of different grades of spinal gliomas, significance differential expression of 416 DEGs were identified in grades IV vs II (Supplementary Table S1).

Identification of neuroactive ligand-receptor interaction signaling pathways. The most enriched
pathways between the astrocytoma grades were identified. A total of 198 differential pathways were identified in IV vs II, including Neuroactive ligand-receptor interaction (11 genes), Cellular senescence (8 genes), cGMP-PKG  www.nature.com/scientificreports www.nature.com/scientificreports/ signaling pathway (7 genes) and so on (See Supplementary Table S3). The gene lists of the Neuroactive ligand-receptor interaction, Cellular senescence, and MAPK signaling pathway were shown in Table 3.
Functional annotations of the DEGs in IV vs II showed that Neuroactive ligand-receptor interaction accumulated the highest number of dysregulated genes, which suggested the close association with spinal glioma progression (Table 3). among the 11 genes, GRM5, GALR1, NPY1R and AGTR1 has been reported playing a role in tumors' development and progression. Therefore, these candidate genes could be used for future research.
Detection of gene fusions. Gene fusions can be identified by searching paired-end reads for two ends mapping to different genes, or for reads containing sequences from two different genes. Of the 44 fusion transcripts from spinal glioma tissues, 11 were in-frame/frameshift and 33 were out of frame (Table 4). Of the in-frame fusions/frameshift, 9 originated from the fusion of sequences located on the same chromosome,and 2 arose from sequences derived from different chromosomes. The majority of gene fusions were from chromosome 1, 7 and 17. 6 out of 11 (54.5%) fusion genes were confirmed by real-time PCR.
Of the 11 fusion genes, GATSL2-GTF2I, were detected in all of four gradeII samples. Other 10 fusion genes were only detected in one sample of grade II or grade IV.

Discussion
In this study, molecular alterations of spinal gliomas were reported by transcriptome analysis, including changes in gene expression, gene fusions, and the mutational landscape of grade II and IV tumors. Totally, 74 DEGs were identified in IV vs II. The identification of the enriched pathways of the DEGs provides information about the cellular processes, which was affected from one spinal glioma grade to the other. Therefore, the most enriched pathways across the astrocytoma grades were compared comprehensively. A total of 187 differential pathways were identified in IV vs II, including neuroactive ligand-receptor interactions (5 genes), hippo signaling (2 genes), Glycosaminoglycan biosynthesis and degradation (10 genes), and pathways in apoptosis (2 genes). These quantitative results suggest that key signaling pathways become increasingly related with the dysregulation of spinal glioma progression. Functional annotations (Supplementary Table S2) of the DEGs in IV vs II showed that neuroactive ligand-receptor interaction pathways accumulated the highest number of dysregulated genes, suggesting their association with spinal glioma progression (  www.nature.com/scientificreports www.nature.com/scientificreports/ in neuroactive ligand-receptor interactions were detected in IV vs II. These quantitative results further support the notion that these pathways become increasingly dysregulated with the progression of glioma malignancy. Neuroactive ligand-receptor signaling pathways accumulated the most DEGs in group's IV vs II, Which suggested the association with spinal glioma progression. Similarily, the enrichment of Neuroactive ligand-receptor interactions was also reported in previous glioma studies [10][11][12] . Among them, GALR1 and GRM5 were further investigated due to their significant dysregulation in IV vs II. Among the relevant candidate genes involved in tumor progression, GALR1 encodes galanin receptor 1 (GALR1), which is most likely responsible for the GAL binding observed in glioblastomas, idicating its influence on GBM differentiation and growth 13 . What's more it is also suggested that GALR1 mutations are responsible for the overexpression of GALR1 in grade IV tumors. GRM5 encodes the metabotropic glutamate receptor 5 that mediates post-synaptic NMDA receptor (NMDAR) currents 14 . It was reported that GRM5 was aberrantly expressed in brain gliomas 15 , suggesting its involvement in both brain and spinal glioma progression. The identification of genes associated with spinal gliomas provides new therapeutic targets and facilitates the development of biomakers for early tumor screening.
In our profiling of the 9 spinal glioma samples, 11 in-frame/frameshift fusions were identified. Notably, these identifications were increased with tumor progression from grade II to IV. This prevalence was consistent with the previously reported study in brain gliomas 16 , which means that the proportion of fusions in grades IV are higher than that in grade II. KIAA1549-BRAF fusions have been identified 17 and it was considered a causes to MAPK activation in pilocytic astrocytoma [17][18][19] . Besides, those previously reported studies, we have identified a number of novel fusions, of which GATSL2-GTF2I were detected in all of grade II samples, suggesting that the gene fusions could be related with the abnormal gene expression observed in spinal gliomas. Most of the identified fusion transcripts involved gene sequences have not been studied in GBMs. Therefore, characterizing the function of these fusions may unravel novel biological mechanisms of spinal glioma progression.
In conclusion, a complex landscape of molecular alterations in spinal gliomas across different tumor grades was revealed in this study. This advances our understanding of the progression of these tumors in the Chinese population. Further investigations of the network of these genes will further identify the characterization of their underlying mechanisms during the development of aggressive spinal gliomas.

Patients and samples. Tissues were collected from patients undergoing treatment at the Beijing Tiantan
Hospital and Beijing Tsinghua Changgung Hospital from 2012 to 2017. Informed consent was obtained from study participants according to institutional guidelines. Tissue samples were snap-frozen in liquid nitrogen. Tumor grades were diagnosed by two experienced pathologists. The clinical information of each patient is shown in Table 1.
All procedures involving human tumor specimens were performed in accordance with the ethical standards of our local research committee and the Helsinki declaration The study was approved by the Human Research Ethical Committee of Beijing Tsinghua Chenggung Hospital, China.
Histopathological diagnosis and mutation detection. Each histopathological slide of the patient sample was reviewed by at least one experienced neuropathologist. The H3K27M, BRAF (V600E), ATRX (all mutations that may affect protein function), TERT (C228T, C250T), IDH1 (R132), IDH2 (R172) mutation status were sequenced by Genome Analysis Toolkit (GATK) which can be used to identify H3K27M mutation type.
RNA-seq and quality control. RNA library construction and sequencing experiments were conducted at Anhui Anlongen Co., Ltd., (Hefei, China). Libraries were sequenced on an Illumina HiSeq 2000 platform using 150-bp paired-end sequencing. Generated images were converted into nucleotide sequences using a base-calling  Table 4. Fusion genes identified in grade II and grade IV spinal cord gliomas. Note: Left Break point: The genome position where the upstream gene break point is located. Right Break point: The genome position where the downstream genes break points is located. Verified By RT-PCR: "Yes" mean that the fusion genes were verified by RT-PCR sucessfully, and "No" mean that the fusion genes were not verified by RT-PCR successfully.