Characterization of Immune-Related Genes and Immune Infiltration Features in Epilepsy by Multi-Transcriptome Data

Introduction

Epilepsy, is a highly heterogeneous neurological disease in nature, is characterized by recurrent unprovoked epileptic seizures on account of abnormal excessive or highly synchronous discharge of brain neurons, which are connected with varying degrees of change in the cellular and/or gene translation levels.^1,2 The main clinical manifestations are recurrent convulsions, which seriously afflict the work, daily activities, and physical and mental health.³ According to estimates from the WHO, nearly 50 million people worldwide have active epilepsy with continuing seizures that need treatment.⁴ So far dozens of antiepileptic drugs (AEDs) have been exploited, which all share mechanisms of action rooted in the modulation of function of neurotransmitters and/or of ion channels, exert excellent efficacy on epileptic seizures.⁵ However, approximately one-third of clinic-based patients remain untreated and often suffer from drug-resistant epilepsy (DRE).¹ Hence, it is of vital importance to investigate the molecular characteristics and mechanism of occurrence and development of epilepsy to provide new strategies for the diagnosis and treatment of epilepsy.

In 2017, ILAE (the International League Against Epilepsy) divided epilepsy into six categories, to provide a framework to circumvent such constraints for basic research as well as clinical practice, based on the following etiologies: 1) structural etiology; 2) genetic etiology; 3) infectious etiology; 4) metabolic etiology; 5) immune etiology, and 6) unknown etiology.⁶ At present, the pathogenesis of epilepsy has been widely explored, which mostly focuses on neurotransmission and ion channels (structural etiology). The imbalance of excitatory and inhibitory neural transmission and the formation of aberrant epileptic networks are considered to be the common mechanisms of various epileptic seizures. With the investigation of the pathogenesis of epilepsy, substantial evidence implicates that immune response and inflammatory mediators contribute to the process of epileptogenesis and decrease the threshold for individual seizures.^7,8 Research confirmed that inflammatory cell infiltration and molecular regulation mechanisms in damaged neuronal tissues are the key factors for epileptogenesis. Inflammatory factors as interleukins, prostaglandins, platelet-activating factor, CD40 and CD44, matrix metalloproteinase-9, Toll-like receptors, and chemokines may create a microenvironment conducive to seizures through different mechanisms, thus increasing the susceptibility of epilepsy to seizures.⁹ Nevertheless, the underlying molecular mechanisms were still not been systematically understood, particularly regarding the effects of immunogenomic and immune cell in epilepsy. With the advent of public, large-scale high throughput genomic assays in recent years, a large number of pieces of evidence have identified biomarkers (immune-related genes) for disease surveillance and treatment with great accuracy and speed.^10,11 Analyzing related characteristics in epilepsy patients may aid the exploration of epilepsy pathogenesis and treatment.

A study by Boer et al¹² found that changes in gene expression (CCL2 and serpinA3) yield insights into new candidate genes that may contribute to seizure onset, representing new targets for potential therapeutic development. Niesen et al¹³ work firstly revealed that A deeper characterization on the dysregulated genes (C1QB and CALCRL) involved in brain tumor-induced epilepsy may shed some light on the management of epilepsy due to brain tumors. Rawat et al¹⁴ provides insight into the relevance of altered blood gene expression patterns in understanding epilepsy and its etiologic classes. Recently, Thrupp et al¹⁵ uncovered that single-nucleus RNA sequencing (snRNA-seq) is not suited for detecting cellular activation in microglia in human disease due to the low sensitivity of snRNA-seq to detect many activation genes. However, they only investigated several genes or single transcriptome data to uncover epileptogenesis. Importantly, almost all studies did not explore the correlation between hub genes, immune cells, and enriched molecular pathways (immune-related pathways). Our current study is the first to comprehensively characterize the immune landscape based on an integrated multi-transcriptome expression profile in epilepsy, in order to the discovery of novel biomarkers, thus uncovering underlying pathogenesis and prospective therapeutic targets. In the present study, we analyzed three epilepsy datasets (GSE16969,¹² GSE32534¹³ and GSE143272¹⁴) from GEO to acquire DEIRGs between epilepsy and normal samples. RF and protein–protein interaction network were constructed to screen core genes, then which was validated in another dataset (GSE31718) and 60 clinical samples via quantitative real-time polymerase chain reaction (qRT-PCR). Subsequently, CIBERSORTx and ssGSEA were performed to reflect immune cell infiltration and to identify the discrepant immune cell subtypes. GSVA and ssGSEA were conducted to determine the pathways that are significantly enriched during normal and epilepsy. Pearson correlation analysis is used to reveal the correlation between hub genes, immune cells, and enriched molecular pathways, thus shedding light on the underlying molecular mechanisms. Finally, The Connectivity Map database (CMap) was utilized to screen candidate small molecules for alleviating epileptic seizures or treating epilepsy.

Materials and Methods

Sample Selection, Data Acquisition and Processing

The GEO (http://www.ncbi.nlm.nih.gov/geo) database was used to collect the gene expression profiles of epilepsy. The inclusion criteria were the following: expression profiling by array; epilepsy; organism: Homo sapiens; samples include normal and epilepsy samples. According to the above inclusion criteria, a comprehensive search was conducted in GEO databases from inception to August 10, 2021, resulting in four GEO datasets, GSE16969,¹² GSE32534,¹³ GSE143272¹⁴ and GSE31718, being included in the analysis. The GSE16969 and GSE32534 series on the same platform GPL570 were merged as the discovery cohort, named GSEM datasets, in our study with removing the batch effect by “sva” packages, and GSE143272 series on the GPL10558 platform (Illumina HumanHT-12 V4.0 expression beadchip) was also used as discovery cohort. The GSE31718 with clinical information (pre-surgery epileptic frequency) on the GPL6254 platform (Phalanx Human OneArray) was used as validation cohort. The basic information of these microarray datasets are listed in Table 1. Additionally, ImmPort (The Immunology Database and Analysis Portal database) was utilized to select immune-related genes (IRGs) set in immunology research, containing 2489 IRGs. All data were normalized with the edgeR package or Limma package in the R computing environment.

Table 1 Basic Information of Gene Expression Profiling

Clinical Specimens

Peripheral blood mononuclear cell (PBMC) were collected from 60 clinical samples, including 30 epilepsy samples and 30 healthy controls. The study was reviewed and approved by the institutional review board (Ethics Committee) of the Shunde Hospital, Southern Medical University (the First people’s hospital of Shunde).

Evaluation of Immune Status

Microenvironment score, including immune score and stromal score, was calculated by ESTIMATE algorithm using the estimateR package on GSEM datasets and GSE16969 between normal and epilepsy tissues. The definition of immune scores is degree of immune cell infiltration. The definition of stromal scores is degree of stromal cell infiltration.

DEIRGs and Functional Enrichment Analysis

A variance analysis by NetworkAnalyst online-Gene Expression Table (https://www.networkanalyst.ca/) was performed on standardized data, GSEM datasets and GSE143272 datasets between normal and epilepsy tissues, respectively. p-values were adjusted with the Benjamini–Hochberg method. An adjusted p < 0.05 and |log2 FC (fold-change)| >1 between two groups were set as the thresholds for screening differentially expressed genes (DEGs). The overlapping DEIRGs were retained from the intersection of three data sets, DEGs of GSEM datasets, DEGs of GSE143272 datasets and IRGs set, with the UpSetR package. The location of DEIRGs on 23 chromosomes was analyzed using the RCircos package. Subsequently, DEIRGs were subject to gene ontology (GO), including the cellular component (CC), molecular function (MF), and biological process (BP) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment Analysis using R package “clusterProfiler”. And we used cnetplot R package to visualize these results.

Quantitative Real-Time PCR (qRT-PCR)

Following the manufacturer’s protocol, Trizol (Invitrogen) was used to extract total RNA from Peripheral blood mononuclear cell (PBMC). Reverse transcription of RNA using RevertAid RT Reverse Transcription Kit (Thermo Scientific). Quantitative PCR was performed using PowerUp™ SYBR™ Green Master Mix (Thermo Scientific). The results are standardized with GAPDH. Quantitative reverse transcription-PCR was conducted using the ABI 7500 real-time PCR system (Applied Biosystems, Foster City, CA, USA). Fold change was determined as 2-ΔΔCt in gene expression. Gene-specific PCR primers are listed in Supplementary Table 1.

Identification and Validation of Hub Genes

RF (Random forest), a tree-based ensemble comprised of tree-structured classifiers, was established to select feature variables via package “randomForest” with minimum error regression trees. The importance of variables was ranked using IncNodePurity. A PPI (protein–protein interaction) network of concordant DEIRGs was built by the STRING (Search Tool for the Retrieval of Interacting Genes/Proteins) online database (http://string-db.org/; version 11.5),¹⁶ and the nodes number of DEIRGs were counted. Then, we used Cytoscape software v.3.7.1 to visualize the PPI network based on STRING analysis result and screen core genes by the CytoHubba plugin.¹⁷ The real hub genes were obtained from the intersection of the result of Betweenness (top 15), Degree (top 15), MCC (top 15), PPICount (top 20), DiffGs (top 20), and RF (top 20). Then, another epilepsy dataset (GSE31718) was utilized to validate the hub genes. In addition, we compute risk score, which was generated through a linear combination of coefficients from logistic regression and the expression levels, based on hub genes in the GSE31718 database. ROC analysis with the pROC package was carried out to evaluate diagnostic performance and correlation analysis is applied to elucidate the relationship between risk score of hub genes and epileptic frequency.

Evaluation of Immune Cell Infiltration by CIBERSORTx and ssGSEA

To evaluate the relative abundance of immune infiltrates, CIBERSORTx (https://cibersort.stanford.edu/),¹⁸ which transformed the normalized gene expression matrix into the composition of infiltrating immune cells, is a kind of deconvolution algorithm with a 1000 permutation count. We filtered out samples with CIBERSORTx output of p-value > 0.05 for the accurate forecast of immune cell composition. Then, a bar graph was drawn with the package “ggplot2” to visualize the content of 22 types of infiltrating immune cells in each sample and a correlation heatmap was drawn with the “corrplot” package to visualize the correlation of immune cell subtypes. Importantly, variance analysis of immune cells between normal and epilepsy tissues was visualized by drawing violin diagrams. On the basis of the expression of metagenes that are representative of specific immune cells, the ssGSEA, using R package “GSVA”, was introduced to quantify the relative infiltration of immune cell types. We focused on the metagene set of 16 immune cell subtypes, including tumor-infiltrating lymphocyte (TIL), regulatory T cells (Treg), T helper cells, CD8+ T cells, type-2 T helper cells (Th2 cells), type-1 T helper cells (Th1 cells), natural killer cells (NK cells), macrophages, follicular helper T cells (Tfh), plasmacytoid dendritic cells (pDCs), dendritic cells, activated dendritic cells (aDCs), immature dendritic cells (iDCs), B cells, mast cells, and neutrophils, which were widely researched and accepted.¹⁹ Supplementary Table 2 will provide the metagene set. To determine differential immune cell subtypes between the two groups (p-value < 0.05), the Wilcoxon test with two-tailed was utilized to analyze the immunoscores. And we adopt vioplotR package to visualize the result.

GSVA and ssGSEA

GSVA, which converts genes from a sample matrix into predefined gene sets without a priori knowledge of experiment design, is a non-parametric unsupervised approach. The KEGG gene sets (c2.cp.kegg.v7.4.symbols.gmt), which was downloaded from the Molecular Signatures Database (MSigDB) (http://softwarebroadinstitute.org/gsea/index.jsp),²⁰ were used to estimate the variation of pathway activity in each sample. The significantly enriched pathways in KEGG gene sets were set at p-value < 0.05 and enrichment score change > 1.0. Additionally, ssGSEA, generate an enrichment score to signify the levels of absolute enrichment of a metagene set within a given dataset in each sample, was applied to evaluate the enrichment degree of immune-related pathways (cytolytic activity, T cell co-inhibition, checkpoint, T cell co-stimulation, antigen-presenting cell (APC) co-stimulation, antigen-presenting cell (APC) co-inhibition, major histocompatibility complex (MHC) class-1, type-1 IFN response, inflammation-promoting, para-inflammation, human leukocyte antigen (HLA), type-2 IFN response, and cytokine cytokine receptor interaction (CCR)¹⁹) in current immunology research. Supplementary Table 2 will provide the metagene set.

Pearson Correlation Analysis

To elucidate the potential molecular mechanism of occurrence and progression of epilepsy, Spearman correlation analyses were performed to reveal the correlation between hub genes, immune cells, and pivotal molecular pathways. A p<0.05 would be considered statistically significant.

Identification of Candidate Small Molecules Drugs

To identify potential agents for epilepsy treatment, we uploaded 51 DEIRGs, which were divided into down-modulated and up-modulated groups, into the CMap database (http://www.broadinstitute.org/cmap/). Then, the connectivity score and enrichment scores were computed, which represent similarity range from −1 to 1. A negative enrichment score (closer to −1) demonstrates that a compound can alleviate or even reverse the disease status. Whereas, a positive enrichment score (closer to +1) indicates that a compound can cause or exacerbate the disease status.

Statistical Analysis

R software (R version 3.6.1) was utilized to conduct the statistical analysis. Statistical significance was set at a two-sided p<0.05 except for where a certain p-value has been given.

Results

Evaluation of Immune Status

Whether in GSEM datasets or in GSE16969 datasets, compared with normal samples, immune score and stromal score were all higher in epilepsy samples (Figure 1), which implied that immune elements play an important role in the occurrence of epilepsy.

Figure 1 The comparison of the microenvironment score in normal samples and epilepsy samples. (A) Immune score in GESM sets. (B) Stromal score in GESM sets. (C) Immune score in GSE16969 sets. (D) Stromal score in GSE16969 sets. *p < 0.05, and **p ≤ 0.01.

DEIRGs and Functional Enrichment Analysis

After pre-processing, the distribution of expression data is consistent among all data sets. Thus, samples from each dataset are comparable. In GSEM series, we identified 1936 DEGs (Supplementary Table 3), in GSE143272 series, we identified 1410 DEGs (Supplementary Table 4). As shown in the UpSet plot, 51 DEIRGs overlapped in the three datasets (GSEM series, GSE143272 series, and IRGs set) (Figure 2A). Figure 2B displays 51 DEIRGs on the location of chromosomes was shown. The GO term enrichment analysis for DEIRGs shows the top 10 significant clusters of enriched sets (Figure 3A), such as positive regulation of chemokine production, positive regulation of T cell activation, regulation of chemokine production, chemokine production, regulation of T cell activation, and positive regulation of lymphocyte activation and so on. As for the (KEGG) biological pathway enrichment analysis of DEIRGs, they were mainly involved in T cell receptor signaling pathway, Th17 cell differentiation, Natural killer cell mediated cytotoxicity, B cell receptor signaling pathway and PD-L1 expression and PD-1 checkpoint pathway, etc (Figure 3B). The detailed results of GO and KEGG enrichment analysis were listed in Supplementary Tables 5 and 6.

Figure 2 (A) UpSet plot presents the intersection of three datasets. (B) The location of differentially expressed immune-related genes (DEIRGs) on 23 chromosomes.

Figure 3 Functional enrichment analysis of differentially expressed immune related genes (DEIRGs) in GSEM + GSE143272 datasets. (A) Gene ontology (GO) analysis on biological process (BP), cellular component (CC), and molecular function (MF). (B) Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways.

Identification and Validation of Hub Genes

To select out convincing core genes, random forest was built for feature selection and precisely calculate the importance of each feature. As shown in Figure 4A, when regression tree =52, the model error is minimum, thereby the corresponding model is constructed with sorting diagram that shows the importance of variables-DEIRGs (Figure 4B). Based on the STRING online database, we obtained the linkage between genes and calculated the number of nodes for each gene (PPICount) (Figure 5A). The PPI network was constructed by Cytoscape software and Cytoscape’s cytoHubba plugin was utilized to identify hub genes. The detailed results of cytoHubba are shown in Supplementary Table 7. According to 3 commonly used classification methods (Betweenness (top 15), Degree (top 15), MCC (top 15)) in cytoHubba, PPICount (top 20), DiffGs (top 20), and RF (top 20), which represent aspects of biological function and machine learning algorithm, we take the intersection of six methods to acquire four hub genes (CSF1R, IL6R, TLR2, and TNFRSF1A). Similar results in the GSE31718 dataset were seen by comparing the mRNA levels of the hub genes between normal and epilepsy groups (Figure S1). Additionally, the diagnostic ability of the four hub genes combined yielded an AUC of 0.92 (95% CI 0.8–1.0) (Figure 6A), indicating that combined hub genes had a high diagnostic ability. Importantly, combined hub genes were significantly positively correlated with epileptic frequency (R=0.77, p =0.0097) (Figure 6B), hinting that hub genes maybe correlated progression of epilepsy.

Figure 4 Random forest (RF) algorithm to screen hub genes in GSEM + GSE143272 datasets. (A) Distribution diagram of regression tree and error. (B) The top 30 most important variables ranked by IncNodePurity.

Figure 5 (A) The top 30 DEIRGs ranked by the number of nodes based on the PPI network. (B) UpSet plot present the intersection of six methods (Betweenness (top 15), Degree (top 15), MCC (top 15), PPICount (top 20), DiffGs (top 20), and RF (top 20)).

Figure 6 (A) The receiver operating characteristic (ROC) curves of the diagnostic efficacy from combined hub genes in GSE31718 dataset. (B) Correlation between combined hub genes and epileptic frequency. (C–F) Verification of hub genes via qRT-PCR on blood PBMCs - Comparison of gene expression levels of 4 hub genes between controls and epilepsy samples. *p < 0.05, and **p ≤ 0.01.

qRT-PCR

To further validate the expression of the 4 hub genes, we performed qRT-PCR in 60 clinical specimens. Compared with healthy controls, CSF1R, IL6R, TLR2, and TNFRSF1A was significantly up-regulated in epilepsy samples (Figure 6C–F). The detailed expression levels of 8 hub genes and clinical information are provided in Supplementary Table 8.

Immune Cell Infiltration Analysis

To investigate the differential immune landscape in normal and epilepsy tissues, microarray data from the GSE143271 dataset was subject to CIBERSORTx to estimate the fraction of 22 types of immune cells. Figure S3 displays the distribution of 22 types of immune cells in each sample and Figure S4 shown the correlation of immune cell subtypes. Compared with normal tissues, activated CD4 memory T cells (p =0.001), resting CD4 memory T cells (p =0.043), T cells gamma delta (p =0.033) and activated mast cells (p =0.025) were more abundant in epilepsy tissues (Figure 7A). This finding indicates that T cells were significantly richer in epilepsy patients. In addition, we adopt ssGSEA, another cell-type quantification method, to quantify the enrichment score of immune cell types. In GSE143271 dataset, dendritic cells (p =0.001), iDCs (p =0.001), macrophages (p<0.001), and mast cells (p =0.01) were significantly higher in epilepsy samples than in normal sample, whereas NK cells (p =0.003) and Th1 cells (p =0.004) were relatively lower (Figure 7B). In GSEM dataset, compared with normal samples, B cells (p =0.003), dendritic cells (p =0.02), neutrophils (p =0.003), pDCs (p =0.003), T helper cells (p =0.011), Th2 cells (p =0.003), TIL (p =0.003), and Treg (p =0.02) were more abundant in epilepsy samples, whereas CD8+ T cells (p =0.05), iDCs (p =0.02) and Th1 cells (p =0.006) were relatively lower (Figure 7C). Notably, ssGSEA found that Th1 cells was more abundant in normal tissues in two different datasets; yet, Th2 cells was more abundant in epilepsy tissues (one dataset shows statistical significance, another shows tendency).

Figure 7 Comparison of infiltrating immune cells between normal samples and epilepsy samples. (A) GSE143272 sets based on CIBERSORTx. (B) GSE143272 sets based on ssGSEA. (C) GESM sets based on ssGSEA.

Correlation Between Hub Genes and Immune Cells

To explore whether our hub genes were related to immune infiltration in epilepsy, Spearman correlation analyses were carried out. In GSE143271 dataset, hub genes were significantly related to B cells, CD8+ T cells, NK cells, Th1 cells, and TIL (Figure 8A). In GSEM dataset, hub genes were significantly related to B cells, neutrophils, pDCs, T helper cells, Th1 cells, Th2 cells, TIL, and Treg (Figure 8B). Based on the above results, we focused on the relational between hub genes and Th1 cell, Th2 cells. In GSE143271 dataset, IL6R (R=−0.4, p<0.001), TLR2 (R=−0.46, p<0.001), and TNFRSF1A (R=−0.44, p<0.001) were significantly negatively correlated with Th1 cell (Figure 9). In GSEM dataset, CSF1R (R=−0.87, p<0.001), IL6R (R=−0.66, p =0.016), TLR2 (R=−0.82, p =0.0011), and TNFRSF1A (R=−0.66, p =0.0016) were significantly negatively correlated with Th1 cell (Figure 9). Encouragingly, In GSE143271 dataset, CSF1R (R=−0.17, p =0.042), IL6R (R=−0.22, p =0.0097), TLR2 (R=−0.33, p<0.001), and TNFRSF1A (R=−0.31, p<0.001) were significantly negatively correlated with Th1/Th2 signatures (Figure 10). Similar results in the GSEM dataset were seen, CSF1R (R=−0.79, p =0.0023), IL6R (R=−0.88, p<0.001), TLR2 (R=−0.83, p<0.001), and TNFRSF1A (R=−0.71, p =0.0088) (Figure 10). Figure S5 visualize the relational between hub genes and Th2 cells.

Figure 8 Correlation between hub genes and infiltrating immune cells. (A) GSE143272 sets. (B) GESM sets. Yellow represents common between GESM and GSE143272. *p < 0.05, **p ≤ 0.01, and ***p ≤ 0.001.

Figure 9 Correlation between hub genes and Th1 cells in GSE143272 sets and GESM sets. (A) Association between CSF1R and Th1 cells in GSE143272 sets. (B) Association between CSF1R and Th1 cells in GESM sets. (C) Association between IL6R and Th1 cells in GSE143272 sets. (D) Association between IL6R and Th1 cells in GESM sets. (E) Association between TLR2 and Th1 cells in GSE143272 sets. (F) Association between TLR2 and Th1 cells in GESM sets. (G) Association between TNFRSF1A and Th1 cells in GSE143272 sets. (H) Association between TNFRSF1A and Th1 cells in GESM sets.

Figure 10 Correlation between hub genes and Th1/Th2 cells in GSE143272 sets and GESM sets. (A) Association between CSF1R and Th1/Th2 cells in GSE143272 sets. (B) Association between CSF1R and Th1/Th2 cells in GESM sets. (C) Association between IL6R and Th1/Th2 cells in GSE143272 sets. (D) Association between IL6R and Th1/Th2 cells in GESM sets. (E) Association between TLR2 and Th1/Th2 cells in GSE143272 sets. (F) Association between TLR2 and Th1/Th2 cells in GESM sets. (G) Association between TNFRSF1A and Th1/Th2 cells in GSE143272 sets. (H) Association between TNFRSF1A and Th1/Th2 cells in GESM sets.

GSVA and ssGSEA

To screen biological differences between normal tissues and epilepsy, GSVA was performed in GSEM databases. The results showed that CCR, apoptosis, Fc gamma R-mediated phagocytosis, gap junction, JAK−STAT signaling pathway, leukocyte transendothelial migration, and Toll−Like receptor signaling pathway were enriched in the epileptic samples (Figure 11). In addition, ssGSEA was utilized to explore the predefine immune-related pathway in epileptogenesis. In GSEM databases, the ssGSEA results showed that CCR, checkpoint, MHC class I, T cell co-inhibition, T cell co-stimulation, and type II IFN response were involved in the mechanism of epilepsy (Figure 12A). In GSE143271, the ssGSEA results showed that APC co-inhibition, APC co-stimulation, CCR, para-inflammation, and type I IFN were mainly enriched in epilepsy (Figure 12B).

Figure 11 Heatmap displays differences in pathway activities between normal samples and epilepsy samples in GESM datasets. Red indicates epilepsy samples, while blue indicates normal samples.

Figure 12 Comparison of the immune-related pathway between normal samples and epilepsy samples. (A) GESM sets. (B) GSE143272 sets.

Correlation Between Hub Genes and Pivotal Molecular Pathways

To determine the molecular pathways potentially involved in epileptogenesis, Pearson correlation analysis was performed to uncover the correlation between hub genes and enriched molecular pathways. In GSEM databases, based on enrichment score of GSVA, hub genes were significantly related to CCR, apoptosis, Fc gamma R-mediated phagocytosis, focal adhesion, JAK−STAT signaling pathway, and leukocyte transendothelial migration (Figure 13A), according to enrichment score of ssGSEA, hub genes were significantly related to CCR, checkpoint, MHC class I, T cell co-inhibition, and type II IFN response (Figure 13B). In GSE143271 databases, based on the enrichment score of ssGSEA, hub genes were significantly correlated with CCR, checkpoint, APC co-inhibition, APC co-stimulation, para-inflammation, T cell co-stimulation, T cell co-inhibition, and type I IFN response (Figure 13C).

Figure 13 Correlation between hub genes and pivotal molecular pathways. (A) GESM sets based on GSVA. (B) GESM sets based on ssGSEA. (C) GSE143272 sets based on ssGSEA. *p < 0.05, **p ≤ 0.01, and ***p ≤ 0.001.

Identification of Related Active Small Molecules

10 most significant small molecules agents mitigating epileptic seizures are listed in Table 2. Among these small molecules drugs, wortmannin (enrichment score = −0.981), securinine (enrichment score = −0.897) and parthenolide (enrichment score = −0.867) were correlated with a highly significant negative fraction and possess the potential to alleviate or even reverse epileptic seizures.

Table 2 List of the 10 Most Significant Small Molecule Drugs That Can Reverse Epilepsy

Discussion

Integrating multi-transcriptome datasets, we identified 51 DEIRGs. Based on RF and PPI network, 4 DEIRGs (CSF1R, IL6R, TLR2, and TNFRSF1A) were identified as hub genes. Results of qRT-PCR validated that higher expression levels of CSF1R, IL6R, TLR2, and TNFRSF1A in epilepsy compared to control sample. Immune infiltration analysis by CIBERSORTx displayed immune cells that are significantly richer in epilepsy, including activated CD4 memory T cells, resting CD4 memory T cells, T cells gamma delta and activated mast cells. Noteworthily, ssGSEA found that Th1 signatures was more abundant in normal tissues; yet, Th2 signatures was signatures was more abundant in epilepsy tissues. According to GSVA and ssGSEA, in spite of combined GEO or GSE143272, CCR were all significantly enriched in epilepsy. Additionally, hub genes were significantly correlated with infiltration score of Th1/Th2 signatures and enrichment score of CCR in multi-transcriptome data. Finally, CMap matched the 10 most important candidate molecules as potential therapeutic drugs to reverse epilepsy.

An increasing body of evidence indicates that several components of the immune system, including alteration of the function or number of immune cells, the release of pro-inflammatory cytokines and chemokines, and the expression of inhibitory receptors or their ligands, are key players in epilepsy occurrence and progression.^21–23 A study by Tang et al found that PD-1 levels in cerebrospinal fluid (CSF) and serum were more pronounced in subjects with intractable status epilepticus, which can serve as a potential clinical diagnostic biomarker for intractable epilepsy.²¹ Gnatek et al work firstly revealed that acetylcholine (ACh) acts directly as a modulator of the brain to suppresses brain innate immune response, and that upregulation of acetylcholinesterase (AChE) following status epilepticus (SE) is related to enhanced immune response, facilitating the epileptogenic process.²²

Gales et al observed that there was a mixture of both B (CD20 +) and T (CD3+) lymphocytes distributed around the blood vessels and dispersed in the brain parenchyma in hippocampal sclerosis-related temporal lobe epilepsy.²³ Interestingly, CD4 positive T cells were more dominant than CD8 positive cells. The phenomenon exerts neurotoxicity and cause neuroglial injury. These studies support the involvement of immune and inflammatory processes in epilepsy. However, the underlying molecular mechanisms remain not been systematically understood, especially with regard to the immunogenomic and immune infiltration effects of epilepsy. A meta-analysis of microarray data sets from animal studies implicating immune system-related genes in all phases of epileptogenesis.²⁴ Hence, it is reasonable to characterize immune-related genes and immune infiltration characteristics in epilepsy to aid the exploration of epileptogenesis.

To our knowledge, this is the first comprehensive study to characterize the immune landscape based on an integrated multi-transcriptome expression profile in epilepsy, leading to the discovery of novel biomarkers, thus uncovering underlying pathogenesis and prospective therapeutic targets. Before diving into immunology etiology, we explore the differential immune score in normal and epilepsy tissues.

Compared to normal subjects, immune scores and stromal scores were all higher in epilepsy samples in different datasets. It was a hint that immune elements were involved in the occurrence of epilepsy. On the basis of DERs of multiple transcriptome data and IRGs set, we acquired 51 DEIRGs. GO term and KEGG enrichment analysis indicate that T cell and chemokine production are the most likely immune components involved in pathological processes and immunoregulatory mechanisms of epilepsy. To select out convincing key genes, RF, represent machine learning algorithm, and PPI network (cytoHubba), signify biological function, were incorporated to identify four hub genes (CSF1R, IL6R, TLR2, and TNFRSF1A). Importantly, Results of qRT-PCR validated that expression level of CSF1R, IL6R, TLR2, and TNFRSF1A were significantly up-regulated in epilepsy compared with control sample in clinical specimens, which in line with the results of bioinformatics analyses derived from the GEO datasets. Interestingly, four hub genes combined yielded a superior diagnostic capacity and were closely related to epileptic frequency. It implied that hub genes may be involved in the initiation, maintenance, and progression of epilepsy.

4 IRGs (CSF1R, IL6R, TLR2, and TNFRSF1A), are belongs to cytokines receptor, which are mediators of the immune system and via binding to corresponding ligands coordinate activities of the acquired and innate immunity. CSF1R, is a class III receptor tyrosine kinase family, which plays a key role in signal transduction through binding ligands CSF-1 and interleukin (IL)-34 to modulate monocyte/macrophage lineage cells.²⁵ CSF-1/CSF-1R axis overactivation may result in aberrant expression of proinflammatory cytokines including interleukin-6 (IL-6) and tumor necrosis factor-α (TNF-α), among others, which can promote the pathogenesis of inflammatory diseases. There is recent evidence that the elimination of microglia by CSF-1R inhibition could curtail plaque formation in an in vivo model of Alzheimer’s disease.²⁶ IL6R, which is an inflammatory cytokine receptor, is only expressed on specific cells including some leukocytes, some epithelial cells and hepatocytes, which are therefore IL-6 target cells.²⁷ The complex of IL-6R and IL-6 bind to a second receptor subunit, gp130, which thereupon dimerizes resulting in the activation of SHP2/Ras/MAPK, YAP, JAK/STAT signaling pathways. The process has been named IL-6 trans-signaling, whose disruption plays a prominent role in inflammation diseases. Currently, IL-6R inhibition, tocilizumab, is approved for the treatment of juvenile idiopathic arthritis, giant cell arteritis and rheumatoid arthritis. Moreover, recent research reveals that an open-label, randomised, multicentre, phase two trial, tocilizumab significantly reduced the risk of neuromyelitis optica spectrum disorder (NMOSD) relapse, and is a safe and effective treatment for NMOSD.²⁸ TLR2, as a member of the Toll-like receptors (TLRs) family, binds the cellular components of the host, thus participated in innate immunity and serves as a bridge between innate and acquired immunity.²⁹ It transmits extracellular signals to activate transcription factors, such as AP-1 and nuclear factor-κB, which contributes to the production of proinflammatory cytokines and chemotactic cytokines, thus regulating the balance between Th1 cells and Th2 cells. Research finding that TLR 2, is activated and up-modulated by high mobility group box-1 (HMGB1), results in surgery-induced neuro-inflammation and impairment of memory and learning.³⁰ TNFRSF1A, is a tumor necrosis factor receptor superfamily member, act as one of the primary transmembrane signal transduction receptors for tumor necrosis factor-α (TNF-α) that is a specific pro-inflammatory cytokine, and plays a critical role in cell fates, induction of cytokines and chemokines.³¹ Its binding with TNF-α, named TNF-α signaling pathway, can activate MAPK/NF-κB signaling pathways and involved in inflammatory diseases. It has been reported that stimulation of cerebellar fastigial nucleus can down-modulate the expression of TNFRSF1A in the rat model, thus slightly reducing neuronal apoptosis.³²

Immune infiltration analysis via CIBERSORTx found immune cells that are more abundant in epilepsy, such as activated CD4 memory T cells, resting CD4 memory T cells, and T cells gamma delta. Recently, Ouédraogo et al observed that compared to healthy controls, CD4 T-cell populations were more abundant in the peripheral blood of focal onset epilepsy subjects, in line with our results.³³ Author data support that epilepsy is linked to an expansion of the CD4 T cell subset, which is skewed toward a proinflammatory Th17/Th1 phenotype. Additionally, Xu et al for the first time demonstrate that proinflammatory IL-17-producing T cells gamma delta are concentrated in the epileptogenic zone, and their populations positively related to seizure severity,³⁴ which are corroborated by animal experiments in the brains of a mouse model of SE induced by the kainic acid displaying similar activation of adaptive and innate immune responses, in accordance with our analysis data. Remarkably, ssGSEA found that Th1 signatures was more abundant in normal subjects in different datasets; yet, Th2 signatures was signatures was more abundant in epilepsy subjects. Hence, a hypothesis is deduced that the imbalance of Th1/Th2 signatures is associated with epileptogenesis. Encouragingly, hub genes were significantly correlated with infiltration score of Th1/Th2 signatures, which indirectly support that Th1/Th2 imbalance plays a crucial role in the pathogenesis of epilepsy. Under normal circumstances, there is cross-regulation between Th1 and Th2 cells by cytokines, that is, Th1 cells via secreting IL2 or FIN inhibit Th2 cells, while Th2 cells through secreting IL-4, IL-6, or IL-10 inhibit Th1, which plays a role in maintaining the normal immune function. Whereas, when the function is abnormal, Th1/Th2 balance often appears a bias toward one side, which is called “Th1/Th2 shift”, that may disrupt the dynamic balance of human cytokines network and contribute to the occurrence and progression of diseases, such as rheumatoid arthritis, tumors, asthma and recurrent respiratory infections.³⁵ In our current immunology research of epilepsy, Th1/Th2 balance may shift from Th1 to Th2, namely Th1-type signatures relatively decreased while Th2-type signatures relatively increased, which hints that humoral immune maybe plays a dominant role in the etiology of epilepsy. To sum up, Immune cell infiltration, T cell subsets, in particular, may be compactly linked to the pathogenesis of epilepsy and need further experiments to prove those phenomena. Additionally, regulators of T cell subsets,³⁶ or microglia/macrophages with biologic agents after epilepsy has become a focus of attention in recent years, and the discovery of new drugs related to T cell subsets or microglia/macrophages polarization may make the realization of immunotherapy for epilepsy possible.

GSVA results showed that CCR, apoptosis, chemokine signaling pathway, JAK−STAT signaling pathway, Toll−Like receptor signaling pathway, and so on were involved in the epileptic subjects, in accord with the results of GO term and KEGG enrichment analysis. They are common and classic signaling pathways that modulate immune responses and inflammatory mediators, and their overactivation leads to the onset and development of disease. ssGSEA results demonstrate CCR was significantly involved in epilepsy and enrichment score of it was significantly correlated with hub genes in multi-transcriptome data, which indirectly confirm that Th1/Th2 shift is likely to contribute to the occurrence of epilepsy. The research found that chemokines (CCL3 and CCL4) were highly up-modulated in chronic temporal lobe epilepsy (TLE), which was observed in the early stages of epileptogenesis in animal models. The underlying mechanisms may be directly or directly affect neuronal excitability to prolonged seizures and enhance immune response, facilitating the epileptogenic process.³⁷ There is growing evidence that inflammatory processes in the brain may be critical for both the development of epileptic foci and the continuation of epileptic relapse.^38,39 Ion channels and Neurotransmitter receptors are regulated by immune mediators,⁴⁰ and during seizures, autonomic discharges alter the distribution of circulatory peripheral immune cell populations.⁴¹ Seizures destabilize the (BBB) (blood–brain barrier) and activate microglia, which in turn promote elevated levels of pro-inflammatory cytokines, IL-1 and TNF in the central nervous system.^9,42 In addition, changes in peripheral blood luminal TNF, IL-6, and IL-1 levels have been reported in epilepsy and are related to recent seizure activity.⁹ Taken together, CCR is probably one of the main immune signaling pathways that contributes to the onset of epilepsy and need further experiments to explore the specific mechanism.

In spite of the remarkable sense, it is inevitable that limitations also existed in our research. First, study patients varied widely from GEO data, including patients with tuberous sclerosis complex (TSC), patients with primary brain tumors, patients with refractory epilepsy. However, our results (hub genes, immune infiltration cells and pivotal molecular pathways) take the intersection of different GEO databases. Additional high-quality studies with single type of epilepsy need to validate our results. Second, though we integrated different epileptic microarray datasets, the sample size remained limited, in part because not all epileptic patients underwent brain surgery or peripheral blood whole blood sequencing. Large sample size can improve the statistical power to detect robust immune-related genes and immune infiltration characteristics in epileptogenesis. Third, patients with epilepsy included in our analysis were not guaranteed to be free of other diseases. Whereas, the influence of other diseases on our results cannot be fully resolved because the original data set did not offer complete details of other comorbidities/diseases. Fourth, Based on bulk RNA-Seq data, CIBERSORTx deconvolution algorithm and ssGSEA with metagenes may not accurately identify immune cell subpopulations although different methods and different data sets validate each other. It is necessary to use single-cell RNA-Seq methods or fluorescence-activated cell sorting to verify our results.

Fifth, no further in vivo experiments to validate these results (hub genes, immune infiltration cells and pivotal molecular pathways). Loss of function and overexpression studies in vitro, as well as in animal models, will help to further identify the exact role of hub genes in the regulation of the inflammatory response and related pathogenic signaling in epilepsy.

Conclusion

4 IRGs (CSF1R, IL6R, TLR2, and TNFRSF1A) were closely associated with the pathogenesis of epilepsy, which maybe by afflicting Th1/Th2 signatures balance and CCR pathway involved in the onset of epilepsy. Our data may provide a new basis for understanding the role of the immune landscape in the treatment of epilepsy.