Key protein-coding genes related to microglia in immune regulation and inflammatory response induced by epilepsy

1.   Introduction

Epilepsy is a chronic neurological disorder that directly affects an estimated 68 million people worldwide [1]. Current clinical treatments manly involve the use of antiepileptic drugs and/or surgical resection of the seizure focus, whereas, they are associated with multiple adverse effects, and few patients are seizure-free after surgery. Thus, the intractability of epilepsy has prompted researchers to maintain a sustained interest in its biological mechanisms and the development of new clinical therapeutic targets.

The different roles of inflammatory responses and immune signals in the development of epilepsy have been highlighted in several studies. Brain injuries associated with central nervous system infections may lead to early occurrence of seizures, which may result in epilepsy [2]. Concurrently, seizures trigger inflammatory responses. Several inflammatory factors, including IL-1β, IL-8, IL-12p70 and MIP-1β, have been detected in surgically resected epileptogenic cortex obtained from children with refractory epilepsy [3]. Animal studies have also found increased levels of pro-inflammatory cytokines in animal models of epilepsy and that seizure frequency correlates with the activation of these cytokines [4–6]. Thus, a number of small-molecule compounds targeting inflammatory signaling pathways, especially for COX-2/PGE₂/EP2 axis, have been extensively studied in different animal models of epilepsy. Treatment with EP2 antagonists after seizures showed overall significant anti-inflammatory and neuroprotective effects in three rodent models of status epilepticus (SE) induced by pilocarpine, kainic acid, and diisopropyl fluorophosphate, respectively, resulting in improved behavioral and functional deficits caused by prolonged seizures [7]. Alternatively, Studies on the use of various COX-2 inhibitors in SE animal models have produced controversial results duo to they inhibit some beneficial effects that might be mediated by other COX-2-derived prostanoids, suggesting a possible consequence of COX-2 induction - early neuroprotection and later neurotoxicity [8].

Lately, growing studies have shown that cytokine release is profoundly dependent on microglia activation, which is prevalent in human epilepsy and animal models, though their phenotypes have exhibited heterogeneity. Boer et al. studied the distribution of microglia in 20 specimens of focal cortical dysplasia(a major cause of intractable epilepsy) by immunocytochemistry and found a specific and sustained increase in the numerical density of HLA-DR-positive activated microglia in areas of abnormal development, significantly higher than in normal-appearing control cortex obtained at autopsy [9]. Benson et al. detected the changes of microglial pro-inflammatory M1 and anti-inflammatory M2 marker expression in different stages of different mouse models. M1 and M2 markers showed increased expression 3 days but not 21 days after pilocarpine-induced status epilepticus. Though, In the late chronic phase, some M1/M2 markers, IL-1β, TNFα, Arg1, Ym1, and CD206, resurged in the kainate, but not pilocarpine model, indicating a complex role of microglia in the epileptic brain [10]. Accordingly, discovering new biomarkers of microglia activation is essential for further studies on its relationship with neuroinflammation and epileptogenes.

Given the popularity of high-throughput microarray and RNA-sequencing (RNA-Seq) analyses as research tools, the use of bioinformatics methods to improve the understanding of transcriptome changes in the onset and development of epilepsy has increased. Kalozoumi et al. reported that glial-mediated inflammatory responses play a key role during epileptogenesis in a model of acute epilepsy induced by kainic acid (KA), with hippocampal genome-wide expression analysis at 6, 12 and 24h post-injection [11]. Luo et al. explored the expression and correlation of long non-coding RNAs (lncRNAs) and microRNAs (miRNAs) in a rat model of temporal lobe epilepsy (TLE). Results showed that 54 lncRNA, 36 miRNA, and 122 mRNA molecules that were dysregulated in experimental groups compared to normal controls [12]. Recently, Chen et al. performed a comprehensive analysis of six TLE microarray datasets using robust rank aggregation to identify co-expressed key genes across different epilepsy models and platforms. These authors have shown that the differentially-expressed genes (DEGs) in the acute, latent, and chronic phases of TLE were significantly enriched in inflammatory response, while the pivotal genes identified in each stage were associated with different biological processes [13]. However, these results were based on data from sequencing analysis of animal model samples, which do not correspond completely to the changes in transcriptome signature during human epileptogenesis. Identification of critical protein-protein interaction (PPI) modules and key regulatory genes within them is crucial for elucidating the pathogenesis of epilepsy.

In the present study, we collected brain tissues from two neurologically healthy controls and five epileptic patients, screened the DEGs in both groups, and applied bioinformatics methods to identify the most important biological process and hub genes associated with microglia. In addition, real-time quantitative polymerase chain reaction (RT-qPCR) was used to validate the expression pattern of several hub genes in a mouse model of epilepsy.

2.   Materials and methods

2.1.   Sample collection

Brain tissue was isolated within 12 hours after death from two subjects without a history or signs of neurological or neuropsychiatric illness, and provided by the Department of Forensic Medicine, Chongqing Medical University. Surgical specimens were obtained from five secondary epileptic patients caused by brain maldevelopment who had undergone surgical resection at the Department of Neurosurgery, Children's Hospital of Chongqing Medical University. Details related to sample collection, transportation and preservation are described in the Supplementary Materials. Clinical characteristics of the patients and controls are summarized in Table 1. For more information, see Supplementary Materials Table S1. The Ethics Committee of the Children's Hospital Affiliated to Chongqing Medical University approved the present study protocol involving the use of human tissue. Written informed consent was signed by the participants or their legal representatives prior to the start of the project.

Table 1.  Characteristics of controls and epileptic patients.

Identifier	Gender	Age (years)	Cause of death
C1	Female	71	Renal failure
C2	Male	38	Car accident
Identifier	Gender	Age (years)	Surgical approach
E1	Male	4	Resection of frontal lobe lesions
E2	Male	9	Resection of frontal lobe lesions
E3	Male	5	Resection of parietal lobe lesions
E4	Male	2	Resection of temporal lobe epilepsy lesions
E5	Male	11	Resection of temporal lobe epilepsy lesions

---

 | Show Table

DownLoad: CSV

Adult male C57BL/6J mice were provided by the Experimental Animal Center of Chongqing Medical University and housed in specific pathogen free (SPF) environment. The KA-induced chronic epileptic model was established, as previously described [14]. Briefly, deeply anesthetized mice were injected with 1.0 nmol of KA (Sigma Aldrich, St. Louis, USA) dissolved in 50 nL saline into the right dorsal hippocampus (anteroposterior = -1.8 mm, mediolateral = -1.8 mm, dorsoventral = -1.9 mm) with bregma as reference. Control mice were injected 50 nL saline in the same manner. At least one behavioral spontaneous recurrent motor seizure observed during 4 weeks post-injection was considered a successful model of chronic epilepsy. All mice were sacrificed 4 weeks after surgery. All experimental procedures involving animals were reviewed and approved by the Laboratory Animal Ethics Committee of Chongqing Medical University.

2.2.   RNA sequencing

Total RNA was extracted using TRIzol Reagent (Invitrogene, USA), according to the manufacturer's instructions. Subsequent library construction and Illumina sequencing were performed by Shanghai Lifegenes Technology Co., Ltd. All downstream analyses were based on clean data, which were obtained from raw data by removing reads containing adapters, reads containing ploy-N and low-quality reads. FPKM, Fragments Per Kilobase of exon model per Million mapped fragments, is commonly used to estimate the abundance of gene expression [15]. FPKMs of the mRNAs in each sample were calculated by aligning paired-end clean reads to the reference genome utilizing HISAT2 v2.1.0., and the FPKM for each gene was computed by summing the FPKMs of transcripts in each genome. Genes with FPKM ≥ 1 were selected for subsequent analysis to remove low expression genes.

2.3.   Bioinformatics analysis

The DESeq R package (1.8.3) was used to analyze the differentially expressed mRNAs between the epilepsy and control groups. Adjusted Р-values < 0.05 and |log2(foldchange)| ≥ 1 were set as the threshold forsignificantly differential expression.

To elucidate interactions among the differentially expressed genes, a PPI network was constructed using STRING tool (https://www.string-db.org/) and visualized using Cytoscape software(http://cytoscape.org/) to reveal the molecular mechanisms underlying epilepsy. The Molecular Complex Detection (MCODE; Version:1.4.2) plugin in Cytoscapse software was used to screen out the key modules [16]. Parameters were set with a degree cutoff of 2, node score cutoff of 0.2, K-core of 4, and a maximum depth of 100.

Gene Ontology (GO) and KEGG enrichment analyses were performed using the online tool DAVID 6.8 (https://david.ncifcrf.gov/). Genes in the key modules were analyzed using the following parameters: OFFICIAL-GENE-SYMBOL, species: Homo Sapiens. GOTERM_BP, GOTERM_CC, GOTERM_MF and KEGG_PATHWAY datasets were downloaded. BP, CC, and MF refer to biological process, cellular component, and molecular function, respectively. GO and KEGG terms with P-values < 0.05 were considered significantly enriched. The hub genes were distinguished by cytoHubba in Cytoscape [17]. The MCC algorithm was used to sort the genes in the network, and top 10 nodes were selected as hub genes, which were further annotated by ClueGO and CluePedia App [18,19].

2.4.   RT-qPCR

Total RNA was extracted from the right hippocampus tissue of mice in the control and epilepsy group induced by KA. Ten hub genes in key modules from RNA-seq were validated by RT-qPCR using the CFX96 Real-Time System (Bio-Rad, USA). Glyceraldehyde-3-phosphate dehydrogenase (GAPDH) served as an internal control. The primer sequences were designed using Primer Premier 5.0, and are listed in Table 4. Experiments were performed in triplicates. Relative gene expression was analyzed using the 2^−ΔΔCt method.

2.5.   Statistical analysis

Statistical analyses were performed using GraphPad Prism 7.0 (GraphPad Software, La Jolla, CA, USA). Two-tailed unpaired t-test was used to the compare two groups. Differences were considered statistically significant at P-values < 0.05.

3.   Results

3.1.   DEGs between controls and epileptic group

The results of the DEG analysis are presented in volcano plot and heatmap (Figure 1). As shown in Figure 1a, a total of 1426 DEGs were detected, including 975 upregulated genes and 451 downregulated genes. Clustering analysis of genes with the same or similar expression patterns may help identify the functions of unknown genes or unknown functions of known genes. Thus, we conducted a hierarchical clustering analysis on the top 100 differentially expressed mRNAs based on the FPKM values in the two groups (Figure 1b). The different color areas represent different clusters of grouping information. Genes in the same group have similar expression patterns and may have similar functions or be involved in the same biological processes. In addition, we summarized ten genes with the most significant differences in Table 2.

Figure 1.  The volcano plot and heatmap of differentially expressed genes. (a) Volcano plot of the DEGs. (b) Heat map of top 100 DEGs in the epileptic group compared with the control group.

DownLoad: Full-Size Img PowerPoint

Table 2.  The top ten DEGs.

Gene name	Description	log FC	adj. P-Value	Up/Down
FOSL2	FOS like 2, AP-1 transcription factor subunit	2.930792	7.05E-56	Up
SERPINF1	serpin family F member 1	4.146166	2.81E-33	Up
OLFML3	olfactomedin like 3	3.385027	2.91E-24	Up
TSPYL2	TSPY-like 2	2.173859	3.86E-24	Up
CCDC88A	coiled-coil domain containing 88A	-1.887247	1.18E-23	Down
TRABD2A	TraB domain containing 2A	2.949808	1.05E-21	Up
ALMS1	ALMS1, centrosome and basal body associated protein	-1.339524	4.71E-20	Down
PACSIN1	protein kinase C and casein kinase substrate in neurons 1	2.853895	4.01E-19	Up
STIM1	stromal interaction molecule 1	1.857381	1.12E-18	Up
C1QA	complement component 1, q subcomponent, A chain	4.279842	1.12E-18	Up

---

 | Show Table

DownLoad: CSV

3.2.   PPI network construction and key modules identification

To further explore the crosstalk between the DEGs, top 500 DEGs were selected to construct a PPI network based on the STRING database. The network consists of 389 nodes and 1480 edges (Figure 2a). The size of the node represents the degree of connectivity, and the color of the node represents the gene expression pattern. The red nodes are upregulated genes, while blue nodes are downregulated genes. Multiple upregulated genes were at the center of the network. A small part of the PPI network with highly connected areas has a higher probability of participating in biological regulation, and the hub genes among the subnets are often of clinical importance [20]. To find biologically essential subnets and corresponding hub genes, the MCODE [16] plugin was used to investigate the whole network. Four subnets were eventually derived from the source network (Figure 2b–e). Table 3 shows the detailed information of each subnet. As the subnet with the highest score, cluster 1, which contains 34 nodes and 287 edges, was selected as the key module for further research.

Figure 2.  The protein-protein interaction network. (a) PPI network for the top 500 DEGs. (b) subnet 1. (c) subnet 2. (d) subnet 3. (e) subnet 4. Red nodes represent upregulated genes. Blue nodes represent downregulated genes.

DownLoad: Full-Size Img PowerPoint

Table 3.  Details of four subnets.

Cluster	Score	Nodes	Edges	Node IDs
1	17.394	34	287	LAPTM5, CSF1R, CD53, IL10RA, SPI1, AIF1, GPSM3, ADRA2C, PDYN, SSTR3, FPR1, GPR183, TAS2R20, CX3CR1, HRH3, LPAR5, P2RY13, P2RY12, CSF2RB, MCHR1, NCF2, CX3CL1, ADCY2, CYBB, ITGB2, LY86, TLR7, PLEK, HCK, C1QA, C1QC, FCGR1A, CTSS, CD86
2	5	5	10	CAMK2A, SLC6A7, KCNA1, GAD2, GABRD
3	4.5	13	27	TACR3, C3AR1, KCNAB2, HTR2C, SNAP25, FOLR2, GNA15, ATP6V1D, CCK, HTR2A, F13A1, VSIG4, HRH1
4	4	8	14	MERTK, JUP, TMEM119, HPSE, NEU1, TOLLIP, CEP290, GPR34

---

 | Show Table

DownLoad: CSV

3.3.   Functional enrichment analysis of key module

To clarify the function of genes in key module, GO annotation and KEGG pathway analyses were performed with DAVID 6.8 and visualized by bioinformatics online tool (http://www.bioinformatics.com.cn/). Figure 2 illustrates the top 10 significantly enriched GO and KEGG terms (P < 0.05). Items related to biological process were mainly enriched in inflammatory response, innate immune response and G-protein coupled receptor signaling pathway (Figure 2a). Meanwhile, the top three enriched terms of cellular components were plasma membrane, integral component of plasma membrane, integral component of membrane (Figure 2b). With regard to molecular functions, the most important functions were concentrated in receptor activity, superoxide-generating NADPH oxidase activity, G-protein coupled purinergic nucleotide receptor activity (Figure 2c). The most significantly enriched pathways were associated with Staphylococcus aureus infection, Osteoclast differentiation, Neuroactive ligand-receptor interaction (Figure 2d). Overall, these findings suggest that the immunoregulatory and inflammatory responses occupy a dominant position in the pathogenesis of epilepsy, which is consistent with previous studies [21].

Figure 3.  Functional enrichment analysis of DEGs in cluster 1 based on GO and KEGG. (a) BP: biological process; (b) CC: cellular component; (c) MF: molecular function; (d) KEGG: signaling pathway.

DownLoad: Full-Size Img PowerPoint

3.4.   Screening hub genes in the key subnet

We obtained the hub genes in the module related to immunoregulatory and inflammatory responses using the cytoHubba plugin in Cytoscape. The MCC algorithm was used to sort the genes in the network, and top 10 nodes were selected, including CX3CR1, CX3CL1, GPR183, FPR1, P2RY13, P2RY12, LPAR5, ADRA2C, PDYN, SSTR3 (Figure 4a). A darker color indicates a higher score, signaling a potentially more important role in the network. In addition, we explored the GO term associated with immune system process of these genes using ClueGO and CluePedia application in Cytoscape, revealing that they were all related to regulation of microglial cell migration and macrophage migration (Figure 4b).

Figure 4.  The hub genes. (a) Protein-protein interaction network consisting of 10 hub genes. The shade of the node color indicates its importance in the network. (b) GO terms associated with immune system process. Each node represents a term, and the connection between the nodes reflects the correlation between the terms.

DownLoad: Full-Size Img PowerPoint

3.5.   Validation by real-time quantitative PCR

To verify the results of RNA-seq analysis, 10 hub genes were chosen for analysis using RT-qPCR. All genes had homologous genes in the mice. Gene name, description, expression of sequencing results, and primer information are shown in Table 4. Figure 5 shows the expression of each gene in epileptic mice. The expression of nine genes was consistent with the sequencing results, except for FPR1, which is down-regulated in epileptic mice compared with the control group, in contrast to the outcome of RNA-seq. This finding may be due to the inconsistent age in human specimens for sequencing and animal samples for verification, or due to species differences between humans and mice.

Table 4.  The top10 hub genes and related information obtained from sequencing results.

Gene name	Description	Up/Down	Primer
CX3CR1	C-X3-C motif chemokine receptor 1	Up	F:TTCCCATCTGCTCAGGACCTC R:CAGACCGAACGTGAAGACGA
CX3CL1	C-X3-C motif chemokine ligand 1	Up	F:ACGAAATGCGAAATCATGTGC R:CTGTGTCGTCTCCAGGACAA
GPR183	G protein-coupled receptor 183	Up	F:GGATGACAGGAAAAAGGCACC R:CCAGTGGTCACAGTCAATATCAG
FPR1	formyl peptide receptor 1	Up	F:GCTGTTGGAAAGTTCAGGAGTC R:CCAGAACGATGTAGCCAGCA
P2RY13	purinergic receptor P2Y13	Up	F:TGGGTTGAGCTAGTAACTGCC R:TTGTCCCGAGCATCAGCTTT
P2RY12	purinergic receptor P2Y12	Up	F:AGGGTCACAGTGCAAGAACA R:TGGAACTTGCAGACTGGCAT
LPAR5	lysophosphatidic acid receptor 5	Up	F:CAGAGGTCCCAAGAAGGTCTCC R:CAGAAGCAAGACATGGGCAC
ADRA2C	adrenoceptor alpha 2C	Up	F:TGCACCTGTGTGCCATTAGT R:GTAGAACGAGACGAGAGGCG
PDYN	prodynorphin	Up	F:CTCCAGAACTGCCATAGGGG R:TGGGGATGAATGACCTGCTTAC
SSTR3	somatostatin receptor 3	Up	F:CCAAGTGCCAGAAAATGGCTG R:AACCAGCCCTAACCTGAAGGG

---

 | Show Table

DownLoad: CSV

Figure 5.  Validation of 10 hub genes by qPCR. Black bars represent gene expression in the control group and gray bars represent gene expression in the epilepsy group. * represents P < 0.05; ** represents P < 0.01; *** represents P < 0.001; **** represents P < 0.0001.

DownLoad: Full-Size Img PowerPoint

4.   Discussion

In the present study, we sequenced the transcriptomes of seven samples derived from human brain tissue, screened differentially expressed genes and performed a series of bioinformatics analyses. We found that genes in the key networks were significantly enriched for immune and inflammatory processes. Ten hub genes were then selected from the key module and a gene ontology annotation analysis was performed based on immune system process. The present findings have shown that these genes are most closely related to the migration of microglia and macrophages, which may identify new immunomodulatory compounds in the development and progression of epilepsy.

An updated review article has examined the diversity of microglia functions, including multiple aspects of microglial activation in the hippocampal sclerotic area, microglial role in epileptogenesis without inflammatory activation, anti-epileptic role of microglia, microglia-neuron interaction and perspectives on glial communication in the epileptic brain [22]. Several hub genes have been associated with microglia. CX3CL1 is a member of the chemokine gene family, also known as fractalkine, which is expressed abundantly in the brain and localized principally to neurons. Numerous studies have reported that CX3CL1 binds to its unique receptor CX3CR1, regulates the activation of microglia, and appears to function as a signal molecule from neurons to microglia [23–27]. A few studies have shown the role of CX3CL1/CX3CR1 axis signaling in the pathogenesis of epilepsy. Xu et al. reported increased expression of CX3CL1 in patients with TLE and rat models [28]. Roseti et al. have shown that CX3CL1 is responsible for positively regulating the function of the GABA A receptor in human TLE brain tissue when it was expressed in Xenopus oocytes, indicating that the CX3CL1/CX3CR1 axis is a critical regulator of potentially neuroprotective microglia-neuron interactions [29]. Additionally, Eyo and colleagues have found that a deficiency of the CX3CL1 receptor reduced the microglial process convergence toward neuronal axons and dendrites, and blocking fractalkine signaling was associated with increased seizure phenotypes [30].

GPR183 and FPR1are functionally similar in that they are both receptors expressed in the central nervous system, and their endogenous ligands direct cell migration. Rutkowska and colleagues found that GPR183 co-expressed with 7α, 25-OHC, which was considered the primary endogenous ligand for GPR183 [31], in human and murine astrocytes [32]. Wanke et al. reported that EBI2 is highly expressed in multiple sclerosis (MS) lesions and promotes early central nervous system migration of encephalitogenic CD4 T cells [33]. A recent study revealed high expression of EBI2 in astrocytes and microglia in MS brain plaques, indicating that EBI2 signaling is involved in the recovery of demyelination [34]. FPR1 encodes a high affinity receptor for N-formyl-methionyl peptides, which are powerful neutrophil chemotactic factors [35–37]. FPR1 has also been reported to induce glial cell activation in a mouse model of MS [38] and to promote the migration and differentiation of rat neural stem cells [39]. However, these two genes have not been reported in epilepsy.

P2RY12 and P2RY13 are members of the metabotropic purinergic P2Y receptor family. P2RY12 was confirmed by Sasaki et al. to be selectively expressed on microglia [40], and has been considered as a marker to distinguish microglia from other brain cells and peripheral immune cells [41]. Evidences shows it may be necessary for microglial translocation [42], response and migration to central nervous system injury [43], and communication with neurons after status epilepticus [44]. Recently, Mo and his colleagues demonstrated that microglia promote aberrant neurogenesis and increased immature neuronal projections following seizures through its P2Y12 receptor by using current selective microglial targeting approaches [45]. P2RY13 is one of the most recently identified P2Y receptors. Studies on P2Y13 receptor knockout mice have highlighted specific functions in nervous system, including pain transmission, neuroprotection, regulation of microglial morphology, and intercellular calcium wave communication [46–49]. Lysophosphatidic acid receptor 5 (LPAR5) signaling has been mainly studied in the context of neuropathic pain [50,51], but Plastira and colleagues have recently shown that LPAR5 affects microglia biology and induces a distinct pro-inflammatory and migratory signature [52,53], and identified small-molecule antagonists that can interfere with LPA-induced pro-inflammatory signaling cascades in microglia [54].

Although the remaining three genes lack evidence to support their connection with microglia, they remain closely related to many neurological symptoms. The adrenergic α2C receptor (ADRA2C) is an inhibitory modulator of the sympathetic nervous system, which regulates neurotransmission at lower levels of neural activity in the cerebral cortex [55]. Meanwhile, the protein encoded by PDYN is a preproprotein that is proteolytically processed to form the endogenous opioid peptide dynorphin, which is a ligand for the kappa-type of opioid receptor (KOR) [56]. Extensive evidence suggests that dynorphin plays an important role as an endogenous anticonvulsant in epileptogenesis and epilepsy via KORs [57–61]. Preclinical data suggests a high potential of KOR agonists to control seizures [62]. SSTR3 encodes a member of the somatostatin receptor protein family with a wide range of physiological effects on neurotransmission, antiproliferative signaling and apoptosis [63,64]. Therefore, its related diseases mainly include pituitary adenoma and glioma [65–68].

Still, there were limitations to our research. First, we used a small number of samples to reflect the differential genes between epileptic brain and normal tissue, which may not be comprehensive. The uneven distribution of participants in terms of gender and age may also affect the accuracy of our results. Additionally, we selected KA-induced chronic epilepsy as an animal model to validate the reliability of sequencing data, which might not be entirely consistent with the process of epileptogenesis in humans. Thus, continuing to collect surgical specimens of epileptic brain and appropriate control tissues to achieve sufficient biological replication will be the focus of our next studies. We aim to obtain reliable biomarkers for in-depth mechanism exploration by combining the results of present study.

5.   Conclusions

In conclusion, by sequencing the transcriptome of normal and epileptic brains, screening for differentially expressed genes, and further analyzing and mining the corresponding PPI network, we obtained seven pivotal genes closely related to microglia activation and migration. These genes reveal the important role of microglia in the immune regulation and inflammatory response mediated by epilepsy, which may provide new candidates for the study of the molecular mechanisms of seizures and lay the theoretical foundation for further research.

Acknowledgments

This study was supported by the National Natural Science Foundation of China [No.81971217].

Conflict of interest

The authors declare there is no conflict of interest.