ANXA1 as a Prognostic and Immune Microenvironmental Marker for Gliomas Based on Transcriptomic Analysis and Experimental Validation

The tumor microenvironment (TME) plays an important role in the growth and invasion of glioma. This study aimed to analyze the composition of the immune microenvironment in glioma samples and analyze the important differentially expressed genes to identify novel immune-targeted therapy for glioma. We downloaded transcriptomic data of 669 glioma samples from The Cancer Genome Atlas database. CIBERSORT and ESTIMATE methods were used to calculate the proportion of tumor-infiltrating immune cells and ratio of immune and stromal components in the TME. The differentially expressed genes (DEGs) were screened by comparing the genes expressed by both stromal and immune cells. Annexin A1 (ANXA1) was determined to be an important prognostic indicator through the common overlap of univariate Cox regression analysis and protein–protein interaction network analysis. The proportion of tumor-infiltrating immune cells, calculated by CIBERSORT algorithm, had a significant difference in distribution among the high and low ANXA1 expression groups, indicating that ANXA1 could be an important immune marker of TME. Furthermore, ANXA1 level was positively correlated with the histopathological factors and negatively related to the survival of glioma patients based on the analysis of multiple databases. Finally, in vitro experiments verified that antagonizing ANXA1 expression promoted cell apoptosis and inhibited the invasion and migration capacities of glioma cells. Therefore, ANXA1 due to its immune-related functions, can be an important prognostic indicator and immune microenvironmental marker for gliomas. Further studies are warranted to confirm ANXA1 as a potential immunotherapeutic target for gliomas.


INTRODUCTION
Glioma is the most common primary intracranial tumor, characterized by high recurrence, easy invasion, and high mortality (Ostrom et al., 2014). Currently, the treatment for gliomas includes surgical resection, chemotherapy, and radiotherapy; however, the prognosis remains poor (Darefsky et al., 2012;Johnson et al., 2012;Koshy et al., 2012). An improved understanding of the pathogenesis and identification of the key molecular biomarkers may help in improving diagnostic accuracy and finding novel therapeutic targets for glioma to achieve better clinical outcomes.
A large body of evidence has shown the importance of tumor microenvironment (TME) in the development of tumors (Caponnetto et al., 2020;Jung et al., 2020;Yi et al., 2020). The microenvironment of glioma is composed of tumor cells, immune cells, stromal cells, and their various secreted factors. Among them, the factors secreted by tumor or immune cells, such as growth factors, chemokines, and pro-inflammatory, and anti-inflammatory factors, constitute a microenvironmental network that interact with each other to collectively regulate the process of tumorigenesis. Therefore, TME is the basis of tumor pathogenesis and is an important therapeutic target (Wood et al., 2014;Schulz et al., 2019). Previous studies have shown that the stromal cells in the TME are conducive to the tumor growth and invasion, and have an anti-tumor immune effect (Guo and Deng, 2018;Giraldo et al., 2019). Furthermore, immune cells, such as T cells and M1 and M2 macrophages, play a key role in the anti-tumor response (Bamodu et al., 2019;Heymann et al., 2019). Several immunotherapies have been shown to be effective, including immune checkpoint inhibition and therapeutic antibody applications. Some processes to serve as immunotherapy targets have been identified, such as programmed death-ligand 1 expression (Kumar et al., 2017;Koh et al., 2020), DNA mismatch repair deficiency (Overman et al., 2017;Baretti and Le, 2018), and tumor mutation burden (Chan et al., 2019;Jang et al., 2020). However, only a fraction of patients benefited from these strategies (Qi et al., 2020). Therefore, there is need to identify novel therapeutic targets to develop improved therapies with higher efficacy.
Therefore, how to accurately evaluate the dynamic composition of immune and stromal cells in TME of glioma is an important issue to address. With the development of genomic analysis, we can calculate the fraction of immune cells and stromal cells in each cancer sample and analyze the relationship between the two. We performed the ESTIMATE and CIBERSORT algorithms aiming to compute the proportion of tumor-infiltrating immune cells (TIC) and the ratio of immune and stromal components of glioma. In addition, we screened differentially expressed genes (DEGs) to identify novel prognostic indicators for immune-targeted therapy for glioma. We found annexin A1 (ANXA1) by evaluating the immune microenvironment of the glioma samples from The Cancer Genome Atlas (TCGA), and further verified the finding using Chinese Glioma Genome Atlas (CGGA) databases and in vitro experiments.

Cell Culture and siRNA Transfection
Human glioma cell lines, U87 MG and U251 were obtained from the Shanghai Institute of Biosciences and Cell Resources Center (Chinese Academy of Sciences, Shanghai, China). U87 MG and U251 cells were cultured in dulbecco's modified eagle medium (DMEM) with 10% fetal bovine serum (FBS) in 6-well plates. All cells were incubated in a humidified cell incubator with 5% CO 2 at 37 • C. The cells were cultured until 70-90% confluence before transfection. Transfection was performed using Lipofectamine 3,000 reagent as per the manufacturer's protocol. The cells were divided into siRNA group, negative control (NC) group and normal group. The siRNA group was transfected with ANXA1 siRNA and the NC group was transfected with NC siRNA for 24 h, the normal group did not receive any intervention.

Western Blot Analysis
Cancer cell lines were seeded in 6-well plates at a density of 5 × 10 4 cells per well and incubated over-night. Then siRNA group was transfected with ANXA1 siRNA and the NC group was transfected with NC siRNA for 24 h, the normal group did not receive any intervention. After replaced with fresh medium, the cells were continued to be incubated for 24 h. The cells were washed with phosphate-buffered saline (PBS) three times and harvested protein using RIPA lysis buffer with 1% PMSF. The protein concentration was determined using the BCA protein assay. The proteins were isolated using sodium dodecyl sulfate-polyacrylamide gel electrophoresis and transferred to a polyvinylidene fluoride membrane. After blocking it with 5% skimmed milk for 2 h and then the membrane was incubated with primary antibody (1:1,000) overnight. On the next day, the secondary antibody (1:5,000) was incubated for 2 h and the target proteins were visualized using enhanced chemiluminescence. Each experiment was done in triplicate and repeated three times independently.

Flow Cytometry
Cancer cell lines were seeded in 6-well plates at a density of 5 × 10 4 cells per well and incubated over-night. Then the siRNA group was transfected with ANXA1 siRNA and the NC group was transfected with NC siRNA for 24 h, the normal group did not receive any intervention. After replaced with fresh medium, the cells were continued to be incubated for 24 h. After washing the cells with PBS three times, they were digested using trypsin and collected in a 15 ml centrifuge tube. The resuspended cells were centrifuged at 1,000 g for 5 min. The supernatant was discarded, and cells were lightly resuspended in 195 µl Annexin V-FITC binding solution. Five microliters of annexin V-FITC was added and mixed gently. After incubating at room temperature in the dark for 20 min, flow cytometry analysis was performed to detect apoptosis on a FACSCalibur (BD Biosciences; Baltimore, MD, United States). Data were analyzed through the Flowjo software. Each experiment was done in triplicate for three independent experiments.

Wound Healing Assay, Invasion Assay, and Colony Formation Assay
Cell migration was detected and evaluated by wounding healing assays. Briefly, Cancer cell lines were seeded in 6-well plates at a density of 5 × 10 4 cells per well and incubated over-night. Then the siRNA group was transfected with ANXA1 siRNA and the NC group was transfected with NC siRNA for 24 h, the normal group did not receive any intervention. After replaced with fresh medium, the cells were continued to be incubated for 24 h. Scrapes were made with 200 µl sterile tips when the 6-well plates were 100% confluent. Finally, the scrape width was observed after 12 and 24 h under a microscope in different group. Invasion assay was performed to assess the ability of cell invasion. Cell  cultured in DMEM medium without 10% FBS overnight in the upper chamber in 8-µm chambers of 24-well plates containing solidified matrigel. About 800 µl of DMEM medium with 10% FBS was added to the lower chambers. The chamber was washed with PBS to remove the matrigel after 24 h. Finally, the cells were fixed with methanol for 15 min, washed with PBS three times, and the number of invasion cells were manually counted under a microscope after 0.1% crystal violet staining. Colony formation assay was performed to assess the proliferative capacity of the cells. Cell treatment is as described above. Cancer cell lines were digested by trypsin and suspended again after transfected with ANXA1 siRNA and NC siRNA. Then the cells (200 cells/well) were cultured in 6-well plates for 14 days. Finally, the cells were fixed using methanol for 15 min, washed three times with PBS, and counted using a microscope after 0.1% crystal violet staining. Colonies with more than 50 cells were considered for manual calculation. Each experiment was done in triplicate for three independent experiments.

Bioinformatics Analysis of TCGA and CGGA Database
We downloaded the mRNA expression profiles and clinical information of low-grade glioma (LGG) and glioblastoma (GBM) samples available on TCGA 1 dataset. mRNA expression profiles were downloaded as Fragments Per Kilobase of transcript per Million mapped reads (FPKM). Molecular subtype and treatment details were downloaded from previous study (Ceccarelli et al., 2016). Then, the "sva" package was utilized for the normalization of RNA expression profiles and to remove the batch effects between TCGA-LGG samples and TCGA-GBM samples (Leek et al., 2012). Furthermore, we downloaded 749 glioma samples and clinical information from the CGGA, 2 which were used as a validation dataset. mRNA expression profiles were downloaded as Expectation Maximization (RSEM). The "sva" package was performed to remove the batch effects between different glioma samples.

Calculation of Immune and Stromal
The immune and stromal scores of each glioma sample from TCGA were calculated using "limma" and "estimate" packages in R language 3.6.1, to evaluate the proportion of immune and stromal components in the TME. The higher the score, the larger the proportion of the corresponding component in the TME. According to the median score achieved, the samples were divided into high and low score groups.

Identification of DEGs and Survival Analysis
Differential expression analysis was performed to explore the differences in gene expression between the high and low score groups. Differentially expressed genes (DEGs) were analyzed using "limma" package in R and screened by the comparison between the high-score group and the low-score group. DEGs were visualized using "pheatmap" package in R. DEGs with |log 2 fold change| > 2 and false discovery rate (FDR) < 0.01 were considered significant. Survival analysis was performed using "survival" and "survminer" packages, to verify the difference in survival rate between the high and low score groups. The survival curve was plotted using the Kaplan-Meier method, and log rank test was used assess the statistical significance. A p-value < 0.05 was considered statistically significant.

Protein-Protein Interaction Network Construction and Cox Regression Analysis
Protein-protein interaction (PPI) network of DEGs was constructed using the STRING database 3 and Cytoscape software (Otasek et al., 2019). Nodes with confidence of interactive relationship greater than 0.95 were screened for network construction. Univariate Cox regression analysis was performed using the "survival" package to screen the genes that were closely associated with the prognosis.

Differential Expression Analysis Between TIC and ANXA1 Expression
The TIC profile of gliomas was estimated using the CIBERSORT computational method available in R. Subsequently, the tumor Frontiers in Cell and Developmental Biology | www.frontiersin.org 6 August 2021 | Volume 9 | Article 659080 samples were screened for quality and only those with p < 0.05 were selected for subsequent analysis. The visualization of TIC profile was constructed using the "vioplot" and "barplot" package. The differential expression analysis between ANXA1 expression and immune cells was performed using "limma" package of R language.

Clinical Outcomes and Clinicopathological Features With ANXA1 Expression
Survival analysis of clinical outcomes was performed using "survival" and "survminer" packages, to verify the difference in survival rate between the high and low ANXA1 expression groups. The survival curve was plotted using the Kaplan-Meier method, and log rank test was used assess the statistical significance. A p-value < 0.05 was considered statistically significant. The relationship of clinicopathological features and ANXA1 expression were visualized using "pheatmap" package in R.
database. The overall survival and disease-free survival were also analyzed using GEPIA. The protein expression level of ANXA1 analysis was performed using Human Protein Atlas (HPA). 6

Statistical Analysis
Statistical analysis and charts preparation were carried out using R (version 3.6.3) and SPSS (version 25). Student's t-test was used for comparison between the two groups, and variance analysis was used for data comparison between multiple groups. p-values < 0.05 were considered statistically significant.

Data Filtering and Processing
The study design is showed in Figure 1.

Immune and Stromal Scores Are Associated With Prognosis of Glioma Patients
To evaluate the relationship between the Immune scores and Stromal scores and survival prognosis, survival analysis was performed for the Immune and Stromal scores. The ESTIMATE algorithms were used to calculate the proportion of the number of immune and stromal components in glioma samples from TCGA. According to the median score of each sample, the patients were divided into high and low score groups. The Immune (Figure 2A) and Stromal ( Figure 2B) scores were found to be negatively associated with the overall survival rate (p < 0.001). This suggests that the immune and stromal components of the TME can be used to predict the survival of glioma patients.

Identification of DEGs Depended on the Immune and Stromal Scores
The DEGs were determined by differential expression analysis of immune cells and stromal cells in the TME compared between high-and low-score samples and the 100 DEGs with the most significant p-values was visualized using heatmap analysis ( Figure 2C). A total of 679 and 4,689 DEG were screened from the Immune score group and Stromal score group, respectively. The Venn diagram analysis showed that 532 upregulated genes and 131 downregulated genes overlapped in the Immune score and Stromal score groups ( Figure 2D).

Intersectional Analysis of PPI Network and Univariate Cox Regression
To further indicate its potential mechanism, we used Cytoscape software to construct a PPI network based on the STRING database. The PPI network of DEGs is shown in Figure 3A, and the histogram of the first 30 genes sort by the node numbers is displayed in Figure 3B. Univariate Cox regression analysis is performed to identify the top 30 significant genes ranked by the p-value among of the 633 DEG in glioma samples ( Figure 3C). A Venn diagram between the first 30 leading nodes of the PPI network and the top 30 genes obtained from the univariate Cox regression analysis, were displayed, and only ANXA1 was found in the overlapped region from the above analyses ( Figure 3D). These results showed that the ANXA1 cloud be a prognostic indicator of glioma.

ANXA1 Is Associated With the Proportion of TIC
To further verify the association of ANXA1 expression with the TME, the proportion of TIC was analyzed using CIBERSORT computing methods, and 22 types of immune-cell distribution levels in glioma samples were built (Figure 4). Immune cells of T cells CD8, Monocytes, Macrophages M0 and Macrophages M1 have significant differences between different ANXA1 expression groups in TCGA LGG samples and CGGA LGG samples ( Figures 4A,C), while immune cells of dendritic cells resting has significant differences in TCGA GBM samples and CGGA GBM samples (Figures 4B,D). The kinds of TICs correlated with ANXA1 expression in TCGA and CGGA database were described in the Table 2. These results further indicated that ANXA1 was involved in the immune activity of TME and the maintenance of cell-mediated immunity in the LGG.

Different ANXA1 Expression Profiles Had Distinct Clinical Outcomes and Clinicopathological Features
To verify the relationship between ANXA1 expression and survival prognosis, all glioma samples were divided into high ANXA1 expression and low ANXA1 expression groups based on the median ANXA1 expression. Glioma patients with higher ANXA1 expression, including LGG and GBM samples, had a significantly shorter survival rate than those with a lower expression from the TCGA database (P < 0.05) (Figures 5A,B); The survival analysis results of CGGA were consistent with those of TCGA (P < 0.01) (Figures 5C,D). Further we compared the clinicopathological features between the high and low ANXA1 expression group. The high ANXA1 expression group was significantly correlated with higher grade (P < 0.001), IDH wild status (P < 0.001) and 1p19q codeletion (P < 0.001) in LGG samples ( Figure 5E); In the GBM samples, the high ANXA1 expression group is significantly correlated with older FIGURE 9 | Wound healing assay, Transwell assay, and colony formation assay. (A) The migration ability of siANXA1 group was significantly reduced compared with NC and normal group (100X). (B) The cell cloning ability of siANXA1 group was also significantly reduced compared with NC and normal group. (C) The invasion ability of siANXA1 group cells was reduced compared with that of NC and normal group (100X). **p < 0.01; ***p < 0.001.
(P < 0.05) and IDH wild status (P < 0.001) from TCGA database ( Figure 5F). In the CGGA database, the high ANXA1 expression group also was significantly correlated with IDH wild status (P < 0.001) and 1p19q codeletion (P < 0.001) in LGG samples ( Figure 5G); Meanwhile, it also was significantly correlated with older (P < 0.01) and IDH wild status (P < 0.001) in GBM samples ( Figure 5H). These findings clearly indicate that ANXA1 expression in TME is negatively associated with the prognosis of glioma. In particular, the expression of ANXA1 was found to increase with the progression of tumor grade and malignancy.

Mutation and Prognostic Analysis of ANXA1 Expression in Multiple Databases
To further investigate the role of ANXA1 expression in glioma, we performed differential expression and mutation analyses of ANXA1 using GEPIA (gepia.cancer-pku.cn) and cBioPortal (see text footnote 4) databases. According to the GEPIA database, the results of differential expression analysis revealed an increased expression of ANXA1 in the tumor samples compared with normal samples in LGG samples ( Figure 6A) and GBM samples ( Figure 6B). The overall survival and disease-free survival rate of patients with low ANXA1 expression were significantly better than those with high expression in LGG samples (P < 0.001) ( Figure 6C) and GBM samples ( Figure 6D). The purpose of mutation analysis is to better identify the mutation status and sites of ANXA1 in glioma, and mutation sites are common targets of gene therapy. The mutation type of ANXA1 in GBM was mainly amplification, while in low-grade gliomas was mainly mutation (Figure 6E). The mutation analysis of ANXA1 showed that the mutation frequency of ANXA1 were 3 and 7% in LGG and GBM samples, respectively ( Figure 6F). Mutation frequency of ANXA1 was higher in GBM samples than in the LGG samples and the common mutation sites were Q23 and R303C ( Figure 6G). These results further validated that ANXA1 was closely related to prognosis and could serve as a therapeutic target in glioma.

Correlation of ANXA1 Expression With Immune Checkpoint Markers of Glioma
In the CGGA dataset, we analyzed relation between the ANXA1 expression and immune checkpoint markers (CD274, CD276, LAG3, CD86, TIMP1, and CHI3L1). We found that the ANXA1 expression was positively correlated with these markers (Figure 7A) and strongly related with CD274 and CD276 expression (Figures 7B,C). Then we download the low-grade and high-grade glioma samples of the same patients from HPA database. The result showed that ANXA1 and CD274 expression (patient id:122; 3,226) were positive in high grade gliomas, and there is a certain correlation between their expressions (Figures 7D,E). ANXA1 and CD276 (Patient id: 3,174;3,241) also were overexpressed in high grade gliomas (Figures 7F,G). These results indicated ANXA1 influenced the tumor microenvironment by regulating these tumor immune genes.

Silencing ANXA1 Expression Promotes Cell Apoptosis by Targeting BCL-2 and Bax
The TME can influence the cell proliferation and death of tumor cells, and apoptosis is a common form of cell death (Su et al., 2015;Yaacoub et al., 2016). Therefore, we further studied the relationship between ANXA1 and cell apoptosis. In the CGGA dataset, we analyzed relation between the ANXA1 expression and apoptosis markers (BAX, BCL-2, CASP3, CYCS, TP53, and PAPR). We found that the ANXA1 expression was positively correlated with these markers ( Figure 8A) and strongly related with BAX expression (Figure 8B). To investigate the apoptotic pathways involved in cell death due to inhibition of ANXA1 expression, we performed western blot to analyze the change in expression levels of the protein involved in cell apoptosis caused by silencing ANXA1 expression in U87 and U251 cell lines. The results showed that in cells transfected with siANXA1, bcl-2 expression was downregulated and Bax was upregulated, which regulated the apoptosis (P < 0.05) (Figures 8C,D). Then the flow cytometry assays were performed to assess the cell apoptosis. the flow cytometry analysis also revealed a significantly higher rate of apoptosis in the siANXA1 group (P < 0.01) (Figures 8E,F). These results verified the hypothesis that inhibition of ANXA1 expression induced cell apoptosis in U87 and U251 cells.

Silencing ANXA1 Reduced Migration, Invasion, and Proliferation in Glioma Cell Lines
Cell migration and invasion are important steps in tumor growth and therefore we evaluated the effects of siANXA1 on migration and invasion in U87 and 251 cell lines. Wound healing assay showed that the migration ability of the siANXA1 group was significantly reduced compared to that of the normal and NC group ( Figure 9A). Furthermore, the colony formation assay indicated that ability to form colonies was significantly reduced in the siANXA1 group compared to that in the normal and NC group ( Figure 9B). Finally, the Transwell assay showed that the siANXA1 group was less invasive than normal and NC group (P < 0.05) ( Figure 9C). These experimental results suggested that the inhibition of ANXA1 expression could reduce the metastatic potential of gliomas, which further verified that ANXA1 expression was related to malignancy in gliomas.

DISCUSSION
Despite the combination of surgery, radiotherapy, and chemotherapy, patients with glioma show poor clinical prognosis. Therefore, we need to find more effective treatment methods to improve survival in these patients. Many studies have shown that the TME plays an important role in the growth and invasion of glioma (Chang et al., 2016;Qian et al., 2018;Meng et al., 2019). Previous studies have shown that therapies targeting the tumor microenvironment have achieved good results, such as anti-programmed death-ligand 1 therapy (Ruan et al., 2019;Ene et al., 2020) and T cell immunotherapy (Brown et al., 2016;Choi et al., 2019). Moreover, immunotherapy has been shown to be effective in other cancers such as lung carcinoma and melanoma (Topalian et al., 2012). Upon evaluating the glioma samples in the TCGA database, our results showed that the immune components in TME contributed to the prognosis of glioma. Additionally, we found that besides the proportion of immune and stromal components of TME, tumor grade and histopathologic characteristics were significantly related to the malignancy of glioma. These results highlight the importance of exploring the interaction between tumor cells and immune cells and researching new therapeutic targets for the development of more effective treatments. The discovery of lymphatic system in the central nervous system has provided a new theoretical basis and opened up novel prospects for immunotherapy for brain tumors (Louveau et al., 2015). Multiple immune-related signaling pathways, such as the phosphatidyl inositol-3-kinase/Akt pathway, IGF pathway, and programmed death-ligand 1 signaling, individually or collectively impact the TME of glioma (Kravchenko et al., 2015;Wang et al., 2016;Cai et al., 2018). Based on transcriptome sequencing analysis of the glioma samples from TCGA database, we found that increased ANXA1 expression was significantly associated with advanced clinicopathological features and poor prognosis in patients with glioma. These results indicated that ANXA1 could be a prognostic indicator and an immunotherapy marker for TME in glioma.
ANXA1 is a membrane-localized protein that binds phospholipids and has an anti-inflammatory activity. It plays important roles in the innate immune response as an effector of glucocorticoid-mediated responses and a regulator of the inflammatory processes (Arcone et al., 1993;Leoni et al., 2015). It contributes to the adaptive immune response by enhancing the signaling cascades that are triggered by T-cell activation, regulates differentiation and proliferation of activated T-cells, promotes the differentiation of T-cells into Th1 cells, and negatively regulates their differentiation into Th2 cells (D'Acquisto et al., 2007). ANXA1 has been shown to be abnormally expressed in the immune microenvironment of cancer (Moraes et al., 2017(Moraes et al., , 2018, and antagonizing ANXA1 expression has often achieved good results (Hatakeyama et al., 2011;Gastardelo et al., 2014). Therefore, we further analyzed the relationship between ANXA1 expression and TME. The CIBERSORT computing method for the evaluating the proportion of TIC revealed a positive correlation between T-cell activation and ANXA1 expression in LGG glioma samples, suggesting that ANXA1 could be involved in the maintenance of cell-mediated immunity in the TME.
ANXA1 is differentially expressed in a variety of tumors and a number of genome-sequencing studies have shown that ANXA1 is highly expressed in gliomas (Ruano et al., 2008;Naryzhnyi et al., 2014;Mallawaaratchy et al., 2015;Qiu et al., 2020), but intensive research on ANXA1 function in glioma was lacking. To explore the underlying mechanism of the ANXA1 in TME, we investigated the relationship between ANXA1 and several important immune checkpoint markers through the HPA database. The result showed that NXA1 expression was positively correlated with CD274 and CD276 expression, which indicated ANXA1 could influenced the tumor microenvironment by regulating these tumor immune genes. Then we investigated the role of ANXA1 in cell proliferation and growth through in vitro experiments. The in vitro experiments showed that suppressing ANXA1 expression could significantly inhibit the proliferation and motility of glioma cells and promote cell death via apoptosis. Overall, these findings have proved that ANXA1, a key gene in the TME, can be an effective target for treatment of glioma.
These results suggest that targeted tumor microenvironment therapy for glioma can achieve satisfactory results. ANXA1 inhibition needs to be explored further in vitro and in vivo to better understand the efficacy of targeting this protein in suppressing growth and tumor progression of glioma. In future, continued screening of core genes in the glioma immune microenvironment can help us develop novel targets more accurately and provide a theoretical basis for the pathological mechanism of glioma.

DATA AVAILABILITY STATEMENT
The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found in the article/supplementary material.

AUTHOR CONTRIBUTIONS
ZL, EY, XL, and HW designed the study and wrote the initial draft of the manuscript. MW, JC, CY, HZ, JR, KW, YZ, and LH contributed to data analysis. QZ and SY reviewed and edited the manuscript. All authors read and approved the manuscript.