- 1Department of Neurosurgery, Xiangya Hospital, Central South University, Changsha, China
- 2National Clinical Research Center for Geriatric Disorders, Xiangya Hospital, Central South University, Changsha, China
- 3Department of Neurology, Xiangya Hospital, Central South University, Changsha, China
- 4Department of Neurosurgery, Affiliated Nanhua Hospital, Hengyang Medical College, University of South China, Hengyang, China
- 5Department of Rehabilitation Medicine, Hunan Provincial People’s Hospital, Hunan Normal University, Changsha, China
- 6Department of Psychiatry, The Second People’s Hospital of Hunan Province, The Hospital of Hunan University of Chinese Medicine, Changsha, China
- 7Department of Cerebrovascular Surgery, The Second People’s Hospital of Hunan Province, The Hospital of Hunan University of Chinese Medicine, Changsha, China
- 8Department of Clinical Pharmacology, Xiangya Hospital, Central South University, Changsha, China
- 9Department of Neurology, Hunan Aerospace Hospital, Changsha Medical University, Changsha, China
Background: The tumor immune microenvironment significantly affects tumor occurrence, progression, and prognosis, but its impact on the prognosis of low-grade glioma (LGG) patients with epilepsy has not been reported. Hence, the purpose of this study is to explore its effect on LGG patients with epilepsy.
Methods: The data of LGG patients derived from the TCGA database. The level of immune cell infiltration and the proportion of 22 immune cells were evaluated by ESTIMATE and CIBERSORT algorithms, respectively. The Cox and LASSO regression analysis was adopted to determine the DEGs, and further established the clustering and risk score models. The association between genomic alterations and risk score was investigated using CNV and somatic mutation data. GSVA was adopted to identify the immunological pathways, immune infiltration and inflammatory profiles related to the signature genes. The Tumor Immune Dysfunction and Exclusion (TIDE) algorithm and GDSC database were used to predict the patient’s response to immunotherapy and chemotherapy, respectively.
Results: The prognosis of LGG patients with epilepsy was associated with the immune score. Three prognostic DEGs (ABCC3, PDPN, and INA) were screened out. The expression of signature genes was regulated by DNA methylation. The clustering and risk score models could stratify glioma patients into distinct prognosis groups. The risk score was an independent predictor in prognosis, with a high risk-score indicating a poor prognosis, more malignant clinicopathological and genomic aberration features. The nomogram had the better predictive ability. Patients at high risk had a higher level of macrophage infiltration and increased inflammatory activities associated with T cells and macrophages. While the higher percentage of NK CD56bright cell and more active inflammatory activity associated with B cell were present in the low-risk patients. The signature genes participated in the regulation of immune-related pathways, such as IL6-JAK-STAT3 signaling, IFN-α response, IFN-γ response, and TNFA-signaling-via-NFKB pathways. The high-risk patients were more likely to benefit from anti-PD1 and temozolomide (TMZ) treatment.
Conclusion: An immune-related gene signature was established based on ABCC3, PDPN, and INA, which can be used to predict the prognosis, immune infiltration status, immunotherapy and chemotherapy response of LGG patients with epilepsy.
Introduction
As the most common brain malignant tumor, patients with glioma are often complicated with epilepsy (Goldstein and Feyissa, 2018; Jones et al., 2019; Liu B. et al., 2019). Compared with high-grade glioma, patients with low-grade glioma (LGG) have a significantly higher incidence of epilepsy (Huberfeld and Vecht, 2016), which may be related to the slower growth rate and higher IDH-1 mutation frequency of LGG (Cohen et al., 2013; Chen et al., 2018). The effect of epilepsy on the prognosis of LGG is currently unclear. Several studies have suggested that epilepsy is a positive prognostic factor in LGG patients (Blümcke et al., 2004; Krzyszkowski et al., 2004; Danfors et al., 2009). However, Lote et al. (1998) found that epilepsy had no effect on the prognosis after analyzing 379 LGG patients. This indicates that there are other factors affecting the prognosis of LGG patients with epilepsy.
In addition to tumor cells, non-tumor cells including immune cells are also important components of tumors. Immune cells infiltrating tumor tissues, such as peripheral macrophages, microglia, granulocytes, T lymphocytes, and myeloid-derived suppressor cells (Gieryng et al., 2017), create a tumor immune microenvironment that affects tumor development, metastasis, and response to treatment (Ma et al., 2018). However, its influence on LGG patients with epilepsy has not been established.
Hence, we explored the immune microenvironment of LGG patients with epilepsy and evaluated its impact on prognosis through the integration analysis of multi-omics data. We developed a gene signature related to immune infiltration, which can be used to predict the prognosis and therapeutic response of LGG patients with epilepsy.
Materials and Methods
Data Acquisition
The RNA transcript data (Workflow Type: HTSeq-Counts), corresponding clinical information, DNA methylation (Platform: Illumina human methylation 450), somatic mutation (Workflow Type: MuTect2), and copy number variation (CNV) data of 473 LGG patients (WHO grade II and grade III) were downloaded from the TCGA database1. Among them, there were 297 patients with epilepsy and 176 patients without epilepsy. The sample screening process was shown in Supplementary Figure 1.
Inference of Immune Infiltration Level and Immune Cell Abundance
The immune score was calculated by the “ESTIMATE” algorithm (version 1.0.13) (Yoshihara et al., 2013) to infer the level of immune infiltration. The proportions of 22 immune cells were inferred by the “CIBERSORT” algorithm (Newman et al., 2015) using LM22, a leukocyte gene signature matrix2.
Screening of Prognostic Signature Genes
The “DESeq2” R package (version 1.30.1) (Love et al., 2014) was used to search for differentially expressed genes (DEGs), |log(2) fold change| >2 and p < 0.05 as a threshold. The Cox regression analysis was adopted to find the DEGs with independent prognostic value. Finally, the DEGs with the most prognostic value and their regression coefficients (β) were obtained by the LASSO regression analysis.
Consensus Clustering Analysis
The consensus clustering analysis was conducted using the “ConsensusClusterPlus” R package (version 1.54.0) (Wilkerson and Hayes, 2010). The clustering run was permuted by varying the category number k from 2 to 10. The optimal k value was determined based on good consistency within the clusters and a relatively small incremental change in the area under the CDF curve.
Establishment and Evaluation of Risk Score Model
The risk score model was established as follows: . The effectiveness of the risk model was evaluated by comparing the prognosis of high- and low-risk patients and calculating the AUC value of the ROC curve. The correlation between risk score and clinicopathological characteristics was further analyzed. The R package “rms” (version 6.1-0) was used to construct a nomogram.
Analysis of Genetic Alterations
The somatic mutation data were processed by the “maftools” R package (version 2.6.05) (Mayakonda et al., 2018), and the genes with the highest mutation frequency were screened and presented. Genomic Identification of Significant Targets in Cancer (GISTIC) 2.0 (Mermel et al., 2011) was used to analyze CNV data in reference to the Human genome reference consortium h19 derived from GISTIC 2⋅0. The threshold copy number at alteration peaks was obtained from GISTIC analysis.
The GSVA Analysis
The gene set variation analysis (GSVA) (Hänzelmann et al., 2013) was used to establish each sample’s immune signature and the correlation analysis was further carried out. The correlogram and heatmap were drawn using the “corrgram” (version 1.13) and “pheatmap” R packages (version 1.0.12), respectively. The immunological pathways associated with signature genes were identified by GSVA based on the hallmark gene sets (Liberzon et al., 2015).
Therapy Response Prediction
The patients’ response to immune checkpoint blocking therapy was predicted by the Tumor Immune Dysfunction and Exclusion (TIDE) (Jiang et al., 2018) algorithm. The patients’ sensitivity to temozolomide (TMZ) therapy was predicted using the data derived from the GDSC database3. The “pRRophetic” R package (version 0.5) (Geeleher et al., 2014) was used for the estimation of the half-maximal inhibitory concentration (IC50), which represents the drug response.
Statistical Analysis
The statistical analysis was finished in R software (version 3.5.3). Statistical significance between groups was determined using the Student’s t-test. Pearson rank was adopted in the correlation analysis. The “survival” R package (version 3.2-7) and Kaplan-Meier method were used for survival analysis. The AUC value of ROC curve was calculated by the “survival ROC” R package (version 1.0.3). P < 0.05 was considered to be statistically significant.
Results
The Level of Immune Infiltration Affected the Prognosis
The characteristics of all included patients were summarized in Supplementary Table 1. The workflow designed for the study was shown in Figure 1. According to the immune score and whether they had epilepsy, all patients were divided into four groups: C1(high immune score, seizure), C2 (high immune score, non-seizure), C3 (low immune score, seizure), and C4 (low immune score, non-seizure). Here, the immune scores that were >0 were considered high while immune scores that were <0 were considered low (Wen et al., 2020). Through survival analysis, we found that a high immune score suggested a poor prognosis in LGG patients with epilepsy (p = 0.0025, Figure 2A). The baseline characteristics of LGG patients with epilepsy were shown in Supplementary Table 2. We further compared the abundance of 22 immune cells between high and low immune score groups in LGG patients with epilepsy. The results showed that the patients with high immune score had a higher percentage of macrophages (M2), and a higher proportion of naive B cells, activated dendritic cells, neutrophils, activated memory CD4+ T cells, and resting memory CD4+ T cells also appeared in the high immune score group (p < 0.05, Figure 2B). While the patients in the low immune score group had higher percentages of memory B cells, naive CD4+ T cells, resting mast cells, CD8+ T cells, and follicular T helper cells (p < 0.05, Figure 2B). Further Cox regression analysis indicated that activated memory CD4+ T cells and activated dendritic cells were negative factors for prognosis (p < 0.05, Supplementary Table 3).
Figure 2. (A) Comparison of survival rate of low-grade glioma (LGG) subgroups differentiated by seizure and immune score. (B) Comparison of 22 immune cell components between high and low immune score groups in LGG patients with epilepsy. *p < 0.05, **p < 0.01, ***p < 0.001.
Three Prognostic Signature Genes Were Identified
A total of 12278 DEGs were found between the high and low immune score groups in LGG patients with epilepsy. Furthermore, 51 DEGs with independent prognostic value were identified by Cox regression analysis (p < 0.05). Finally, three genes (ABCC3, INA, PDPN) with the most prognostic value were obtained using the LASSO regression (Figures 3A,B) and their coefficients (ABCC3: 0.0314, INA: −0.0018, PDPN: 0.0957) were shown in Figure 3C. Consequently, a prognostic immune-related gene signature was established. We found that the high immune score group had higher expression levels of ABCC3 and PDPN, while the low immune score group showed higher expression of INA (p < 0.001, Figures 3D–F). The high expression of ABCC3 (HR 3.6, p < 0.001) and PDPN (HR 3.73, p < 0.001) was associated with poor prognosis (Figures 3G,H), while the high expression of the INA (HR 0.34, p < 0.001) was a positive prognostic factor (Figure 3I). This indicates that PDPN and ABCC3 are risk genes, and INA is a protective gene for LGG patients with epilepsy.
Figure 3. Identification of signature genes. (A) The cross-validation plot of LASSO regression [the dashed lines signify the optimal log (λ) value]. (B) LASSO coefficients of 51 DEGs with independent prognostic value. (C) LASSO coefficients of the signature genes. (D–F) Expression difference of ABCC3, PDPN, INA between high and low immune score groups. (G–I) Kaplan-Meier plot showing different prognosis based on the differential expression of ABCC3, PDPN, and INA. (J) Correlation between DNA methylation and signature genes expression. ***p < 0.001.
We also explored the relationship between the methylation levels of signature genes and their expression levels in LGG patients with epilepsy. We found that the expression of ABCC3 (R = −0.51, p < 0.001) and PDPN (R = −0.56, p < 0.001) genes was negatively correlated with methylation levels (p < 0.001), while the INA (R = 0.35, p < 0.001) was positively correlated (Figure 3J and Supplementary Figure 2). This suggests that the hypomethylation of ABCC3 and PDPN is present in LGG patients with epilepsy who have high immune score and poor prognosis.
Clustering Model Based on Signature Genes Identified Two Subgroups With Distinct Outcomes
The samples of LGG patients with epilepsy were grouped using consensus clustering analysis to preliminarily explore the prognostic value of signature genes. In consensus clustering analysis, k = 2 was determined as the optimal selection, which had the small incremental change in the area under CDF curve while keeping the maximal consensus within clusters (Supplementary Figures 3A–C). Therefore, we could obtain two patient subgroups, cluster1 and cluster2. And the principal component analysis (PCA) could obviously distinguish the two clusters, which also proved the rationality of the clustering (Supplementary Figure 3D). A significant difference in prognosis was found between cluster1 and cluster2. The overall survival (OS), disease specific survival (DSS), and progression free interval (PFI) of patients in cluster1 were all better than those in cluster 2 (p < 0.0001, Supplementary Figure 3E).
Risk Score Model Based on Signature Genes Had Excellent Predictive Ability in Prognosis
We established a risk score model based on the signature genes. According to the median value of patients’ risk score calculated by the risk score model, the LGG patients with epilepsy were divided into high- and low-risk groups. The prognosis of patients was significantly different between the two groups, and the OS, DSS, and PFI of patients with high risk were significantly worse (p < 0.0001, Figures 4A–C). The AUC values of 3-year and 5-year ROC curves of the risk model were 0.846 and 0.759 in OS, 0.863 and 0.760 in DSS, 0.734 and 0.711 in PFI, respectively (Figures 4D–F). The above results show that the risk model has a good prognosis prediction ability.
Figure 4. Evaluation of the risk score model. (A–C) Comparison of clinical outcomes between high- and low-risk groups using OS, DSS, and PFI as the endpoints. (D–F) ROC curves and AUC values of risk model.
Risk Score Was Correlated With Clinicopathological and Genomic Aberration Features
We conducted a subgroup analysis to explore the correlation between clinicopathological features and risk score. We found significant differences in risk scores within the subgroups distinguished by MGMT promoter status, IDH status, subtype, WHO grade, and 1p19q codel status (p < 0.001), while no significant differences were found within gender and age subgroups (p > 0.05). Patients with unmethylated MGMT promoter, 1p/19q non-codel, classical, and mesenchymal subtype, IDH-wildtype, and higher WHO grade had higher risk scores (Supplementary Figure 4).
Somatic mutations were observed in 137 (94.48%) of the high-risk samples and in 142 (97.93%) samples of the low-risk group. The mutation frequency of IDH1, CIC, and IDH2 genes in high-risk patients was lower than that in low-risk patients (IDH1: 65% vs. 88%, p < 0.001; CIC: 10% vs. 33%, p < 0.001; IDH2 0% vs. 9%, p < 0.001), while high-risk patients had a higher frequency of mutations in EGFR, NF1, and PTEN (EGFR: 14% vs. 0%, p < 0.001; NF1: 11% vs. 1%, p < 0.001; PTEN 10% vs. 1%, p < 0.001) (Figures 5A,B and Supplementary Figure 5).
Figure 5. Analysis of genomic alterations. (A,B) The somatic mutation profiles in high-risk (A) and low-risk (B) groups. (C,D) The analysis of CNV. The deleted (blue) and amplified (red) chromosomal regions in high- panel (C) and low-risk panel (D) groups; q value = 0.25 as the threshold (the green line).
In the analysis of CNV data, several significant amplification regions containing multiple oncogenes, such as 12q14.1(CDK4), 7p11.2 (EGFR), 4q12 (PDGFRA), and 1q32.1 (PIK3C2B), were identified in high-risk patients (Figure 5C). In addition, focal deletion peak 9p21.3 (CDKN2A, CDKN2B) was detected in patients with high risk (Figure 5C). In the low-risk group, the focal amplification and deletion peaks also appeared, but they had significantly lower G values (Figure 5D).
The Risk Score Was an Independent Predictor for Prognosis
We comprehensively analyzed the effect of risk score and clinicopathological characteristics on a clinical outcome to assess whether the risk score has an independent prognostic value. In univariate Cox regression analysis, age, WHO grade, IDH status, subtype, and risk score were identified as survival-related factors (p < 0.001, Supplementary Table 4). Further multivariate Cox regression analysis indicated that age and risk score were independent prognostic factors. Age ≥ 45 (HR 3.29, p < 0.001) was associated with a worse prognosis, while low risk score (HR 0.40, p = 0.026) was a benign prognostic factor (Figure 6A).
Figure 6. The nomogram predicting the OS of LGG patients with epilepsy. (A) Forest map showing the results of multivariate Cox regression analysis. (B) The nomogram was used by summing the points of each variable. (C) The calibration curves showing the predicted and actual observed OS rates. (D) K-M plot showing the difference in OS between high- and low-risk groups. (E) The ROC curves and AUC values of the nomogram.
In addition, we constructed a nomogram that integrated the clinicopathological features and risk score to improve clinical feasibility (Figure 6B). The predicted 3/5-year OS were close to the actually observed ones, as shown in the calibration curve (Figure 6C). According to the median value of patients’ points calculated by the nomogram, we could also obtain two groups of patients, the high-risk group and the low-risk group. Survival analysis suggested that high-risk patients had worse OS than low-risk patients (p = 0.00037, Figure 6D). The AUC values of ROC curves from the nomogram were 0.917 at predicted 3-year OS and 0.783 at 5-year OS, both of which were better than those obtained by using the risk model alone (Figure 6E). Therefore, the nomogram has better predictive power compared with the risk model.
Immunological Pathways, Immune Infiltration and Inflammatory Profiles Related to Signature Genes
We explored the immunological activities related to the signature genes. We first carried out the GSVA based on specific marker genes of 24 immune cells and seven metagenes representing different types of inflammation and immune responses (Rody et al., 2009). Further correlation analysis showed that macrophages infiltration was positively correlated with risk score (p = 0.002; Figures 7A,B and Supplementary Table 5). While B cells, T cells, T helper (Th) cells, follicular helper T cell (TFH), Th17 cells, Th1 cells, Th2 cells, gamma delta T Cells (Tgd), central memory T cell (Tcm), effector memory T cell (Tem), NK CD56dim cells, NK CD56bright cells, mast cells, dendritic cell (DC), activated DC (aDC), immature DC (iDC), cytotoxic cells, and neutrophils infiltration were negatively correlated with risk score, especially the NK CD56bright cells (p < 0.05; Figures 7A,B and Supplementary Table 5). In terms of inflammatory activities, MHC-II, STAT1, LCK, MHC-I, HCK, and Interferon were positively correlated with risk score, while IgG, a marker for B cells, was negatively correlated (p < 0.001; Figures 7C,D and Supplementary Table 6).
Figure 7. The GSVA results. (A,C) Correlograms showing the correlation between immune cell infiltration, inflammation activities, and risk score. (B,D) Heatmaps displaying the activities of immune cells panel (B) and inflammation panel (D).
We also conducted GSVA to identify the immunological pathways associated with the signature genes. Compared with the low-risk group, the IL6-JAK-STAT3 signaling, IFN-α response, IFN-γ response, and TNFA-signaling-via-NFKB pathways were significantly active in patients with high-risk (p < 0.05, Supplementary Figure 6).
Risk Score Was Correlated With Immunotherapy and Chemotherapy Response
Cancer growth and progression are associated with immunosuppression, and it is currently believed that the immune checkpoint pathways are the major mechanism of tumor immune resistance (Pardoll, 2012). The immune checkpoint blockade has become a promising immunotherapy method for various tumors including glioma (Pardoll, 2012; Hodges et al., 2017). In this study, we found that immune checkpoint molecules CD724 (PD-L1), CTL4, LAG3, and PDCD1 (PD1) were expressed at higher levels in high-risk patients (p < 0.01, Figure 8A). We further predicted the patients’ response to anti-CTLA4 and anti-PD1 therapy and found that high-risk patients were more likely to respond to anti-PD1 therapy than low-risk patients (p = 0.01, Figure 8B). No significant difference was observed in response to anti-CTLA4 therapy between the two groups (p > 0.05, Figure 8B). Moreover, a higher tumor mutation burden (TMB) was observed in high-risk patients (p < 0.01, Figure 8C).
Figure 8. Prediction of response to immunotherapy and chemotherapy. (A) Expression of immune checkpoints in high- and low-risk groups. (B) The analysis of response to immunotherapy. (C) Comparison of TMB between the two groups. (D) The estimated IC50 value of temozolomide (TMZ) in the two groups. **p < 0.01, ***p < 0.001.
Meanwhile, we also predicted the patients’ sensitivity to TMZ therapy, which is currently the most commonly used chemotherapeutics for glioma. We trained a predictive model using cell line data derived from the GDSC database. With this predictive model, the IC50 value of each patient’s TMZ was estimated. We found that the estimated IC50 value in low-risk patients was higher compared with high-risk patients (p < 0.001, Figure 8D). This indicates that patients with high risk may be more sensitive to TMZ therapy.
Discussion
In this study, we explored the effect of the immune microenvironment on the prognosis of LGG patients with epilepsy. We found that the immune score was associated with survival, and patients with a high immune score had a poor prognosis. Moreover, the difference of immune cell infiltration between different immune score groups was observed, and a higher proportion of M2 macrophages appeared in high immune score patients. The M2 macrophage infiltration subverts the host’s adaptive immune response and fosters a tumor milieu ripe for angiogenesis, migration, and metastasis in the glioma (Michelson et al., 2016; Liu et al., 2018). These findings indicate that the immune microenvironment affects the prognosis of LGG patients with epilepsy.
Further, three prognostic DEGs (ABCC3, INA, and PDPN) were identified between high and low immune score groups in LGG patients with epilepsy. Hence, a prognostic immune-related gene signature was established. The ATP binding cassette subfamily C member 3 (ABCC3) gene belongs to the ABCC gene family encoding the ABC transporters and participates in transporting ATP-dependent substances across lipid membranes (Chen and Tiwari, 2011; Mao et al., 2019). A high expression level of the ABCC3 gene is associated with various tumors’ poor prognosis, such as gastric cancer, breast cancer, and pancreatic cancer (Adamska et al., 2019; Mao et al., 2019; Wang et al., 2020). Podoplanin (PDPN) is a glycoprotein receptor that affects cell behavior by interacting with other proteins and has various physiological functions (Krishnan et al., 2018). The PDPN is crucial for the development of the lymphatic system, lungs, and heart (Astarita et al., 2012). It is also a malignant tumor promoter, which promotes tumor migration, invasion, lymphangiogenesis, and angiogenesis (Astarita et al., 2012; Quintanilla et al., 2019); Its up-regulation in tumor tissues is correlated with the poor prognosis of mesothelioma, melanoma, and squamous cell carcinoma (Watanabe et al., 2020). Alpha-internexin (INA) gene encodes a neurofilament interacting protein and is involved in neuronal development and axonal outgrowth (Ducray et al., 2011; Rajmohan et al., 2020). The reduced/loss of its expression is related to tumor invasion and suggests a poor prognosis in patients with gastroenteropancreatic neuroendocrine neoplasms (GEP-NENs) (Wang et al., 2018). In this study, the high expression of ABCC3 and PDPN genes was found in the high immune score group and was related to the poor prognosis of LGG patients with epilepsy. And the high expression of INA gene was observed in patients with low immune score and suggests a better prognosis. We also found that the expression level of signature genes was regulated by DNA methylation in LGG patients with epilepsy. DNA methylation is associated with gene silencing and tumorigenesis, and tumors are often accompanied by hypomethylation of oncogenes and hypermethylation of tumor suppressor genes (Morgan et al., 2018). This also provides new potential targets for the regulation of glioma epigenetics.
We established a risk score model based on the above signature genes, and the high risk-score suggested a poor prognosis. The risk score was determined as a factor with an independent prognostic value. The 1p/19q non−codel, IDH-wildtype, classical and mesenchymal subtype, MGMT promotor unmethylated, and higher WHO grade, which are the malignant clinicopathological features in glioma (Verhaak et al., 2010; Weller et al., 2012), were associated with a higher risk score. In addition, the frequency of IDH1, CIC, and IDH2 mutations suggested a positive prognosis (Appin and Brat, 2014; Miller et al., 2017; Liu Z. et al., 2019) were lower in high-risk patients, while the malignancy driving mutations (EGFR, NF1, and PTEN) (Appin and Brat, 2014; Asif et al., 2019) were more frequent. And the deletion peaks of tumor suppressor genes (CDKN2A, CDKN2B) and amplification peaks of oncogenes (EGFR, CDK4, PDGFRA, PIK3C2B) (Rao et al., 2010) were also detected in patients with high risk. These results suggest that this gene signature is a reliable prognostic predictor for LGG patients with epilepsy. Additionally, the nomogram constructed by integrating the clinicopathological features and risk score had a better prognostic prediction ability, and further improves the clinical practicality.
The relationship between the immune microenvironment and gene signature was further explored. The high- and low-risk patients had different immune cell infiltration and inflammatory profiles. Macrophage infiltration level was higher in high-risk patients, while low-risk patients had higher B cells, T cells, Tgd, Th1 cells, T helper cells, Th2 cells, TFH, Tem, Tcm, Th17 cells, NK CD56bright cells, NK CD56dim cells, mast cells, iDC, aDC, DC, cytotoxic cells, and neutrophils infiltration, especially the NK CD56bright cells. Macrophages are recruited to the glioma environment and create a supportive stroma for neoplastic cell invasion and expansion through the release of various growth factors and cytokines, and the number of macrophages increases with increasing malignancy grade (Hambardzumyan et al., 2016). The CD56bright cells are a functional subset of NK cells, which are characterized by low cytotoxicity and high cytokine production, and play an immunoregulatory role (Shevtsov and Multhoff, 2016). The infiltration of NK CD56bright cells is negatively correlated with the malignant degree of gastric carcinomas, lung tumors, and papillary thyroid cancer (Gogali et al., 2013). However, the effect of NK CD56bright cells on glioma has not been sufficiently explored. The infusion of IL15, which is required for survival, expansion, and activation of NK cells, significantly expands the subpopulation of CD56 bright NK Cells compared with the CD56dim NK Cells (Dubois et al., 2017), another functional subset of NK cells with high cytotoxic activity and low cytokine release properties (Shevtsov and Multhoff, 2016). Moreover, infusion of IL-15 in glioma model mice significantly increases the infiltration of NK cells into tumor and reduces tumor growth (Garofalo et al., 2015). Combined with the results of this study, we believe that NK CD56bright cells play an important role in glioma immunity, and its infiltration improves the prognosis. In addition, the high-risk patients had more active inflammatory responses related to T cells and macrophages, while patients with low risk had more active inflammatory response associated with B cells. We further explored the mechanism of signature genes involved in the regulation of the immune microenvironment, and found that the signature genes participated in the regulation of IL6-JAK-STAT3 signaling, IFN-α response, IFN-γ response, and TNFA-signaling-via-NFKB pathways, which are onco-immunological pathways (Wang et al., 2009; De Simone et al., 2015; Mostofa et al., 2017; Silginer et al., 2017; Zhu et al., 2020).
Immune checkpoint blockade is a promising treatment for glioma patients. However, not all patients can benefit from it (Hodges et al., 2017). Therefore, we wonder whether this gene signature can be used as a predictive marker for the clinical response to immune checkpoint blocking therapy in LGG patients with epilepsy. The results showed that high-risk patients were more likely to benefit from anti-PD1 treatment. This may be related to their high expression of PD1/PD-L1, and the active IL6-STAT3 signaling and IFN-γ response pathways in high-risk patients can promote the expression of PD1/PD-L1 (Zhang et al., 2016; Qian et al., 2018; Litak et al., 2019). Also, a higher TMB was found in high-risk patients, and high TMB is currently considered to indicate a relatively better response rate to immunotherapy (Hellmann et al., 2018). We also found that high-risk patients might be more sensitive to TMZ therapy, which may be related to their active IFN pathway down-regulating MGMT level and increasing the sensitivity of glioma cells to TMZ (Shen et al., 2015).
In general, we established an immune-related gene signature, which can be used to predict the prognosis of LGG patients with epilepsy, and infer their immune infiltration status, response to immunotherapy and chemotherapy. Moreover, this study also deepens our understanding of the immune microenvironment in glioma. Of course, there are some limitations in this study. First of all, only TCGA data was used for analysis in this paper, and the sample size was relatively limited. It may be necessary to further expand the sample size and verify the results of this study in other databases. Secondly, due to the difficulty in obtaining clinical samples from LGG patients with epilepsy, no further in vitro experiments were conducted in this study to verify the results. In the end, the specific mechanisms of signature genes involved in glioma immunity need further exploration. The p values obtained from the analysis in this paper were summarized in Supplementary Table 7.
Data Availability Statement
Publicly available datasets were analyzed in this study. This data can be found here: https://portal.gdc.cancer.gov/.
Author Contributions
ZX and QC conceived and designed the research. QC and WD drafted and revised the manuscript. ZX, QC, and WD analyzed the data. KL, CL, SH, WY, BY, and HC acquired and interpreted the data. All authors contributed to the article and approved the submitted version.
Funding
This work was supported by Science Foundation of AMHT Group (No. 2020YK10), the Key Research and Development Programs of Hunan Provincal Science and Technology Department (No. 2017SK2081), China Postdoctoral Science Foundation (No. 2018M633002), Hunan Provincial Natural Science Foundation of China (No. 2018JJ3838), Hunan Provincial Health Committee Foundation of China (C2019186), and Xiangya Hospital Central South University Postdoctoral Foundation.
Conflict of Interest
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
The handling editor declared a shared affiliation, though no other collaboration, with several of the authors WD, QC, and WY at the time of the review.
Supplementary Material
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fcell.2021.686909/full#supplementary-material
Footnotes
References
Adamska, A., Domenichini, A., Capone, E., Damiani, V., Akkaya, B. G., Linton, K. J., et al. (2019). Pharmacological inhibition of ABCC3 slows tumour progression in animal models of pancreatic cancer. J. Exp. Clin. Cancer Res. 38:312. doi: 10.1186/s13046-019-1308-7
Appin, C. L., and Brat, D. J. (2014). Molecular genetics of gliomas. Cancer J. 20, 66–72. doi: 10.1097/PPO.0000000000000020
Asif, S., Fatima, R., Krc, R., Bennett, J., and Raza, S. (2019). Comparative proteogenomic characterization of glioblastoma. CNS Oncol. 8:CNS37. doi: 10.2217/cns-2019-0003
Astarita, J. L., Acton, S. E., and Turley, S. J. (2012). Podoplanin: emerging functions in development, the immune system, and cancer. Front. Immunol. 3:283. doi: 10.3389/fimmu.2012.00283
Blümcke, I., Luyken, C., Urbach, H., Schramm, J., and Wiestler, O. D. (2004). An isomorphic subtype of long-term epilepsy-associated astrocytomas associated with benign prognosis. Acta Neuropathol. 107, 381–388. doi: 10.1007/s00401-004-0833-3
Chen, D. Y., Chen, C. C., Crawford, J. R., and Wang, S. G. (2018). Tumor-related epilepsy: epidemiology, pathogenesis and management. J. Neurooncol. 139, 13–21. doi: 10.1007/s11060-018-2862-0
Chen, Z.-S., and Tiwari, A. K. (2011). Multidrug resistance proteins (MRPs/ABCCs) in cancer chemotherapy and genetic diseases. FEBS J. 278, 3226–3245. doi: 10.1111/j.1742-4658.2011.08235.x
Cohen, A. L., Holmen, S. L., and Colman, H. (2013). IDH1 and IDH2 mutations in gliomas. Curr. Neurol. Neurosci. Rep. 13:345. doi: 10.1007/s11910-013-0345-4
Danfors, T., Ribom, D., Berntsson, S. G., and Smits, A. (2009). Epileptic seizures and survival in early disease of grade 2 gliomas. Eur. J. Neurol. 16, 823–831. doi: 10.1111/j.1468-1331.2009.02599.x
De Simone, V., Franzè, E., Ronchetti, G., Colantoni, A., Fantini, M. C., Di Fusco, D., et al. (2015). Th17-type cytokines, IL-6 and TNF-α synergistically activate STAT3 and NF-kB to promote colorectal cancer cell growth. Oncogene 34, 3493–3503. doi: 10.1038/onc.2014.286
Dubois, S., Conlon, K. C., Müller, J. R., Hsu-Albert, J., Beltran, N., Bryant, B. R., et al. (2017). IL15 infusion of cancer patients expands the subpopulation of cytotoxic CD56 NK cells and increases NK-cell cytokine release capabilities. Cancer Immunol. Res. 5, 929–938. doi: 10.1158/2326-6066.CIR-17-0279
Ducray, F., Mokhtari, K., Crinière, E., Idbaih, A., Marie, Y., Dehais, C., et al. (2011). Diagnostic and prognostic value of alpha internexin expression in a series of 409 gliomas. Eur. J. Cancer 47, 802–808. doi: 10.1016/j.ejca.2010.11.031
Garofalo, S., D’Alessandro, G., Chece, G., Brau, F., Maggi, L., Rosa, A., et al. (2015). Enriched environment reduces glioma growth through immune and non-immune mechanisms in mice. Nat. Commun. 6:6623. doi: 10.1038/ncomms7623
Geeleher, P., Cox, N., and Huang, R. S. (2014). pRRophetic: an R package for prediction of clinical chemotherapeutic response from tumor gene expression levels. PLoS One 9:e107468. doi: 10.1371/journal.pone.0107468
Gieryng, A., Pszczolkowska, D., Walentynowicz, K. A., Rajan, W. D., and Kaminska, B. (2017). Immune microenvironment of gliomas. Lab Invest. 97, 498–518. doi: 10.1038/labinvest.2017.19
Gogali, F., Paterakis, G., Rassidakis, G. Z., Liakou, C. I., and Liapi, C. (2013). CD3(-)CD16(-)CD56(bright) immunoregulatory NK cells are increased in the tumor microenvironment and inversely correlate with advanced stages in patients with papillary thyroid cancer. Thyroid 23, 1561–1568. doi: 10.1089/thy.2012.0560
Goldstein, E. D., and Feyissa, A. M. (2018). Brain tumor related-epilepsy. Neurol. Neurochir. Pol. 52, 436–447. doi: 10.1016/j.pjnns.2018.06.001
Hambardzumyan, D., Gutmann, D. H., and Kettenmann, H. (2016). The role of microglia and macrophages in glioma maintenance and progression. Nat. Neurosci. 19, 20–27. doi: 10.1038/nn.4185
Hänzelmann, S., Castelo, R., and Guinney, J. (2013). GSVA: gene set variation analysis for microarray and RNA-seq data. BMC Bioinformatics 14:7. doi: 10.1186/1471-2105-14-7
Hellmann, M. D., Ciuleanu, T.-E., Pluzanski, A., Lee, J. S., Otterson, G. A., Audigier-Valette, C., et al. (2018). Nivolumab plus ipilimumab in lung cancer with a high tumor mutational burden. N. Engl. J. Med. 378, 2093–2104. doi: 10.1056/NEJMoa1801946
Hodges, T. R., Ott, M., Xiu, J., Gatalica, Z., Swensen, J., Zhou, S., et al. (2017). Mutational burden, immune checkpoint expression, and mismatch repair in glioma: implications for immune checkpoint immunotherapy. Neuro Oncol. 19, 1047–1057. doi: 10.1093/neuonc/nox026
Huberfeld, G., and Vecht, C. J. (2016). Seizures and gliomas–towards a single therapeutic approach. Nat. Rev. Neurol. 12, 204–216. doi: 10.1038/nrneurol.2016.26
Jiang, P., Gu, S., Pan, D., Fu, J., Sahu, A., Hu, X., et al. (2018). Signatures of T cell dysfunction and exclusion predict cancer immunotherapy response. Nat. Med. 24, 1550–1558. doi: 10.1038/s41591-018-0136-1
Jones, D. T. W., Bandopadhayay, P., and Jabado, N. (2019). The power of human cancer genetics as revealed by low-grade gliomas. Annu. Rev. Genet. 53, 483–503. doi: 10.1146/annurev-genet-120417-031642
Krishnan, H., Rayes, J., Miyashita, T., Ishii, G., Retzbach, E. P., Sheehan, S. A., et al. (2018). Podoplanin: an emerging cancer biomarker and therapeutic target. Cancer Sci. 109, 1292–1299. doi: 10.1111/cas.13580
Krzyszkowski, T., Czepko, R., and Betlej, M. (2004). [Prognostic value of epileptic seizures in patients with cerebral gliomas]. Ann. Acad. Med. Stetin. 50, 35–40.
Liberzon, A., Birger, C., Thorvaldsdóttir, H., Ghandi, M., Mesirov, J. P., and Tamayo, P. (2015). The molecular signatures database (MSigDB) hallmark gene set collection. Cell Syst. 1, 417–425. doi: 10.1016/j.cels.2015.12.004
Litak, J., Mazurek, M., Grochowski, C., Kamieniak, P., and Roliński, J. (2019). PD-L1/PD-1 Axis in glioblastoma multiforme. Int. J. Mol. Sci. 20:5347. doi: 10.3390/ijms20215347
Liu, B., Liu, J., Liu, K., Huang, H., Li, Y., Hu, X., et al. (2019). A prognostic signature of five pseudogenes for predicting lower-grade gliomas. Biomed. Pharmacother. 117:109116. doi: 10.1016/j.biopha.2019.109116
Liu, Z., Kuang, W., Zhou, Q., and Zhang, Y. (2018). TGF-β1 secreted by M2 phenotype macrophages enhances the stemness and migration of glioma cells via the SMAD2/3 signalling pathway. Int. J. Mol. Med. 42, 3395–3403. doi: 10.3892/ijmm.2018.3923
Liu, Z., Liu, H., Liu, Z., and Zhang, J. (2019). Oligodendroglial tumours: subventricular zone involvement and seizure history are associated with CIC mutation status. BMC Neurol. 19:134. doi: 10.1186/s12883-019-1362-y
Lote, K., Stenwig, A. E., Skullerud, K., and Hirschberg, H. (1998). Prevalence and prognostic significance of epilepsy in patients with gliomas. Eur. J. Cancer 34, 98–102. doi: 10.1016/s0959-8049(97)00374-2
Love, M. I., Huber, W., and Anders, S. (2014). Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol. 15:550.
Ma, Q., Long, W., Xing, C., Chu, J., Luo, M., Wang, H. Y., et al. (2018). Cancer stem cells and immunosuppressive microenvironment in glioma. Front. Immunol. 9:2924. doi: 10.3389/fimmu.2018.02924
Mao, X., He, Z., Zhou, F., Huang, Y., and Zhu, G. (2019). Prognostic significance and molecular mechanisms of adenosine triphosphate-binding cassette subfamily C members in gastric cancer. Medicine (Baltimore) 98:e18347. doi: 10.1097/MD.0000000000018347
Mayakonda, A., Lin, D.-C., Assenov, Y., Plass, C., and Koeffler, H. P. (2018). Maftools: efficient and comprehensive analysis of somatic variants in cancer. Genome Res. 28, 1747–1756. doi: 10.1101/gr.239244.118
Mermel, C. H., Schumacher, S. E., Hill, B., Meyerson, M. L., Beroukhim, R., and Getz, G. (2011). GISTIC2.0 facilitates sensitive and confident localization of the targets of focal somatic copy-number alteration in human cancers. Genome Biol. 12:R41. doi: 10.1186/gb-2011-12-4-r41
Michelson, N., Rincon-Torroella, J., Quiñones-Hinojosa, A., and Greenfield, J. P. (2016). Exploring the role of inflammation in the malignant transformation of low-grade gliomas. J. Neuroimmunol. 297, 132–140. doi: 10.1016/j.jneuroim.2016.05.019
Miller, J. J., Shih, H. A., Andronesi, O. C., and Cahill, D. P. (2017). Isocitrate dehydrogenase-mutant glioma: evolving clinical and therapeutic implications. Cancer 123, 4535–4546. doi: 10.1002/cncr.31039
Morgan, A. E., Davies, T. J., and Mc Auley, M. T. (2018). The role of DNA methylation in ageing and cancer. Proc. Nutr. Soc. 77, 412–422. doi: 10.1017/S0029665118000150
Mostofa, A. G. M., Punganuru, S. R., Madala, H. R., Al-Obaide, M., and Srivenugopal, K. S. (2017). The Process And Regulatory Components Of Inflammation In Brain Oncogenesis. Biomolecules 7:34. doi: 10.3390/biom7020034
Newman, A. M., Liu, C. L., Green, M. R., Gentles, A. J., Feng, W., Xu, Y., et al. (2015). Robust enumeration of cell subsets from tissue expression profiles. Nat. Methods 12, 453–457. doi: 10.1038/nmeth.3337
Pardoll, D. M. (2012). The blockade of immune checkpoints in cancer immunotherapy. Nat. Rev. Cancer 12, 252–264. doi: 10.1038/nrc3239
Qian, J., Wang, C., Wang, B., Yang, J., Wang, Y., Luo, F., et al. (2018). The IFN-γ/PD-L1 axis between T cells and tumor microenvironment: hints for glioma anti-PD-1/PD-L1 therapy. J. Neuroinflammation 15:290. doi: 10.1186/s12974-018-1330-2
Quintanilla, M., Montero-Montero, L., Renart, J., and Martín-Villar, E. (2019). Podoplanin in inflammation and cancer. Int. J. Mol. Sci. 20:707. doi: 10.3390/ijms20030707
Rajmohan, K. S., Sugur, H. S., Shwetha, S. D., Pandey, P., Arivazhagan, A., and Santosh, V. (2020). Alpha internexin: a surrogate marker for 1p/19q codeletion and prognostic marker in anaplastic (WHO grade III) gliomas. Neurol. India 68, 832–837. doi: 10.4103/0028-3886.293453
Rao, S. K., Edwards, J., Joshi, A. D., Siu, I. M., and Riggins, G. J. (2010). A survey of glioblastoma genomic amplifications and deletions. J. Neurooncol. 96, 169–179. doi: 10.1007/s11060-009-9959-4
Rody, A., Holtrich, U., Pusztai, L., Liedtke, C., Gaetje, R., Ruckhaeberle, E., et al. (2009). T-cell metagene predicts a favorable prognosis in estrogen receptor-negative and HER2-positive breast cancers. Breast Cancer Res. 11:R15. doi: 10.1186/bcr2234
Shen, D., Guo, C.-C., Wang, J., Qiu, Z.-K., Sai, K., Yang, Q.-Y., et al. (2015). Interferon-α/β enhances temozolomide activity against MGMT-positive glioma stem-like cells. Oncol. Rep. 34, 2715–2721. doi: 10.3892/or.2015.4232
Shevtsov, M., and Multhoff, G. (2016). Immunological and translational aspects of NK cell-based antitumor immunotherapies. Front. Immunol. 7:492. doi: 10.3389/fimmu.2016.00492
Silginer, M., Nagy, S., Happold, C., Schneider, H., Weller, M., and Roth, P. (2017). Autocrine activation of the IFN signaling pathway may promote immune escape in glioblastoma. Neuro Oncol. 19, 1338–1349. doi: 10.1093/neuonc/nox051
Verhaak, R. G. W., Hoadley, K. A., Purdom, E., Wang, V., Qi, Y., Wilkerson, M. D., et al. (2010). Integrated genomic analysis identifies clinically relevant subtypes of glioblastoma characterized by abnormalities in PDGFRA, IDH1, EGFR, and NF1. Cancer Cell 17, 98–110. doi: 10.1016/j.ccr.2009.12.020
Wang, H., Lathia, J. D., Wu, Q., Wang, J., Li, Z., Heddleston, J. M., et al. (2009). Targeting interleukin 6 signaling suppresses glioma stem cell survival and tumor growth. Stem Cells 27, 2393–2404. doi: 10.1002/stem.188
Wang, Y., Chen, Y., Li, X., Hu, W., Zhang, Y., Chen, L., et al. (2018). Loss of expression and prognosis value of alpha-internexin in gastroenteropancreatic neuroendocrine neoplasm. BMC Cancer 18:691. doi: 10.1186/s12885-018-4449-8
Wang, Y., Zhang, X., Zhao, B., Xu, Z., and Lv, Y. (2020). Suspension state promotes drug resistance of breast tumor cells by inducing ABCC3 overexpression. Appl. Biochem. Biotechnol. 190, 410–422. doi: 10.1007/s12010-019-03084-0
Watanabe, N., Kidokoro, M., Tanaka, M., Inoue, S., Tsuji, T., Akatuska, H., et al. (2020). Podoplanin is indispensable for cell motility and platelet-induced epithelial-to-mesenchymal transition-related gene expression in esophagus squamous carcinoma TE11A cells. Cancer Cell Int. 20:263. doi: 10.1186/s12935-020-01328-2
Weller, M., Stupp, R., Hegi, M. E., van den Bent, M., Tonn, J. C., Sanson, M., et al. (2012). Personalized care in neuro-oncology coming of age: why we need MGMT and 1p/19q testing for malignant glioma patients in clinical practice. Neuro Oncol. 14(Suppl. 4), iv100–iv108. doi: 10.1093/neuonc/nos206
Wen, P., Gao, Y., Chen, B., Qi, X., Hu, G., Xu, A., et al. (2020). Pan-cancer analysis of radiotherapy benefits and immune infiltration in multiple human cancers. Cancers (Basel) 12:957. doi: 10.3390/cancers12040957
Wilkerson, M. D., and Hayes, D. N. (2010). ConsensusClusterPlus: a class discovery tool with confidence assessments and item tracking. Bioinformatics 26, 1572–1573. doi: 10.1093/bioinformatics/btq170
Yoshihara, K., Shahmoradgoli, M., Martínez, E., Vegesna, R., Kim, H., Torres-Garcia, W., et al. (2013). Inferring tumour purity and stromal and immune cell admixture from expression data. Nat. Commun. 4:2612. doi: 10.1038/ncomms3612
Zhang, N., Zeng, Y., Du, W., Zhu, J., Shen, D., Liu, Z., et al. (2016). The EGFR pathway is involved in the regulation of PD-L1 expression via the IL-6/JAK/STAT3 signaling pathway in EGFR-mutated non-small cell lung cancer. Int. J. Oncol. 49, 1360–1368. doi: 10.3892/ijo.2016.3632
Keywords: glioma, epilepsy, immune, immunotherapy, prognosis
Citation: Cheng Q, Duan W, He S, Li C, Cao H, Liu K, Ye W, Yuan B and Xia Z (2021) Multi-Omics Data Integration Analysis of an Immune-Related Gene Signature in LGG Patients With Epilepsy. Front. Cell Dev. Biol. 9:686909. doi: 10.3389/fcell.2021.686909
Received: 28 March 2021; Accepted: 22 June 2021;
Published: 16 July 2021.
Edited by:
Lei Deng, Central South University, ChinaReviewed by:
Jun Liu, Yuebei People’s Hospital, ChinaJukun Song, Guizhou Provincial People’s Hospital, China
Weifeng Hong, The First Affiliated Hospital of Guangdong Pharmaceutical University, China
Copyright © 2021 Cheng, Duan, He, Li, Cao, Liu, Ye, Yuan and Xia. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
*Correspondence: Zhiwei Xia, xiazhiwei2011@gmail.com
†These authors have contributed equally to this work