All RNA-sequencing data of diffuse glioma patients were obtained from two independent databases: the CGGA dataset (Dataset 1, n = 325 cases) (http://www.cgga.org.cn) and the CGGA dataset (Dataset 2, n = 693 cases) (Zhao et al., 2021). To compare the gene expression patterns of tumor tissues and normal brain tissues, we also collected 20 RNA-seq samples of normal brain tissues from the CGGA database in this study. All WES data of diffuse glioma patients from WHO II-IV were obtained from the CGGA Network (Dataset 3, n = 286). Clinical information of all patients was also collected from the CGGA Network, such as WHO grade (WHO II-IV), histology grade (oligodendroglioma, anaplastic oligodendroglioma, astrocytoma, anaplastic astrocytoma, and glioblastoma, abbreviated as O, AO, A, AA, and GBMs, respectively); gender, age, and overall survival data; progression status (primary and recurrent states); and molecular pathological features (IDH mutation status, MGMT promoter methylation status, and chromosome 1p/19q co-deletion status). This research was approved by the Ethics Committee of Capital Medical University, and all patients provided written informed consent.

To further explore ANXA1 expression in different histological regions of GBM blocks, we obtained Ivy data from the Ivy Glioblastoma Atlas Project–Allen Institute for Brain Science datasets (Dataset 4, n = 270 cases) (Puchalski et al., 2018) (http://glioblastoma.alleninstitute.org/). For this dataset, we collected different histological regions that contain 1) cellular tumor (CT), 2) infiltrating tumor (IT), 3) leading edge (LE), 4) microvascular proliferation (MP), and 5) pseudopalisading cells (PC).

The scRNA-seq data of diffuse glioma patients were obtained from a previous study (Neftel et al., 2019) (https://singlecell.broadinstitute.org/single_cell/study/SCP393/single-cell-rna-seq-of-adult-and-pediatric-glioblastoma). Of them, there are 6863 tumor cells, 754 macrophages, 219 oligodendrocytes, and 94 T cells (Dataset 5). For tumor cells, we also obtained four types of cellular state annotations that recapitulate 1) neural-progenitor–like (NPC-like), 2) oligodendrocyte-progenitor–like (OPC-like), 3) astrocyte-like (AC-like), and 4) mesenchymal-like (MES-like) states.

WES data were mapped to the human reference genome (hg19) using the Burrows–Wheeler Aligner (BWA) tool (Li and Durbin, 2009) with default parameters. Then, SAMtools (Li et al., 2009) and Picard (http://broadinstitute.github.io/picard/) were used to sort the reads by coordinates and mark duplicates. Next, we used the CNVkit software (Talevich et al., 2016) to estimate the CNA status of well-known driver genes in gliomas, such as PTEN, MET, EGFR, and CDKN2A/B. In this study, a copy number gain is identified as log2 (ratio) larger than 0.5, while a copy number loss is identified as log2 (ratio) less than 1.0.

For RNA-seq data, we estimated the abundance of member cell types using the CIBERSORT method (Newman et al., 2015). We uploaded gene expression profiles and ran CIBERSORT software online (https://cibersort.stanford.edu/runcibersort.php) by selecting LM22 (gene signature) and 1000 permutation parameters. As result, we obtained the 22 kinds of cell composition for each sample from gene expression profiles.

For RNA-seq data, we identified the TCGA subtypes for each sample as previously described (Wang et al., 2017). In this pipeline, ssGSEA was performed to obtain the scores of the four signatures for each sample from gene expression profiles. Since the scores of the four signatures were not directly comparable, this pipeline was used to perform a resampling procedure to generate null distributions for each of the four subtypes (1000 permutations). Following this procedure, this method generated random ssGSEA scores for each subtype to provide empirical p-values and scaled ssGSEA scores for the raw ssGSEA scores of each sample. Finally, we assigned the TCGA subtypes for each sample based on the p-values and scaled ssGSEA scores.

The selected glioma samples were collected from the CGGA tissue bank and were supervised by the Beijing Tiantan Hospital Institutional Review Board (KY 2019-143-02). IHC analysis was performed as previously reported (Hu et al., 2018). Briefly, the slides were deparaffinized and boiled in antigen-retrieval buffer. Then, the slides were blocked using endogenous peroxidase with H₂O₂, subsequently blocking non-special sites, and the slides were incubated with primary antibodies against ANXA1 (Cell Signaling Technology #32934, 1:400 dilution) overnight at 4°C. On the second day, the slides were rinsed three times in PBS buffer and incubated with the secondary antibody working solution (PV6000 Beijing Zhongshan Jinqiao Biological Company) for 60 min at room temperature. Last, the IHC images were captured using an Axio Imager 2 microscope (Zeiss). The scores were calculated according to the intensity score multiplied by the areas as follows: The intensity was defined as follows: 0 for no staining, one for weak staining, two for moderate staining, and three for strong staining. The area score was determined as follows: 0 for less than 5% cells positive, 1 for 5–25% cells positive, 2 for 26–50% cells positive, 3 for 51–75% cells positive, and 4 for greater than 75% cells positive.

To investigate the biological functions of the ANXA1 gene, the ANXA1 coexpressed genes were obtained and gene set enrichment analysis (GSEA) (Subramanian et al., 2005) was performed. First, we downloaded the gene sets from the GSEA website (http://www.gsea-msigdb.org/gsea), including the Gene Ontology (GO) biological process, Kyoto Encyclopedia of Genes and Genomes (KEGG), and cancer hallmark. Then, ANXA1 coexpressed genes were obtained by the Pearson expression correlation analysis between ANXA1 and other genes. Finally, we implemented the ClusterProfiler R package to reach this process (Yu et al., 2012).

The R statistical software (v4.0.3) (http://www.r-project.org) was used for statistical calculations and graph rendering. The prognostic value of ANXA1 was estimated by using the Kaplan–Meier analysis and Cox proportional hazard model analysis using the “survival” and “survminer” packages in R. In this study, the Pearson correlation analysis was used to obtain ANXA1 coexpressed genes. In particular, a positive correlation is defined as a correlation coefficient larger than 0.6 and p-value < 0.05, while a negative correlation is defined as a correlation coefficient less than −0.6 and p-value < 0.05. The Wilcoxon test and one-way ANOVA test were used for two and multiple group comparisons, respectively. For all statistical methods, p < 0.05 was considered as a significant difference.

In this study, a total of 1,018 patients with gliomas aged 8–79 years (median ± sd, 42 ± 12 years) were included. The majority of glioma patients were males (59%) and WHO IV (38%), and there were 651 cases of primary gliomas. For these patients, 617 case deaths were recorded, with the median survival of NA (3470-NA), 1208 (1028-1657), and 378 (344-415) for WHO II, WHO III, and WHO IV, respectively. All patients with transcriptomic data were used to analyze ANXA1 expression, and 231 of patients were also performed with WES to investigate genetic changes. The clinical and pathological features of these patients are described in Table 1.

To explore ANXA1’s role in gliomas, we examined its transcriptomic level in different subtypes of gliomas in two batches of RNA-seq data from the CGGA database. We found that the expression values of ANXA1 were significantly higher in GBM patients than in those with normal brain and lower-grade gliomas (WHO II and WHO III) in Dataset 1 (p < 5e-5, Figure 1A). Our further results showed that the ANXA1 expression levels were statistically more abundant in GBMs than in other histology (p < 1e-2, Figure 1B). In addition, due to the genetic and clinical differences between IDH-mutated gliomas and IDH wild-type gliomas, we explored the role ANXA1 played in gliomas with different IDH statuses. The ANXA1 expression was highest in IDH wild-type and lowest in IDH mutation and 1p/19q co-deletion in LGGs (all p ≤ 5e-5, Figure 1C left), while ANXA1 expression was higher in the IDH wild-type than in IDH mutant gliomas in GBMs (p < 5e-11, Figure 1C right). There was a reduced expression of ANXA1 in glioma with IDH mutation based on LGGs and GBMs (all p < 5e-9, Figure 1D). It is well known that the MGMT promoter methylation status is a key biomarker indicating temozolomide (TMZ) chemotherapy sensitivity in gliomas. As a result, we found that patients without MGMT promoter methylation possessed a higher ANXA1 expression level in GBMs, suggesting that ANXA1 may play a resistance role in TMZ therapy of GBMs (p < 5e-3, Figure 1E). Notably, we also found that ANXA1 expression was higher in recurrent LGGs (Figure 1F). The aforementioned results are well validated in independent CGGA RNA-seq data (Figures 1G–L). Consistently, the immunohistochemistry (IHC) experiments of glioma patients (WHO II–IV grade) showed that ANXA1 was the highest in WHO IV patients and lowest in WHO II patients (all p < 0.05, Figures 1M–N). Taken together, these results suggest that the ANXA1 gene acts as an oncogene and may serve as a biomarker for disease progression in gliomas.

To investigate the clinical value of ANXA1, we examined the association between gene expression of ANXA1 and clinical information, including primary/recurrent status, WHO grade, histology, age, gender, well-known molecular status, TCGA subtype, survival, and therapy information. The evaluation of the association between clinicopathological features and the ANXA1 gene was conducted for 1,018 glioma patients from a Chinese cohort. Gliomas in Dataset 1 were ordered by increasing ANXA1 expression (Figure 2A). Our results showed that primary gliomas and LGGs had lower levels of ANXA1 expression (all p < 0.01), suggesting that ANXA1 may play a positive role in malignant progression. Younger patients with glioma had lower expression of ANXA1 (p < 2.14e-10). Gender of patients is not associated with ANXA1 expression. With regard to genomic alterations, IDH mutation, 1p/19q co-deletion, and MGMT promoter methylation indicated lower ANXA1 expression (all p < 0.01). Gliomas with lower ANXA1 expression are more likely to belong to proneural (PN) and neural (NE) subtypes and have a good prognosis, while gliomas with high ANXA1 expression are more likely to belong to mesenchymal (MES) and classical (CL) subtypes and have poor survival (all p < 2.00e-16). Gliomas with chemotherapy and/or radiotherapy tend to have a high ANXA1 expression. The aforementioned results are well validated in independent CGGA RNA-seq data (Dataset 2, Figure 2B). These results indicate that ANXA1 is closely related to clinical behavior.

To investigate the association between ANXA1 expression and genomic alterations, we analyzed the somatic mutations and copy number alteration data from cases with RNA-seq and WES data for this purpose. In total, 231 samples in the entire cohort harbored both RNA-seq and WES data (Dataset 3). Recapitulating previous studies, we confirmed frequency mutation in IDH, TP53, ATRX, CIC, NOTCH1, EGFR, and PDGFRA in this study. According to ANXA1 expression, gliomas were divided into G1 group (low expression, n = 116) and G2 group (high expression, n = 115). Approximately two-thirds of cases in the G1 group carried either an IDH1 mutation or IDH2 mutation.

In addition, cases in the G1 group were enriched in CIC and NOTCH1 mutation that have been well-described in oligodendroglia histology (Figure 3). In contrast, both TP53 and ATRX mutation in cases of the G2 group were 1.25 times higher than those in the G1 group. On the other hand, cases in the G2 group have a much higher mutation frequency of EGFR and PDGFRA than those in the G1 group. Notably, although previously not recognized, mutations in RYR2, IGSF10, BNC2, CADPS2, COL12A1, TRABD2A, and USP34 were significantly enriched in the G2 group. Moreover, we also explored the frequency of copy number alterations in G1 and G2 groups. For CN amplification, G1 had a higher alteration frequency in AHNAK and CD276, while G2 had a high alteration frequency in EGFR, PDGFRA, MET, and TTN. For CN loss, our results showed that deletion in CDKNA2A/B genes in interferon-α family and olfactory receptor family 4 subfamilies mainly occurred in G2 cases. Taken together, glioma subtypes classified by ANXA1 expression showed distinct mutation and CNA features.

ANXA1 expression was heterogeneous in different glioma subtypes. To explore ANXA1’s biological role in gliomas, RNA-seq data were collected. First, we obtained the genes that significantly correlated with ANXA1 expression (Pearson |R| > 0.6 and p < 0.05). Totally, 462 positive and 107 negative coexpressed genes were identified in Dataset 1. Then, we predicted the GO biological process and cancer hallmark of these coexpressed genes. GSEA showed that the coexpressed positive genes associated with ANXA1 were mainly involved in immune-related functions, such as interferon-gamma response and regulation of innate immune response, suggesting a regulatory role in the immune microenvironment in gliomas (Figure 4A). In particular, we found that these coexpressed genes also positively function in apoptosis, epithelial–mesenchymal transition, NF-κB signaling, etc., indicating that ANXA1 may play an important role in regulating cell fate in gliomas. In contrast, we found that coexpressed negative genes of ANXA1 participate in the neuro-basic functions in gliomas, such as synapse structure and organization, regulation of cellular component biogenesis, and neuro-projection morphogenesis and differentiation. Furthermore, GSEA verified that ANXA1 was associated with immune, apoptosis, and neuron function (Figure 4B). For genes in interferon-gamma response, we confirmed that they are associated with ANXA1 expression, and show the differential expressed patterns in glioma subtypes grouped by ANXA1 expression (Figure 4C). In summary, ANXA1 correlates with immune-related function and cancer hallmark and plays a critical role in the glioma ecosystem.

As ANXA1 confers an extended immune status, we sought to further explore the ANXA1 regulatory immune role in the glioma ecosystem. We applied CIBERSORT software on Dataset 1 for estimating the relative abundances of 22 infiltrating immune cells (Newman et al., 2015). These cells mainly include lymphocytes, plasma, myeloid cells, and eosinophils. As a result, the majority of cell types in gliomas are myeloid cells and lymphocytes. In addition, we found that M2 macrophages are significantly enriched in gliomas with high ANXA1 expression (Figure 5A). We also validated that ANXA1 expression exhibited a significant positive correlation with the expression of M2-related genes (all R > 0.6 and p < 2.2e-16), including CD276, CLE7A, CTSA, FN1, IL4R, MMP9, MSR1, TGFB1, and VEGFA (Figure 5B), suggesting that ANXA1 acts a potential regulatory factor for M2 macrophages. In addition, we further collected anatomic transcriptional data in gliomas (Dataset 4), including leading edge (LE), infiltrating tumor (IT), cellular tumor (CT), pseudopalisading cells around necrosis (PAN), and microvascular proliferation (MVP) (Puchalski et al., 2018). Therefore, we found that ANXA1 was significantly under-expressed in CL enriched in the PN TCGA subtype and significantly overexpressed in MVP enriched in the MES TCGA subtype. This result is consistent with previous findings that ANXA1 was highly expressed in MES gliomas. To further explore ANXA1’s role in the tumor environment, we collected single-cell transcriptomic data in gliomas from a previous study (Neftel et al., 2019) (Dataset 5). We found that ANXA1 is highly expressed in macrophages, indicating a potential role for macrophages, especially M2 macrophages (Figure 5D). Moreover, we also noticed that tumor cells are highly expressed in the ANXA1 gene. In the single-cell level, our result showed that tumor cells with high expression of ANXA1 are in the MES cellular state (Figure 5E), indicating that ANXA1 could drive transitions to MES-like states in gliomas as reported in a previous study (Hara et al., 2021). In summary, we found that ANXA1 is highly expressed in M2 macrophages and MES tumor cells.

To further explore the role of the ANXA1 gene in clinical application, we examined the prognostic value in all kinds of subtypes in gliomas. We used the quartile of ANXA1 expression to divide the samples into three groups and explore their prognostic differences (Dataset 1). Gliomas with high expression levels of ANXA1 showed a significant poor prognosis for overall survival (OS) in both gliomas and LGGs (log-rank test, p < 1.0e-4, Figures 6A,B). We also found that ANXA1 expression stratified patients with MGMT promoter methylation into distinct survival groups (log-rank test, p < 1.0e-4, Figure 6C), assuming that patients previously thought to be sensitive to TMZ could be stratified based on ANXA1 expression. Meanwhile, our results suggest that patients previously thought to be resistant to TMZ can be stratified based on ANXA1 expression, and patients with low ANXA1 expression could also have a good prognosis (log-rank test, p < 1.0e-4, Figure 6D). These analyses in Dataset 2 were conducted in parallel (log-rank test, all p < 1.0e-4, Figures 6E–H). Furthermore, we conducted the univariate and multivariable Cox regression analyses in Dataset 1, which implies that ANXA1 expression is an independent predictor for survival prognosis after adjusting for other clinicopathological factors (Table 2). These results suggest that ANXA1 expression is significantly correlated with patient outcome.