Integrative Analysis of DNA Methylation and Transcriptome Identifies a Predictive Epigenetic Signature Associated With Immune Infiltration in Gliomas

Glioma is the most common primary brain tumor with poor prognosis and high mortality. The purpose of this study was to use the epigenetic signature to predict prognosis and evaluate the degree of immune infiltration in gliomas. We integrated gene expression profiles and DNA methylation data of lower-grade glioma and glioblastoma to explore epigenetic differences and associated differences in biological function. Cox regression and lasso analysis were used to develop an epigenetic signature based on eight DNA methylation sites to predict prognosis of glioma patients. Kaplan–Meier analysis showed that the overall survival time of high- and low-risk groups was significantly separated, and ROC analysis verified that the model had great predictive ability. In addition, we constructed a nomogram based on age, sex, 1p/19q status, glioma type, and risk score. The epigenetic signature was obviously associated with tumor purity, immune checkpoints, and tumor-immune infiltrating cells (CD8+ T cells, gamma delta T cells, M0 macrophages, M1 macrophages, M2 macrophages, activated NK cells, monocytes, and activated mast cells) and thus, it may find application as a guide for the evaluation of immune infiltration or in treatment decisions in immunotherapy.


INTRODUCTION
Gliomas are the most common primary tumors of the central nervous system, accounting for 60% of craniocerebral tumors (Ostrom et al., 2019). According to the malignant features, the World Health Organization (WHO) divides gliomas into four grades. With the continuous discovery of new molecular markers, the accurate classification of gliomas has been promoted. Currently, the WHO recommends that molecular markers such as co-deletion of chromosome arms 1p and 19q (1p/19q co-deletion) and isocitrate dehydrogenase (IDH) mutation status be included in the histopathological classification of gliomas (Louis et al., 2016). Surgery combined with chemotherapy or radiotherapy is a common treatment strategy for gliomas, but the prognosis varies greatly (van den Bent et al., 2013;van den Bent, 2014;Buckner et al., 2016). In particular, glioblastoma (GBM), or grade 4 glioma, is the most common and fatal malignant brain tumor, with a short median survival of only 12 to 15 months (Hanif et al., 2017). Therefore, the development of effective biomarkers for risk assessment may help to guide the choice of future targeted treatment strategies.
The complex tumor microenvironment (TME) has a serious influence on the gene expression in cancer tissue, affecting disease progression and clinical outcome (Cooper et al., 2012;Yoshihara et al., 2013). In the progression of gliomas, due to the destruction of the blood-brain barrier, tumor cells secrete a large number of chemokines that recruit immune cells from the peripheral blood, and the proportion of microglia/macrophages in the tumor cavity increases, which is associated with the malignant degree of gliomas (Yi et al., 2011). Therefore, the evaluation of immune infiltration of gliomas plays an important role in monitoring the progression of gliomas, and tumor-infiltrating immune cells have the potential to be used as drug targets to improve the prognosis of glioma patients (Xiong et al., 2018). Immunotherapy suppresses the rejection of tumor cells by the host immune system by inducing, enhancing, or suppressing the immune response, and thus uses the host immune system to kill tumor cells. At present, immune checkpoint inhibitors, such as anti-CTLA-4 monoclonal antibody (ipilimumab) and anti-PD-1 monoclonal antibody (nivolumab), have achieved good results in tumor therapy (O'Day et al., 2010;Rizvi et al., 2015). However, research on immune checkpoint inhibitors in glioma is still insufficient, and new drugs are currently being evaluated in animal models and in clinical trials (Omuro et al., 2018;Crommentuijn et al., 2020;Ladomersky et al., 2020).
DNA methylation plays an important role in maintaining the structural stability of the genome and in regulating gene expression. Epigenetic modification of DNA, including abnormal DNA methylation, especially in the promoter region of genes as such as CDKN2A, has been reported to be associated with tumorigenesis and has great potential as a tumor biomarkers (Malzkorn et al., 2011;Klutstein et al., 2016;Pan et al., 2018) with numerous advantages. First, epigenetic information is more stable than RNA and protein, and is not vulnerable to physical and chemical damage (Issa, 2012). Second, abnormally altered DNA methylation is an early event of carcinogenesis that predates genetic defects and abnormal gene expression. Detection of DNA methylation changes may provide a more timely and accurate assessment of cancer progression (Fleischer et al., 2014;Timp et al., 2014). Third, the development of drugs to reverse epigenetic modification has great potential for cancer treatment (Zhu and Yao, 2009). The large number of samples and data regarding genome-wide methylation sites stored in The Cancer Genome Atlas (TCGA) database provides a source for the identification of identify biomarkers.
Multi-omics integration analysis interrogates previously isolated data relative to degree of DNA methylation and RNA expression, which may provide new insights into the pathogenesis and treatment of gliomas (Ceccarelli et al., 2016;Binder et al., 2019). We explored epigenetic differences between lower-grade glioma (LGG) and glioblastomas (GBM) based on integrated expression profile data and DNA methylation data in the promoter region. Using a series of statistical methods, we constructed a novel epigenetic signature based on eight DNA methylation sites. Furthermore, we evaluated the ability of the signature to predict survival outcomes and its association with immune infiltration in gliomas. Finally, to apply the epigenetic signature to clinical practice, we integrated a series of independent prognostic factors to construct a nomogram.

Datasets Extracted From TCGA and GEO Databases
The DNA methylation data, gene expression data, and corresponding clinical information of glioma patients (including LGG and GBM samples) were obtained from the TCGA database 1 . In addition, we downloaded the LGG dataset (GSE104293) and the GBM dataset (GSE48462) from the GEO database 2 . DNA methylation data were generated using the Infinium Human Methylation 450 BeadChip covering 485,577 DNA methylation sites. For each CpG site, the β value represented the DNA methylation level from 0 (no methylation) to 1 (100% methylation). CpGs in the promoter regions located 2 kb upstream to 0.5 kb downstream from transcription start sites (TSS) covering 145,907 DNA methylation sites were selected for the present study.

Identification of Genes, and DNA Methylation Sites, and the Epigenetic Genes Involved in Glioma Progression
Using fold change >2 or <0.5 and a false discovery rate (FDR) <0.01 as the threshold, gliomas samples were analyzed to identify differentially expressed genes between LGG samples and GBM samples using the edgeR package 3 . Based on the thresholds β value >0.2 and FDR <0.01, different DNA methylation sites were identified using the Limma package 4 (Ritchie et al., 2015). The correlation between differentially upregulated genes and the degree of reduced methylation on sites on promotors of downregulated genes, and conversely, the downregulate differentially expressed genes and the degree of methylation of promotors sites was analyzed according to the standard of correlation <−0.3 and FDR <0.01. Finally, we identified epigenetic genes associated with glioma progression including epigenetically-induced genes (low promoter methylation with high gene expression) and epigenetically-suppressed genes (high promoter methylation with low gene expression).

Functional Enrichment Analysis of Epigenetic Genes
To explore the biological implications of epigenetic genes in glioma progression, functional enrichment analyses including gene ontology (GO) hierarchy and Kyoto Encyclopedia of Gene and Genomes (KEGG) pathways were performed using the Database for Annotation, Visualization, and Integrated Discovery 5 using P < 0.05 as a significance threshold.

Construction of a Prognosis Risk Signature Based on DNA Methylation Sites of Epigenetically-Regulated Genes
Glioma samples were randomly assigned to training and validation data sets using the caret package 6 . A training data set was used to develop a risk model, and a validation data set and a LGG dataset (GSE104293) and GBM dataset (GSE48462) were used to verify the model. Univariate' cox regression was used to screen for DNA methylation sites that had a significant effect on prognosis, and the filter condition was set at P < 0.05. Lasso regression can reduce the complexity of the model by adjusting the parameters of the model to avoid over-fitting by the glmnet package 7 (Simon et al., 2011). The degree of complexity adjustment of LASSO regression was controlled by the parameter λ-the greater the λ, the greater the punishment for the linear model with more variables-to obtain a model with fewer variables. Finally, we used multivariate Cox regression analysis to screen for independent prognostic factors, which could be used to develop a risk prognostic model. A risk score formula was defined by the sum of the products of the β value of each DNA methylation site and the correlation coefficient. Based on the median value of risk, patients are divided into two groups: high-and low-risk groups. Finally, we evaluated the prediction performance of the signature in the internal and external validation datasets (LGG: GSE104293 and GBM: GSE48462).

Gene Set Enrichment Analysis
Gene set enrichment analysis (GSEA) was preformed to analyze differences between high-and low-risk groups using the "c2 cp keg v7.0 symbols collection" as the reference gene set (Subramanian et al., 2005). A FDR <0.05 after performing 1,000 permutations was considered to significantly enriched.

ESTIMATE Algorithm
The ESTIMATE algorithm was applied to analyze the expression characteristics of specific genes in immune cells and stromal cells, and to calculate immune and stromal scores to evaluate the degree of invasion of non-tumor cells and tumor purity, and to compare differences between high-and low-risk groups (Bindea et al., 2013).

CIBERSORT Algorithm
CIBERSORT, a deconvolution algorithm, evaluates the gene composition of each immune cell by calculating the expression level of each gene in each immune cell. In other words, the adjusted expression profiles of complex tissues were used to predict the proportion of 22 tumor immune infiltrating cell (TIICs) types (Newman et al., 2015). The r package of CIBERSORT was used to calculate each the proportion of TIIC types in each sample, and non-conforming samples were filtered out when the P-value was <0.05. Using the filtered data, the types of TIIL cells were compared in the high-and low-risk groups.

Construction and Assessment of the Nomogram
To apply the model to the clinic, we integrated a series of indicators to construct a nomogram. We used univariate Cox regression to analyze a series of indicators (age, sex, IDH mutation status, ATRX mutation status, 1p/19q co-deletion status, promoter methylation status of MGMT, radiation therapy history, and glioma subtype) to screen for prognostic factors associated with the overall survival (OS) of patients with glioma, and then multivariate Cox regression was performed to screen for independent prognostic factors. Finally, a nomogram was constructed to predict the 1-, 3-, and 5-year OS rates of glioma patients using the rms package 8 . We constructed receiver operating characteristic (ROC) curves to assess the 1-, 3-, and 5-year prediction performance of the nomogram.

Immunohistochemistry
Glioma paraffin-embedded tissues were collected from the Affiliated Cancer Hospital of Guangzhou Medical University (LGG, n = 16; high-grade gliomas, n = 17). Immunohistochemical staining of the target proteins was performed on paraffinembedded tissues using an anti-SPON2 antibody (Proteintech, #20513-1-AP), anti-IFI44 antibody (Proteintech, #7233-1-AP), anti-CD68 antibody (Abcam, #ab125212), anti-CD206 antibody (CST, #91992S). Then, the slices were stained with diaminobenzidine (DAB) and counterstained with hematoxylin. The scoring of tumor cells was as follows: negative (0), yellowish (1-4), light brown (5-8), and dark brown (9-12). The differences in expression of four genes in low-grade and high-grade gliomas were compared by unpaired Student's t-test, and their association was evaluated by Pearson's correlation analysis. This study was approved by the Ethics Committee of the Affiliated Cancer Hospital of Guangzhou Medical University. Informed consent was obtained from each patient.

Statistical Analysis
All statistical analyses were carried out using R software. The Chi-square test was used to analyze the association between the epigenetic signature and common clinical features. The prediction accuracy of the signature was verified by Kaplan-Meier (K-M) analysis and ROC curve analysis. P < 0.05 was considered statistically significant. The area under the curve (AUC) value was used as the evaluation standard of accuracy.

Data Collection and Difference Analysis
A total of 689 glioma samples, consisting of 534 LGG and 155 GBM samples, with information on 485,577 DNA methylation sites was obtained from TCGA. Because DNA methylation in promoter regions significantly affects gene expression, we analyzed the TSS200 and TSS1500 sites located in transcription start sites, and as a result, 145,907 methylation sites were filtered for subsequent study. As shown in the workflow diagram (Figure 1), compared with the LGG samples, 20,220 differential DNA methylation sites were screened from GBM samples (7087 upregulated DNA methylation sites and 13,134 downregulated DNA methylation sites). In addition, we analyzed differentially expressed genes in glioma patients, including 529 LGG samples and 173 GBM samples. A total of 4,134 different expression genes were obtained between LGG samples and GBM samples based on above screening criteria (2,222 upregulated genes and 1,912 downregulated genes).

Identification of Epigenetic Genes Using Correlation Analysis
The degree of gene expression is regulated by DNA methylation. Hypermethylation inhibits the expression of downstream genes, while hypomethylation promotes the expression of downstream genes. Overall, 592 samples with both expression profile and methylation data were selected for follow-up correlation analysis.
First, we considered the intersection of highly methylation genes and downregulated genes, and the intersection of low methylation genes and upregulated genes. Next, using the standard of correlation <−0.3 and FDR < 0.01, a group of 407 epigenetic genes was identified, including 319 epigenetically-induced genes (low promoter methylation with high gene expression) and 88 epigenetically-suppressed genes (high promoter methylation with low gene expression).

GO and Pathway Enrichment Analysis of Epigenetic Genes
To further explore the function of the screened epigenetic genes, the online software DAVID was used for GO analysis of epigenetic genes. GO analysis results classified epigenetic genes functions into three functional groups: biological processes, cellular components, and molecular function. We focused on the biological processes of epigenetically-regulated genes. The results showed that epigenetically-induced genes were enriched in 'positive regulation of I-kappaB kinase/NF-kappaB signaling' , 'regulation of apoptotic process' , 'innate immune response' , and 'inflammatory response' (Figure 2A), while epigenetically suppressed genes were enriched in 'peripheral nervous system

Establishment of the Epigenetic Signature and Validation
In total, 679 glioma samples with clinical follow-up information were randomly distributed to a training data set (n = 340) and a validation data set (n = 339). Clinical information including age, sex, IDH and ATRX mutation status, 1p/19q status, and treatment conditions are summarized in Supplementary Table 1. Statistical methods were used to build a risk prognosis signature in the training set. DNA methylation sites were identified by univariate Cox analysis, and a total of 629 methylation sites were obtained (FDR <0.01). Next, the selected DNA methylation sites were analyzed by LASSO analysis, and 18 key DNA methylation sites were found (Supplementary Figures 1A,B). Finally, eight DNA methylation sites were screened by multivariate Cox regression, stepwise regression, and screening, which can be used to construct an optimal prognosis signature. Information relative to the eight methylation sites is shown in Table 1 and Figure 3A.
The genes corresponding to the eight DNA methylation sites were spondin 2 (SPON2), non-SMC condensin I complex subunit G (NCAPG), interferon induced protein 44 (IFI44), S100 calcium binding protein A2 (S100A2), collagen type XXII alpha 1 chain (COL22A1), CD200 receptor 1 (CD200R1), interferon gamma receptor 2 (IFNGR2), and derlin 3 (DERL3). The degree of DNA methylation of eight methylation sites and corresponding gene expression in high-and low-risk groups are shown in Figures 3B-I. Risk score = −1.39897577 × β value of cg01784327 -2.456105464 × β value of cg02810967 -1.334605549 × β value of cg07107453 -1.525164679 × β value of cg13997435 −1.174403511 × β value of cg16407323 -2.348936398 × β value of cg17638468 −1.200080898 × β value of cg24865779 − 1.071625831 × β value of cg25940946. The eight identified CpGs were all positive prognostic factors for gliomas. The distributions of the risk score, survival status, and β values of DNA methylation sites of glioma patients in the training data set and validation data set are shown in Figures 4A,B. K-M analysis showed that there was a significant correlation between the risk score and OS (p < 0.001), the AUC was 0.918 ( Figure 4C).
In the validation data set, the survival time of highrisk gliomas patient was significantly reduced, with an AUC of 0.874 ( Figure 4D). The percentage of high-and lowrisk gliomas in different subtypes from TCGA database is shown in Supplementary Table 2 (GBM: Classical, G-CIMP, Mesenchymal, Neural, Proneural; LGG: IDH status, ATRX status, 1p/19q status). In the LGG (GSE104293) and GBM (GSE48462) datasets, the OS of the high-risk group was worse, and the AUC values were 0.948 and 0.851, respectively (Figures 4E,F). In addition, we also analyzed the predictive ability of our epigenetic signature across different statuses of IDH mutation, 1p/19q codeletion, and ATRX and MGMT methylation states. K-M analysis and ROC curve results showed that the survival rates of high-risk   glioma patients were all significantly reduced, which indicated that the model had good independent prognostic ability and reliability ( Supplementary Figures 2A-D). The above results indicated that our epigenetic signature has good sensitivity and specificity in predicting the prognosis of glioma patients.

Association of the Epigenetic Signature With Clinical and Molecular Features in Gliomas
The distribution of clinical and molecular features in high-and low-risk groups and the results of chi-square test are shown in Supplementary Table 3. The heatmap shows the distribution of clinical and molecular features in high-and low-risk groups ( Figure 5A). The results showed that except for sex, there were significant differences in age, glioma type, IDH mutation status, ATRX mutation status, 1p/19q status, MGMT promoter methylation status, and survival status between the high-and low-risk groups. Samples with low risk factors (MGMT promoter methylated, mutant-type IDH, mutant-type ATRX, and 1p/19q co-deletion) had lower risk scores.

Comparison of the Epigenetic Signature With Other Known Prognostic Biomarkers
With the development of genetics and molecular biology, many molecular biomarkers have been developed for gliomas, such as IDH mutation, ATRX mutation, MGMT methylation status, and 1p/19q status, but they have limitations (Ludwig and Kornblum, 2017). In recent years, numerous prognostic signatures have been developed based on multiple-level molecular data.

Zeng et al. (2019) constructed a prognostic model for
LGG using DNA damage repair-related genes. Their prognostic model was composed of six genes defined using a weighted gene co-expression network analysis and COX regression analysis . In addition to using DNA methylation sites, lncRNAs were also used to construct prognostic models (Yin et al., 2018;Li et al., 2019). Although previous studies have also used DNA methylation sites to build a prognostic risk model, there was no integrated analysis of DNA methylation and gene expression, or use of methylation sites of epigenetic genes to construct a signature. The ROC curve analysis showed that our signature consisting of eight DNA methylation sites had the highest accuracy and was more accurate than the signature composed of the eight corresponding genes alone ( Figure 5B).

Establishment of a Nomogram for the Prognosis of Patients With Glioma
Age, sex, IDH mutation status, ATRX mutation status, 1p/19q codeletion status, promoter methylation status of MGMT, radiation therapy history, gliomas type, and risk score were significantly related to OS in the TCGA cohort based on the results from the univariate analysis. Through multivariate analysis of the above factors, age, sex, 1p/19q status, glioma type, and the risk score remained independent and stable prognostic factors (p < 0.05) in the cohort ( Table 2). A prognostic nomogram based on the independent factors (age, sex, 1p/19q status, gliomas type, and risk score) was constructed ( Figure 5C). ROC analysis was performed to evaluate the performance of the nomogram in predicting the prognosis of patients in 1-, 3-and 5-years ( Figure 5D).

Comparison of the Results of the GSEA Analysis Between High-and Low-Risk Groups
We used GSEA analysis to compare the KEGG pathways involved in high-and low-risk groups, the results showed that mismatch repair, focal adhesion, leukocyte transendothelial migration, apoptosis, pathways in cancer, and the p53 signaling pathway were significantly enriched in the high-risk group (Figure 6A). Associations Between the Epigenetic Signature and Immune-Checkpoint Blockade Immunotherapy-Related Signature To study the potential role of the epigenetic signature in immunecheckpoint blockade (ICB) immunotherapy, correlation analysis was carried out with these known immunotherapy-related signatures. The results showed that the epigenetic signature was positively correlated with expression of PD-1, PD-L1, PD-L2, CTLA-4, and TMB ( Figure 6B). In addition, we also analyzed the correlation between the eight methylation sites and immunotherapy-related signatures, as well as the correlation between the eight genes and immunotherapy-related signatures, as shown in Supplementary Figures 3A,B. In summary, although our eight-DNA methylation prognostic signature was constructed to accurately predict the prognosis of gliomas patients, it may also play a potential role in defining the application of ICB immunotherapy for glioma patients.

Epigenetic Signature Was Associated With the Tumor Microenvironment
The immune score, stromal score, and ESTIMATE score in the high-risk group were higher than those in the low-risk group, while the tumor purity in the high-risk group was lower than that in the low-risk group (Figure 6C). K-M analysis showed that a poorer prognosis was associated with higher immune scores, stromal scores, ESTIMATE scores, and lower tumor purity ( Figure 6D). The Immune, stromal, and ESTIMATE scores were positively correlated with the risk score, while tumor purity was negatively correlated with risk score (| correlation| > 0.20, P < 0.05), as shown in Figure 6E. In addition, we analyzed the association between the eight methylation sites, their corresponding genes, and the immune, stromal, and ESTIMATE scores, and tumor purity (Supplementary Figures 4A,B).

Relationship Between the Epigenetic Signature and the Proportion of Immune Cell Infiltration
After samples with P < 0.05 were excluded, the remaining 216 samples were subjected to CIBERSORT analysis,including 192 LGG samples and 24 GBM samples. Twenty-two subtypes of immune cells were obtained ( Figure 7A). A t-test was used to compare differences between the 22 kinds between the immune cell types in the high-risk and low-risk groups, as shown in Figure 7B. The proportion of CD8+ T cells, CD4+ memory activated T cells, gamma delta T cells, M0 macrophages, M1 macrophages, M2 macrophages, and (B) Correlation analysis between the risk score and immune checkpoint expression. (C) Distribution of stromal scores, immune score, ESTIMATE score, and tumor purity among high-and low-risk glioma patients. (D) Glioma patients were divided into two groups based on Stromal scores, immune score, ESTIMATE score, and tumor purity. K-M survival curve show overall survival of the low score group is longer than high score group, as indicated by the log-rank test. (E) Correlation analysis between the risk score and stromal scores, immune score, ESTIMATE score, and tumor purity. Abbreviations: K-M, Kaplan-Meier; ns, P > 0.05, ****P < 0.0001.
neutrophils were higher in high-risk group, while plasma cells, activated NK cells, monocytes, activated mast cells, and eosinophils were higher in the low-risk group (P < 0.05). Pearson's correlation analysis was used to analyze the correlation between the risk score and the different immune cell proportions. CD8+ T cells, CD4+ memory activated T cells, gamma delta T cells, M0 macrophages, M1 macrophages, and M2 macrophages were positively correlated with risk score, while activated NK cells, monocytes, activated mast cells were negatively correlated with risk score (| correlation| > 0.20, P < 0.05), as shown in Figure 7C. The results of the K-M analysis showed that glioma patients with a low proportion of activated NK cells, activated Mast cells and high percentage of M0 macrophages and M1 macrophages had a shorter OS ( Figure 7D). Finally, we analyzed the correlation between the methylation degree of eight DNA methylation sites and the proportion 22 immune cell subtypes, and the correlation between the expression of the eight corresponding genes and the proportion of the immune cell subtypes (Supplementary  Figures 5A,B).
macrophage marker CD206. The immunostaining score of IFI44 was positively correlated with CD68 and CD206 expression ( Figure 7F). These results suggested that SPON2 may be involved in the recruitment of non-M2 macrophages, while IFI44 may participate in the recruitment of M2 macrophages.

DISCUSSION
Glioma is the most common malignant tumor of the central nervous system, accounting for about 80% of diagnosed cases. Despite the implementation of surgical treatment and adjuvant therapy, the prognosis of patients with glioma is still poor, resulting difficulty to cure and death (Hemmati et al., 2003). Investigating the differences between GBM and LGG can improve the understanding of the glioma development and the development of treatment strategies of brain tumors (Ceccarelli et al., 2016;Binder et al., 2019). Risk classification of gliomas patients is helpful to guide the choice of treatment, while molecular markers have the advantages of providing prognostic information and exploring the mechanism of tumor progression. We integrated multiple DNA methylation sites to construct a more sensitive and specific signature. Using the methylation sites of the screened epigenetic difference genes, a prognostic signature based on eight methylation sites was developed by applying a series of statistical approaches. K-M analysis showed that the OS of high-and low-risk groups could be easily stratified, and ROC analysis verified that the signature had accurate predictive ability. The clinical molecular characteristics and risk scores were analyzed using Cox regression analysis. The prognosis of GBM, older age, male sex, absence of 1p/19q co-deletion, and a high-risk score for gliomas showed worse prognosis, and all these features were independent prognostic factors (HR > 1 and P < 0. 05). We used the above factors to construct a nomogram able to predict the 1-, 3-, and 5-year prognosis of patients with glioma. There are more immunosuppressive factors than proinflammatory factors expressed in gliomas, suggesting that the microenvironment of gliomas is mainly immunosuppressive (Sokratous et al., 2017). Monocytes from the peripheral blood infiltrate to form microglia-like cells, which participate in the renewal of microglia. Under different signal stimulation conditions, microglia/macrophages M0 in gliomas can develop into M1 microglia/macrophages that promote inflammation and inhibit tumor growth, or M2 microglia/macrophages that inhibit inflammatory response and promote tumor growth, with mainly the infiltration of the latter and namely tumor associated macrophages (TAM) (Durafourt et al., 2012;Shapouri-Moghaddam et al., 2018). Although the tumor immune microenvironment plays an important role in tumor initiation and development, there is no effective signature available that evaluates the immune infiltration of gliomas. In addition, the epigenetic signature has shown limited value in the exploration or evaluation of tumor immune subtypes, especially in the immune infiltration (Aichmüller et al., 2020;Wu et al., 2020).
We integrated publicly available gene expression and DNA methylation data to define epigenetic regulatory mechanisms of some genes essential for immune infiltration. Further, this study has emphasized the potential relationship between CpG methylation markers and indicators of immune infiltration in gliomas. The Infinium Human Methylation450 BeadChip can detected more than 450,000 CpG sites, which may serve as a more suitable alternative marker of immune infiltration. To investigate the correlation between our signature and immune infiltration, we used the following approach. First, the GSEA findings revealed that leukocyte transendothelial migration was significantly enriched in the high-risk group. Second, the ESTIMATE algorithm was used to evaluate the degree of immune and stromal cell infiltration and tumor purity of each glioma patient. Pearson's correlation analysis showed that the risk score was positively correlated with the degree of immune cell and stromal cell infiltration, and negatively correlated with tumor purity. Third, the proportion of 22 immune cells subtypes in each glioma patient was evaluated using the CIBERSORT algorithm. Pearson correlation analysis revealed that the risk score was positively correlated with the infiltration of CD8+ T cells, gamma delta T cells, M0 macrophages, M1 macrophages, and M2 macrophages and was negatively correlated with the proportion of activated NK cells, monocytes, and activated mast cells. The high proportion of M0 macrophages and M2 macrophages was negatively correlated with patient prognosis, whereas a high proportion of NK cells and activated mast cells was positively correlated with prognosis.
NK cells have been reported to exert their cytotoxicity by secreting tumor necrosis factor (TNF) and interferon (IFN) to kill susceptible target cells (Moretta and Moretta, 2004;Raulet, 2004). The possible mechanism of mast cells inhibiting tumor are related to the immune response, but the specific roles of their intracellular bioactive cytokines in tumorigenesis and progression is still difficult to determine. Many experimental results have suggested that mast cells can secrete a variety of cytokines, such as TNF-α, interleukin (IL)-3, IL-4, IL-6, and IL-8 (Gordon et al., 1990;Burd et al., 1995). According to our analysis, it can be inferred that NK cells and mast cells can clear tumor tissue through immune activity. In brief, the epigenetic signature can stratify glioma patients into different immune infiltration subgroups ("low immune infiltration type" and "high immune infiltration type"). Therefore, this novel signature may reflect the immune microenvironment and serves as a prognostic marker in gliomas. We also found that the type of immune infiltration of gliomas was related to the status of IDH mutation or 1p/19q co-deletion. Gliomas with mutated IDH or 1p/19q codeletion presented a lower degree of immune infiltration and better prognosis. Amankulor et al. (2017) found that IDH mutation significantly reduced the infiltration of immune cells such as macrophages, microglia, monocytes, and neutrophils, resulting in the suppression of the tumor-associated immune system. However, it is not clear whether the effect of 1p/19q codeletion status on prognosis is mediated by immune regulation, which is worthy of further study. In addition, we verified the association of SPON2 and IFI44 genes with macrophage recruitment and performed immunohistochemical staining on glioma samples for SPON2, IFI44, the macrophage marker CD68, and the M2 macrophage marker CD206. SPON2, IFI44, CD68, and CD206 were highly expressed in high-grade gliomas, and SPON2 was positively correlated with CD68 rather than CD206, whereas IFI44 was positively correlated with both CD68 and CD206 expression. Recent studies have found that SPON2, a member of the Mindin F-pondin family of conserved secretory ECM proteins, plays an important role in immunity by participating in the initiation of the immune response, and acts as an integrin ligand for inflammatory cell recruitment and T cell priming (Jia et al., 2005;Li et al., 2006). Combined with our experimental results, SPON2 may be involved in the recruitment of non-M2 macrophages in glioma. IFI44 gene, whose expression is induced by alpha/beta-IFN, can block extracellular signal-regulated kinase signals, and causes cell cycle arrest by binding to intracellular GTP (Kitamura et al., 1994;Hallen et al., 2007). Further, IFI44, as an immune-related gene, is involved in immune diseases, but the studies in tumors is insufficient (DeDiego et al., 2019). Our study showed that IFI44 was positively correlated with macrophage infiltration in gliomas. The potential role of IFI44 in tumor formation and development needs further experimental verification.
Previous studies have shown that CD200R1 is also an immune-related gene that encodes the OX-2 membrane glycoprotein receptor, which is mainly distributed on the surface of myeloid cells and T lymphocytes (Gorczynski, 2005). CD200 interacts with its receptor CD200R1 to trigger immunosuppressive signals, resulting in macrophage suppression, regulatory T cell induction, cytokine profile switching from Th1 to Th2, and finally the inhibition tumorspecific T-cell immunity (Gorczynski et al., 1999;Vaine and Soberman, 2014). IFN-γ binds to the IFN-γ receptor (heterodimer composed of IFNGR1 and IFNGR2) and activates a downstream signal pathway that can inhibit tumorigenesis and immunomodulatory effects (Parker et al., 2016;Lin et al., 2017). IFN-γ induces the polarization of M1 macrophages, but not M2 macrophages (Duluc et al., 2009). However, it has not been reported whether the other four genes, NCAPG, S100A2, DERL3, and COL22A1, exert any biological role in the interaction between tumor and immune cells. In the future, we will perform additional experiments to verify the biological function in gliomas of the eight genes identified in this study, which may provide additional treatment targets and improve patient management.
Finally, we studied the monitoring role of DNA methylation biomarkers in ICB immunotherapy. The epigenetic signature was also highly correlated with the expression of CTLA4, PD-L1, and PD-1, suggesting that it may be a potential indicator of cancer immune infiltration and may predict the patients' response to immunotherapy drugs. An individualized treatment regimen can be realized through the stratification of patients, which will significantly improve the effects of immunotherapy and prolong the survival of patients.
In conclusion, this study highlighted the potential role of DNA methylation in assessing risk prognosis and monitoring immune infiltration. Our research innovatively used an approach based on an epigenetic signature to classify gliomas into two categories: the "low immune infiltration type" and "high immune infiltration type". In addition, the eight genes in the signature may be involved in glioma progression and may have a significant potential to improve risk stratification, and regulate immune cell infiltration, and ICB immunotherapy response. If further verified, the epigenetic signature may improve the stratification and immunotherapy options of cancer patients in clinical trials.

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.

ETHICS STATEMENT
The studies involving human participants were reviewed and approved by the Ethics Committee of the Affiliated Cancer Hospital of Guangzhou Medical University. Written informed consent to participate in this study was provided by the participants' legal guardian/next of kin.

AUTHOR CONTRIBUTIONS
JZ and GZ: conception and design. GZ and QZ: administrative support. JY, LL, YS, and DZ: provision of study materials or patients. JY, DH, and GW: collection and assembly of data. JZ, NX, and MY: data analysis and interpretation. All authors are final approval of manuscript.