Systems-level analysis identifies key regulators driving epileptogenesis in temporal lobe epilepsy

Genomics xxx (xxxx) xxx–xxx

Contents lists available at ScienceDirect

Genomics journal homepage: www.elsevier.com/locate/ygeno

Original Article

Systems-level analysis identifies key regulators driving epileptogenesis in temporal lobe epilepsy Yingxue Fua,b, Ziyin Wua,b, Zihu Guoa,b, Liyang Chena,b, Yaohua Mab,c, Zhenzhong Wangb, ⁎ ⁎ Wei Xiaob, , Yonghua Wanga, a

Center for Bioinformatics, College of Life Sciences, Northwest A&F University, Yangling, Shaanxi 712100, China State Key Laboratory of New-tech for Chinese Medicine Pharmaceutical Process, Lianyungang, Jiangsu 222001, China c Lab of systems pharmacology, College of Life Sciences, Northwest University, Xi'an, Shaanxi 710069, China b

A R T I C LE I N FO

A B S T R A C T

Keywords: Epilepsy progression Epileptogenesis Network analysis Key regulators Antiepileptogenic targets

Temporal lobe epilepsy (TLE) is the most prevalent and often devastating form of epilepsy. The molecular mechanism underlying the development of TLE remains largely unclear, which hinders the discovery of effective antiepileptogenic drugs. Here we adopted a systems-level approach integrating transcriptomic profiles of three epileptogenesis stages to identify key regulators underlying epilepsy progression. Associating stage-specific gene meta-signatures with brain cell-specialized modules revealed positive regulation of glial migration and adhesion, cytokine production, and neuron death, and downregulation of synaptic transmission and ion transport during epileptogenesis. We identified 265 key regulators driving these processes and 72 of them were demonstrated associating with seizure frequency and/or hippocampal sclerosis in human TLE. Importantly, the upregulation of FAM107A, LAMB2, LTBP1 and TGIF1, which are mainly involved in nervous system development, were found contributing to both conditions. Our findings present the evolution landscape of epileptogenesis and provide candidate regulators that may serve as potential antiepileptogenic targets.

1. Introduction Epilepsy is a complex neurological disorder characterized by recurrent unprovoked seizures, of which temporal lobe epilepsy (TLE) is the most prevalent form [1]. The term epileptogenesis refers to the gradual process through which normal neuronal networks are altered resulting in the generation of chronic spontaneous seizures [2,3]. The process can be triggered by diverse brain insults, including traumatic brain injury, stroke, infections and prolonged seizures such as status epilepticus (SE), and is typically thought to involve three stages [4,5]. The first is the acute phase after the brain insult, in which a cascade of morphologic and biologic changes occurs in the injured area. This is followed by the latent period with no observed behavioral seizures, and subsequent chronic established epilepsy. Several molecular signaling cascades, such as mTOR, BDNF-TrkB and REST/NRSF pathways, have been demonstrated playing a role in epileptogenesis [6–8]. However, the key regulators driving the entire process of epilepsy progression remain largely unexplored. Applying high-throughput omics technologies to human epileptic tissues offers a great opportunity to unveil the crucial pathways or molecules that can contribute to the formation of seizure-generating ⁎

neuronal circuits. For example, a large-scale transcriptomic profiling of surgically resected hippocampi of TLE patients has recently been generated and used to identify gene-regulatory networks and regulators genetically associated with epilepsy [9,10]. However, there are obvious limitations and challenges in exploring epileptogenesis in human. One issue is that omics studies of human TLE generally lack appropriate control samples of healthy brain. Furthermore, the specimens collected from epilepsy patients are usually at an advanced stage and have been subjected to the treatment of antiepileptic drugs [11]. Alternatively, well-characterized animal epilepsy models which mimic prominent histopathological and electroencephalographic features of human TLE can be employed to examine the key molecular alterations during epilepsy progression [12]. A few reports have studied the genome-wide molecular changes throughout epileptogenesis using animal TLE models [13,14]. While other studies covered time points close to either acute responses to SE or cumulative effects of chronic spontaneous seizures [15,16]. As the modeling approaches and tissue dissection time varies across these studies, a systematic integration analysis of the existing datasets will likely provide a more comprehensive and robust molecular profiling for the epileptogenic process from the early hippocampal injury to the onset of chronic epilepsy.

Corresponding authors at: Center for Bioinformatics, College of Life Sciences, Northwest A&F University, Yangling, Shaanxi 712100, China E-mail addresses: [email protected] (W. Xiao), [email protected] (Y. Wang).

https://doi.org/10.1016/j.ygeno.2019.09.020 Received 10 July 2019; Received in revised form 31 August 2019; Accepted 25 September 2019 0888-7543/ © 2019 Elsevier Inc. All rights reserved.

Please cite this article as: Yingxue Fu, et al., Genomics, https://doi.org/10.1016/j.ygeno.2019.09.020

Genomics xxx (xxxx) xxx–xxx

Y. Fu, et al.

epileptogenesis. The VIPER algorithm was then used for inferring the differential activity of candidate regulators for each stage compared to control group. Key regulators were defined as those that exhibit a significant activity change. Finally, by interrogating the RNA-seq datasets of human TLE patients with different seizure frequencies and hippocampal sclerosis, the validity and clinical relevancy of these key regulators were investigated and proved.

Systems biology-based approaches that utilize network theory to organize transcriptome datasets have been used to prioritize candidate disease genes or to discern transcriptional regulatory programs [17–19]. One method to infer critical genes from gene expression data is the coexpression network analysis, which builds scale-free gene networks based on the pairwise gene expression correlations [20]. Genes with higher similarity scores tend to co-activate in a specific biological condition. Although coexpression analysis can identify genes or gene modules that associate with disease or biological phenotypes, it normally does not infer causality or distinguish between regulatory and regulated genes [21]. The algorithm for the reconstruction of accurate cellular networks (ARACNe) uses an information theoretic approach to eliminate most indirect interactions inferred by co-expression methods, leaving those expected to be regulatory [22]. Although originally being applied to infer the relationship between transcription factors (TFs) and their target genes, the method can also be adapted to infer the indirect transcriptional targets for other kind of regulators, such as signaling proteins [23]. Using these methods, the observed gene expression changes can be placed into a systems context that was related to the underlying disease biology. In this study, we proposed a systems-level framework which integrates gene meta-signatures, gene coexpression network and cellular regulatory network to reveal the evolution landscape of epilepsy progression and distinguish key regulators that govern the transition between different epileptogenesis stages (Fig. 1). The time-specific hippocampal transcriptomic profiles of rodent TLE models were collected and classified into acute, latent and chronic stages of epileptogenesis. These profiles were then subjected to the RankProd meta-analysis and utilized for generating stage-specific gene meta-signatures. Based on these gene signatures, a module association score (MAS) was proposed to depict the expression dynamics of functional modules during

2. Materials and methods 2.1. Data collection and preprocessing We searched the Gene Expression Omnibus (GEO) database using key words “temporal lobe epilepsy”, “TLE” or “MTLE”, and restricted the study type as “Expression profiling by array, or by high throughput sequencing”. The organisms of samples were limited to Homo sapiens, Rattus norvegicus and Mus musculus. After manually checking all resulting datasets, we obtained five microarray datasets of rodent TLE models that covered different time points following SE and two human TLE patient RNA-seq datasets with epilepsy symptom information (seizure frequency or hippocampal sclerosis). Details about these datasets, including accession numbers, platforms and references, were listed in Table 1. For microarray datasets, the series matrix files were downloaded and then subjected to quality assessment using the arrayQualityMetrics package from Bioconductor [24]. Outliers were identified using heatmaps and dendrograms based on inter-array expression distances, and also boxplots and density estimate plots. For samples from GSE27166 and GSE73878, which contain both sides of hippocampus, only the expression profiles of the ipsilateral hippocampi were included for analysis. For RNA-seq datasets, see the section “2.12 Human TLE patient RNA-seq data analysis”.

Fig. 1. Schematic overview of the study design. To identify key molecules driving epilepsy progression, we proposed an analytic framework which comprises four steps. Firstly, the hippocampal transcriptome datasets of rodent TLE models were collected and divided into the acute, latent and chronic stages based on the described tissue extraction time and phenotypes. Meta-analysis was then performed using the RankProd method to evaluate differential gene expression and to generate the gene meta-signature for each stage. Secondly, to elucidate the functional organization of genes under the epilepsy context, gene coexpression network was constructed and used for identifying functional modules. A module association score (MAS) was then defined to quantify a module's association degree with an epileptogenesis stage. Thirdly, to identify key regulators controlling the transition between stages, a gene regulatory network was constructed and the VIPER algorithm was employed to infer regulator activity changes in epileptic samples. Finally, using the RNA-seq profiles of human TLE patients with recurrent seizures and hippocampal sclerosis, key regulators associated with both conditions were screened out. 2

Genomics xxx (xxxx) xxx–xxx

Y. Fu, et al.

[13] [16] [14] NA [15] NA [62]

2.3. Differential gene expression analysis for individual datasets For individual microarray datasets, we used the limma [26] and RankProd (RP) [27] packages from Bioconductor for differential expression analysis between control samples and epileptic samples of different epileptogenesis stages. The limma approach compare groups of samples by fitting gene-wise linear models and applying empirical Bayes methods to identify differentially expressed genes (DEGs). The genes with absolute log2FC (fold change) > 0.5 and adjusted P-value for multiple comparison (false discovery rate, FDR) < 0.05 were considered significantly differentially expressed. RP is a non-parametric statistical method used to detect variables consistently upregulated or downregulated in replicate samples. It provides several advantages over linear modeling, including the biologically intuitive criterion, fewer model assumptions, and increased performance with noisy data. The DEGs were identified based on the same thresholds as limma (absolute log2FC > 0.5 and FDR < 0.05). The list of DEGs identified by the two methods for each dataset was compared using Venn diagrams created by jvenn [28].

+ +

+

GE Bioarray GE Bioarray Affymetrix Array Illumina Beadchip Affymetrix Array Illumina HiSeq Illumina HiSeq Pilocarpine Intrahippocampal KA Electrical stimulation of amygdala Intrahippocampal KA Ditto TLE with low or high seizure frequency TLE with and without hippocampal sclerosis +

+ +

For each epileptogenesis stage, the gene expression matrices of samples belonging to the corresponding stage were integrated from multiple datasets and meta-analysis was performed to evaluate the differential gene expression using the RankProd package [27]. Though the RP method was initially developed to detect DEGs in a single experiment, it is able to integrate datasets from multiple origins and overcome the heterogeneity among them because of the use of ranks instead of actual expression values. The four microarray platforms GPL1261, GPL2896, GPL6247 and GPL6885 for the rodent TLE model datasets contain 21,720, 12,733, 15,124 and 17,125 unique genes, respectively, of which 9139 genes were common across all platforms. The expression values of these common genes in each dataset were then extracted. For multiple probes that correspond to the same gene in a dataset, the probe with maximum mean expression values was retained to represent that gene. The RP method was then applied to the combined datasets of each epileptogenesis stage to assess the differential expression of genes. Meta-DEGs of each stage were determined under the thresholds of absolute log2FC > 0.5 and FDR < 0.05. Besides, we also extracted the whole ranked gene lists of the three stages for further analysis. As RP employs separate ranks for up- and down-regulated genes, we integrated the two rank lists using the following equation (Eq. (1)).

Human

Mouse

KA, kainic acid; NA, not available.

+ + +

+ +

0 14 0 45 0 0 0 Rat

GSE14763 GSE27166 GSE49849 GSE73878 GSE88992 GSE127871 GSE71058

18 24 20 80 17 12 22

Acute phase

Latent period

Chronic epilepsy

Platforms Model type/Disease Number of excluded samples

Epileptogenesis stages

Unsupervised principal component analysis (PCA) was performed to further visualize the correlations among samples belonging to different epileptogenesis stages. All datasets were normalized and log2 transformed, and then analyzed using the prcomp function from the “stats” module in R. PCA methodology captures the inherent gene expression patterns in the data by projecting multivariate data objects onto a lower dimensional space while retaining much of the original variance [25].

2.4. Gene meta-signatures for specific epileptogenesis stages

Number of all samples Dataset accession ID Species

Table 1 List of rodent TLE model microarray datasets used for epileptogenesis analysis and human TLE patient RNA-seq datasets used for validation.

References

2.2. Principal component analysis

rnorm

⎧ 1 − rup − 1 , average log FC > 0 2 ⎪ max(rup ) = ⎨ rdown − 1 , average log 2 FC < 0 ⎪− 1 + max(r down ) ⎩

(1)

The ranked gene lists with the normalized rank (rnorm) and average log2FC were then utilized as gene meta-signatures for the three epileptogenesis stages. 3

Genomics xxx (xxxx) xxx–xxx

Y. Fu, et al.

for multiple gene sets (modules) was performed via Metascape [31]. Pathway and process enrichment analysis was carried out with the following ontology sources: GO Biological Processes, KEGG Pathway, Reactome Gene Sets and CORUM. All genes in the genome have been used as the enrichment background. Redundant terms were clustered into groups based on their similarities and the top 20 scored clusters were used as the final functional annotation for modules.

2.5. Text mining to determine gene-epilepsy and -epileptogenesis associations To determine whether a meta-DEG has potential association with epilepsy or epileptogenesis, we applied text mining technique by searching certain terms in MEDLINE documents (3,811,179 abstracts) using SCAIView-neuro webserver (http://neuro.scaiview.com/neuro//, Version 1.7.3). Specifically, we searched “Epilepsy” or “Epileptogenesis” with the restriction of [Epilepsy], resulting in 69,521 and 8744 abstracts, respectively. We then searched each gene with the [Human Genes/Proteins] restriction, and also gene's co-citation with “Epilepsy” or “Epileptogenesis”, and recorded the number of documents of each search. The association significance of a gene-epilepsy or -epileptogenesis pair was then quantified using hypergeometric test to determine if that gene’s co-citation with epilepsy or epileptogenesis was more frequent than expected by chance. Gene-epilepsy or -epileptogenesis pair with a P-value < 0.05 was considered as an established association.

2.8. Module association score (MAS) with different epileptogenesis stages To evaluate the association degree of modules with a specific epileptogenesis stage, a module association score (MAS) was defined to reflect the overrepresentation degree of a module at the extremes (top or bottom) of the ranked gene meta-signatures. For each module, only genes with the scaled intramodular connectivity > 0.2 were kept to represent the module. This reduced the noise for calculating MAS as genes with lower connectivity within a module typically contribute less to its functions. The MAS and the corresponding significance level were then calculated using the fgsea R package [32] against the gene metasignature of each epileptogenesis stage. The fgsea method implements a special algorithm for fast gene set enrichment analysis. The significance of gene set enrichment was determined using the empirical enrichment score null distributions simultaneously calculated for all the gene set sizes. Module MAS with the adjusted P-value < .05 was considered to be significantly associated with a specific epileptogenesis stage.

2.6. Gene coexpression network construction and module detection To construct a epileptogenesis-context gene coexpression network, we subjected the dataset containing the entire process of epileptogenesis to weighted gene coexpression network analysis (WGCNA) [18,20]. To overcome outlier bias, a robust correlation measure, biweight midcorrelation, was used to quantify the co-expression similarity sij between each pair of genes [29]. Then, a weighted network adjacency matrix A = [aij] was computed by applying a power function on all positive gene correlations, and was set to be zero when two genes have negative (or zero) correlations (Eq. (2)).

sij β sij > 0 aij = ⎧ ⎨ 0 sij ≤ 0 ⎩

2.9. Construction of gene regulatory network Candidate gene regulators were collected from three aspects. 1) Transcription factors or regulators (TFs). Transcription factors were obtained by extracting genes annotated in GO molecular function as GO:0003700, “DNA binding transcription factor activity”. Transcription regulators were the intersection of genes annotated as GO:0003677 “DNA binding” and genes annotated with GO:0140110, ‘transcription regulator activity’ or GO:0003612 “transcription coregulator activity” or GO:0006355, “regulation of transcription, DNA-templated”. 2) Synaptic proteins (SPs). SPs are regarded as genes annotated in GO cellular component as GO:0045202, ‘synapse’ or GO:0030424, “axon” or GO:0030425, “dendrite”. 3) Signaling proteins (Signal), which were built upon genes annotated in GO Biological Process GO:0007165 “signal transduction” and not overlapping with above two gene list. The genes corresponded to these GO terms were extracted using the biomaRt package [33]. To construct the gene regulatory network under the epileptogenesis context, we used the transcriptomic dataset that contains samples of all three epileptogenesis stages, i.e. GSE73878. Those candidate regulators along with the gene expression matrix were then subjected to the ARACNe-AP software [34] for reverse engineering a gene regulatory network. ARACNe was run with default parameters: 0 DPI (data processing inequality) tolerance and MI (mutual information) P-value threshold of 10−8, with 100 bootstrap iterations. Here, the MI threshold corresponded to the P-value threshold of 10−8 was 0.573. Only interactions with MI above this threshold were included in the final network.

(2)

This ensures the connections of all gene pairs have the same direction and reduces the strength of weak correlations while preserving connection strength of highly correlated genes. The connectivity k for n each gene was then defined as ki = ∑ j ≠ i aij and was used for the network analysis. To balance the scale-free topology (i.e. p(k) ~ k-γ) and the sparsity of connections between genes in the network, a set of β values was evaluated to obtain the optimal specificity and sensitivity. To detect modules from the gene coexpression network, the topological overlap matrix (TOM), which reflects the relative interconnectedness between each pair of genes, was calculated. Based on the topological overlap dissimilarity (1 − TOM) between genes, a gene dendrogram was generated using average hierarchical clustering. The dynamic branch cutting method was then used for detecting gene coexpression clusters (modules) in the dendrogram depending on its shape. Module eigengene, which is the first principal component of gene expression, was calculated to summarize the gene expression within a module and to merge modules with high similarities. 2.7. Cell-type and functional enrichment analysis For the cell-type enrichment, marker genes of nine brain cell types were obtained from the single-cell RNA-seq profiles of the mouse cortex and hippocampus [30]. These include three types of neurons (cortical pyramidal neuron, CA1 pyramidal neuron and interneuron), four types of glia cells (astrocyte, oligodendrocyte, microglia and ependymal cell), and the vascular endothelial and mural cells. The number of marker genes of these cells ranges from 155 to 484 genes, leading to a total of 3189 unique marker genes (Supplementary Table S1). Enrichment between modules and the cell-type marker genes was measured using the hypergeometric test with subsequent Benjamini-Hochberg (BH) correction for multiple comparison as implemented in the userListEnrichment function in the WGCNA package [20]. Functional meta-analysis

2.10. Differential activity inference for gene regulators To infer the relative activity of the gene regulators at different epileptogenesis stages, we applied the VIPER algorithm [23] to test for regulon (a group of genes that are regulated by the same regulator) enrichment on stage-specific gene meta-signatures. VIPER uses a probabilistic framework that integrates target mode of regulation (i.e., activated, repressed or undetermined represented by an index ranging from −1 to 1), statistical likelihoods of regulator-target interactions and target overlap between different regulators (pleiotropy). To compute the enrichment of a protein's regulon in ranked gene lists, an analytic rank-based enrichment analysis (aREA) method, which 4

Genomics xxx (xxxx) xxx–xxx

Y. Fu, et al.

Fig. 2. Principal component analysis (PCA) for microarrays of rodent TLE models. The variances captured by the first two PCs are shown along the respective axes. Samples belonging to the acute, latent and chronic stages are colored in red, blue and green, and control samples in gray. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

function implemented in the ‘edgeR’ [37] package. The original count data was then transformed to the log2 scale by using the rlog (regularized log transformation) function in DESeq2. To quantify the association between genes and seizure frequency, a Spearman's rank correlation coefficient was calculated between the log2 transformed counts of genes and the seizure frequency (SF) of patients. Genes with coefficient P-value < 0.05 were considered to be significantly associated with SF. For the hippocampal sclerosis (HS) dataset (GSE71058), the filtering threshold for genes was set to CPM > 0.1 in more than eight samples, which is the minimum number of samples in the two groups (with and without HS). Before performing differential expression analysis, we applied the RUVg function from RUVSeq package [38] to estimate the factors of unwanted variation using a set of negative control genes, i.e., genes that can be assumed not to be influenced by HS. The control genes here were defined as those least significantly DEGs based on edgeR analysis. We then performed differential expression analysis with DESeq2 [39] by considering a design matrix that includes both the HS status and the estimated factors of unwanted variation. Genes with the absolute log2FC > 0.5 and adjusted P-value < .05 were considered as HS-related genes. To make cross-species gene mappings, we standardized gene identifiers from microarray probe identifiers to NCBI Entrez ID identifiers and mapped mouse Entrez ID identifiers to their human ortholog using the biomaRt [33] package.

conduct a statistical analysis based on the mean of ranks, was used. The normalized enrichment score (NES) computed by aREA for each regulon were employed to quantitatively represent the regulators' relative protein activity in an epileptogenesis stage compared to the control group. For each stage, regulators with an absolute differential activity score > 2 and corresponding P value < .05 were defined as key regulators. 2.11. Integration analysis of synaptic signaling pathways Pathway enrichment analysis for key regulators was performed using DAVID v6.8 [35] against the KEGG pathway database [36]. To investigate the synaptic signaling variations at different epileptogenesis stages, an integrated synaptic signaling pathway was constructed. The integrated pathway is composed of key regulators enriched in the synapse-related pathways (i.e., the dopaminergic, cholinergic, glutamatergic, serotonergic and GABAergic synapses), and critical intracellular signaling pathways including MAPK signaling, calcium signaling, cGMP-PKG signaling and Ras signaling pathways. The interactions between these key regulators and their regulatory relationships with biological functions were extracted from the KEGG pathway database. 2.12. Human TLE patient RNA-seq data analysis

3. Results

The matrices of raw gene counts of the two RNA-seq datasets were downloaded from the GEO database. For the dataset of TLE patients with different seizure frequencies (GSE127871), genes with the countper-million (CPM) value > 0.5 in > 50% of the samples were retained for further analysis. The CPM of genes was calculated using the cpm

3.1. Gene meta-signatures associated with different epileptogenesis stages To investigate the molecular profiles underlying epileptogenesis, we 5

Genomics xxx (xxxx) xxx–xxx

Y. Fu, et al.

Fig. 3. Network functional dynamics during epileptogenesis. A. Dendrogram showing clustering of 8384 genes based on the topological overlap dissimilarity of genes in the gene coexpression network. Bottom colour bar indicates the 13 gene coexpression modules (M1 − M13) and their corresponding sizes (i.e. the number of genes in a module). B and C. Heatmaps showing the functional meta-analysis results for modules enriched with glia marker genes (M1, 2, 4 and 5) (B) and modules enriched with neuron marker genes (M3, 6, 7, 8, 9, 10 and 11) (C). The top 20 enriched functional terms are shown as rows and modules shown in columns. The heatmap is colored by the P-values. D. The MASs and corresponding significance levels (−log10(adjusted P-value)) of modules at the three epileptogenesis stages. The size of circle is proportional to the number of genes in each module used for calculating MAS.

epilepsy (Fig. 1 and Table 1). There were four, three and two datasets that contain samples of the acute, latent and chronic stages, respectively. After data normalization and preprocessing, 99 expression profiles were obtained and the individual datasets were then subjected to the PCA analysis. As shown in Fig. 2, epileptic samples belonging to acute and latent stages exhibited good separation from control samples along the first two PCs in different datasets, while samples of chronic epilepsy stage were close to control samples. To evaluate the differential gene expression between control samples and epileptic samples of each epileptogenesis stage, we first

used the time-specific hippocampal transcriptomic data of rodent TLE models from five independent studies (Table 1). These studies contain both rat and mouse TLE models and the modeling approaches were various, including systematic administration of pilocarpine, intrahippocampal KA injection and electrical stimulation of amygdala. The epileptic samples from these datasets covered a wide range of time points after SE, ranging from hours, to days and months. Based on the described tissue extraction time and phenotype, the samples in each dataset were divided into the control group and groups of three epileptogenesis stages, i.e., the acute phase, latent period and chronic

6

Genomics xxx (xxxx) xxx–xxx

Y. Fu, et al.

marker genes, modules M3 and M9 were mainly involved in functions of “regulation of ion transport”, “regulation of vesicle-mediated transport”, “modulation of chemical synaptic transmission”, “signal release” and “behavior” etc. (Fig. 3C). M8 was significantly enriched in “secondmessenger-mediated signaling” and “regulation of system process”, while modules M6, 7 and 10 participated in processes of “regulation of intracellular transport”, “signal release” and “tissue morphogenesis”. The consistency between the cell-type specificity and functional annotation of modules demonstrates the biological significance of the identified modules in the context of epileptogenesis. To further investigate how the expression of these modules was regulated at different epileptogenesis stages, we defined a module association score (MAS) to reflect the degree of which a module was enriched at the top or bottom of the stage-specific gene meta-signatures. All modules were significantly associated with at least one epileptogenesis stage, among which five modules (M1, 8 and M5, 6, 7) were consistently up- or down- regulated in all three stages (Fig. 3D). The two positively associated modules M1 and 8 were mainly related to the inflammatory responses and increased intracellular signaling activity. While the negatively associated modules M5, 6 and 7 may imply an impairment of the synapse function after SE. M4, 3 and 13 exhibited a specific association with the acute, latent and chronic stage, respectively. Module M9 was downregulated in the acute and latent phases but not the chronic epilepsy stage. Overall, the dynamic changes of modules' association with different stages provide us a global landscape of the evolution process of epileptogenesis.

applied the limma method [26] to individual datasets. The DEGs were defined as those that achieved an absolute log2FC > 0.5 and a FDR < 0.05 between control and epilepsy. For the acute and latent stages, the limma method yielded gene lists with very small overlap across the datasets (43 out of 5281 genes for acute stage, and 25 out of 1373 genes for latent stage) (Supplementary Fig. S1A). While for the chronic phase, no DEGs were detected in one of the two datasets. We then applied another differential expression analysis method RankProd (RP) [27] that differs from the linear modeling-based approaches. RP is a rank-based technique that detects genes that consistently appear among the most highly ranked genes (either strongly upregulated or downregulated) in a number of replicate samples. Under the same threshold (absolute log2FC > 0.5 and FDR < 0.05), the RP method depicted a slightly better but still small overlap of DEGs in latent and chronic stages (Supplementary Fig. S1B). These results suggest that direct comparison across individual datasets was not feasible due to the heterogeneity of experimental approaches and profiling platforms. Since the RP technique transforms the actual expression values into ranks, it has the ability to handle variability among datasets from multiple origins [40]. We thus adopted the RP for meta-analysis for each of the three epileptogenesis stages. We obtained 462, 201 and 82 meta-DEGs (absolute log2FC > 0.5 and FDR < 0.05) for acute, latent and chronic stages, respectively (Supplementary Table S2 and Fig. S2). 40 meta-DEGs were dysregulated across all three epileptogenesis stages, among which the upregulation of Gfap (Glial fibrillary acidic protein) and Bdnf (Brain-derived neurotrophic factor) has been widely described as biomarkers of epileptogenesis [41,42]. To investigate whether other meta-DEGs have been related to epileptogenesis or epilepsy processes, we performed text-mining based on the published biomedical literature using the SCAIView-neuro webserver. Among the 40 meta-DEGs, five of them showed established gene-epilepsy associations, and eight genes have been related with epileptogenesis (Supplementary Table S3), suggesting that the meta-DEGs captured by our analysis may be more specifically associated with epileptogenesis. Finally, to provide a full molecular depiction for each epileptogenesis stage, we extracted the whole ranked gene lists ordered by the average log2FC of genes to serve as the stage-specific gene meta-signatures (Supplementary Table S4).

3.3. Identification of gene regulators driving epileptogenesis To discover gene regulators controlling the transition from acute to latent and chronic stages of epileptogenesis, we interrogated the timespecific hippocampal transcriptome profiles using the VIPER algorithm [23] to infer the activity change of regulators in a specific stage. VIPER infers protein activity by systematically analyzing the expression of a protein's regulon, which refers to the transcriptional targets of that protein. The ARACNe technique [34], which detect maximum information path targets, was used to systematically infer regulons, resulting in an gene regulatory interactome of 41,364 interactions between 1493 regulators and 5695 target genes. VIPER then compute the enrichment of a protein's regulon in gene expression signatures based on a probabilistic framework that directly integrates target mode of regulation, regulator-target interaction confidence and target overlap between different regulators. Differential activities of 521 regulators were obtained for each of the three epileptogenesis stages. Though we observed a strong association between regulator's differential activity and the average log2FC of genes derived from meta-analysis, VIPER-based activity inference may point out critical regulators that are not differentially expressed at the mRNA level (Supplementary Fig. S4). For each stage, key regulators were defined as those with absolute differential activity score > 2 and corresponding P value < .05, which represents a significant activity alteration compared to the control group (Fig. 4A and Supplementary Fig. S5). There were 214, 198 and 156 key regulators respectively associated with the acute, latent and chronic stage (Fig. 4B). 43% of the key regulators were dysregulated in all three stages, indicating that these regulators were immediately involved in the epileptogenesis following SE, and exhibited continuing changes extending into the chronic epilepsy period. Modules M1, 3, 6, 7 and 8 have the highest numbers of key regulators (Supplementary Fig. S6). Regulator activities in M1 and M8 were both upregulated during the epileptogenesis, yet their activities had opposite changing trends (Fig. 4C). Whereas M1 regulators were mainly associated with acute response after SE, M8 regulators may play major roles in the latent and chronic epilepsy stages. Regulators in M3, 6 and 7 show constant downregulated activity at all three epileptogenesis stages (Fig. 4C). The module-based changing pattern of regulator's differential activity offers a straightforward illustration of

3.2. Network functional dynamics during epileptogenesis The gene meta-signatures associated with different epileptogenesis stages provide a significant starting point for dissecting the development of epilepsy. However, it is difficult to pinpoint epileptogenic mechanisms without considering the functional organization of these genes. To this aim, an epilepsy-context gene coexpression network was built using WGCNA [20] based on the dataset that contains samples of all three epileptogenesis stages. 13 gene modules (M1 − M13), with their size varied from 44 to 1896 genes, were identified (Fig. 3A). To investigate the cell-based context of these modules, we performed enrichment analysis against marker genes of nine brain cell types including neurons, glial cells and endothelial cells which were derived from single-cell RNA-seq analysis of the mouse cortex and hippocampus [30]. Module M1, the largest module with 1896 genes, was significantly enriched with three types of glial cells, i.e., microglia, astrocyte and oligodendrocyte (BH-corrected P = 4.0E−55, P = 7.9E−7, and P = 2.7E−5, respectively) (Supplementary Fig. S3). Modules M3, 6, 7, 8 and 9 were specifically enriched for pyramidal neurons and/or interneurons (BH-corrected P value range: 3.7E-2 to 2.7E-25). And modules M12 and 13 were respectively enriched with the ependymal cell and mural cell (BH-corrected P = 4.4E-11 and 2.2E-4). Functional meta-analysis of these modules showed that module M1 mainly participated in GO biological processes of “positive regulation of cell migration”, “cytokine production”, “regulation of cell adhesion” and “apoptotic signaling pathway” (Fig. 3B). Besides, M1 and M2, 4, 5 were all enriched with items of “neuron death” and “regulation of cellular response to stress”. Consistent with the enrichment for neuronal 7

Genomics xxx (xxxx) xxx–xxx

Y. Fu, et al.

Fig. 4. Identification of gene regulators driving epileptogenesis. A. Heatmap of the differential activity of 265 key regulators (KRs) at the three epileptogenesis stages compared to the sham control group, with red denoting upregulation and blue denoting downregulation. The modules and types of these regulators are displayed at the left. SP, synaptic protein; Signal, signaling protein; TF, transcription factor. B. Venn plot showing the overlap of KRs in the three epileptogenesis stages. C. Module-based regulator activity dynamics during epileptogenesis indicated by the mean and standard deviation of relative activities of KRs in a module. D. KEGG pathway enrichment analysis for KRs. The y-axis represents the top 15 significantly affected canonical pathways and x-axis the −log10 transformed BH adjusted Pvalues. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

(adjusted P-value = 1.7E−03) which can be activated by nerve growth factor (NGF) and BDNF is an important pathway involved in the survival, development, and function of neurons. Other enriched pathways include the MAPK signaling, Calcium signaling, cGMP-PKG signaling and Ras signaling pathways (adjusted P-value range = 3.4E−04 to 2.3E−03), which are critical intracellular signaling pathways related to multiple cellular functions.

how the regulators were modulated along with the development of epilepsy. To understand the molecular processes affected by these key regulators, we performed enrichment analysis against the KEGG pathways database [36]. The analysis showed that the key regulators were involved in signaling pathways related to chemical synaptic transmission, immune response, growth factor signaling, and pathways related to cell proliferation and death (Fig. 4D). Among the synaptic transmission -related pathways, Dopaminergic synapse was the top enriched pathway (adjusted P-value = 1.4E−05), followed by Cholinergic, Glutamatergic, Serotonergic and GABAergic synapses. Besides, the Retrograde endocannabinoid signaling, which can suppress both excitatory and some inhibitory synapses, was also enriched with the key regulators (adjusted P-value = 1.7E−04). Endocannabinoids and their receptors are altered by epileptic seizures and can in turn control key epileptogenic circuits by inhibiting synaptic transmission in the hippocampus [43]. Multiple immune response-related pathways were also highly enriched, including the Toll-like receptor signaling, Chemokine signaling and TNF signaling pathways (adjusted P-value range = 8.7E−05 to 8.3E−04). The Neurotrophin signaling pathway

3.4. Variations of synaptic signaling at different epileptogenesis stages A thorough knowledge of signaling pathways involved in both acute- and long-term responses to SE is crucial to unravel the origins of epilepsy [42]. To better understand the regulatory mechanism of synaptic transmission underlying epileptogenesis, we constructed an integrated synaptic signaling pathway that composed of the key regulators involved in synapse-related functions (Fig. 5A). A heatmap of the inferred differential activities of these key regulators in the three epileptogenesis stages was shown in Fig. 5B, in which the regulator types and modules were marked on the top. We first examined the relative activity changes of key regulators 8

Genomics xxx (xxxx) xxx–xxx

Y. Fu, et al.

Fig. 5. Variations of synaptic signaling at different epileptogenesis stages. A. The integrated synaptic signaling pathway showing the relationship between KRs and their activity changes during epileptogenesis (red, upregulated KR; blue, downregulated KR; white, other proteins). B. Heatmap showing differential activity of KRs of the synaptic signaling pathway in three epileptogenesis stages. Dysregulated KRs in the pathway are shown as rows and columns show differential activity in the acute, latent and chronic stage, respectively. The types and modules of these KRs are shown on the top. SP, synaptic protein; Signal, signaling protein; TF, transcription factor. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

SHC1 and GAB1, and PLCγ (PLCG2). This result is consistent with previous studies reporting that excessive activation of TrkB caused by SE promotes the development of TLE [45,46]. Activities of phospholipase C (PLC) isotypes (PLCB2, PLCB3, and PLCD4) were also upregulated, which further activates the calcium signaling and protein kinase C (PKC), and downstream transcription factors NFKB1 and RELA. For MAPK signaling, though there were increased activity of upstream Ras (RRAS), Raf (ARAF) and MAP3K8, the JNK1 (MAPK8) activities were downregulated, along with the downregulation of activating transcription factor 2 (ATF2). Finally, the PI3K-AKT-CREB pathway, CaM kinase (CAMK2B) and calcium binding protein SCGN all exhibited decreased activity. In sum, these results demonstrate that the synapse-tonucleus signaling underlying epileptogenesis displayed a complex restructured system related to neuronal hyperexcitation and impaired synaptic plasticity.

located on synaptic membrane. The kainate receptor (KAR) subunit 1 (GRIK1), which belongs to the family of ionotropic glutamate receptors, was strikingly downregulated in the acute and latent stages. The β3 and γ2 subunits of GABAA receptor (GABRB3 and GABRG2) showed opposite changing trends during the epileptogenesis. GABRB3 was significantly downregulated in the acute phase, while GABRG2 was upregulated in the latent and chronic periods. Besides, decreased activity of acetylcholine receptors (CHRNB2 and CHRM3) and serotonin receptor (HTR1A) was also observed. Among these regulators, mutations of GABRB3, GABRG2 and CHRNB2 genes have been reported associating with some familial epilepsy syndromes [44]. Notably, increased activity of glutamate transporters (SLC1A2 and SLC1A3) and decreased activity of GABA transporters (SLC6A1 and SLC32A1) further support the idea that imbalance between excitation and inhibition and altered threshold for neural excitation are underlying epileptic behaviors. The voltage-gated potassium channel Kv3.1 (KCNC1, alpha subunit) was markedly downregulated in the acute stage, implying the inability of the neuron to normally depolarize following SE. The accessory subunits of Kv1 (KCNAB1 and KCNAB2) and Kv4 (KCNIP3) also display stagespecific activity changes. Constant upregulation of BDNF and TrkB (NTRK2) occurred during the entire epileptogenesis, which further activated the adaptor proteins

3.5. Key regulators associated with seizure frequency and hippocampal sclerosis in human TLE One of the direct outcomes of the epileptogenesis is the presence of spontaneous recurrent seizures. To investigate whether the identified key regulators were involved in controlling seizure frequency in 9

Genomics xxx (xxxx) xxx–xxx

Y. Fu, et al.

understand the molecular mechanisms underlying epileptogenesis, it is necessary to take advantage of multiple animal epilepsy models with different origins and phenotypes. In this study, we performed an integration analysis on time-specific transcriptome profiles of various rodent TLE models, which include both chemical and electrical kindling models of epilepsy. Direct comparison of DEGs identified from individual datasets led to very limited information, indicating there need more systematic analysis to detect genes that are consistently up- or down-regulated across datasets of different origins. By applying a rankbased meta-analysis method RP to the gene expression matrix of each epileptogenesis stage, we obtained a set of meta-DEGs related to acute, latent and chronic stages, respectively. Besides, we also extracted stagespecific gene meta-signatures to depict the molecular features of the three epileptogenesis stages at a genome-wide level. Compared to using only DEGs as molecular signatures for a specific phenotype, gene lists with the rank information and relative fold change of all genes provide a more comprehensive molecular representation. Utilizing these gene list-based meta-signatures, a MAS value was proposed to assess the molecular functional dynamics between stages. In comparison with previous transcriptomic studies of epilepsy or epileptogenesis, our analysis recaptured the key genes/regulators associated with specific stages of epileptogenesis. For example, Kalozoumi et al. [15] conducted a genome expression analysis to determine the early mechanisms of epileptogenesis (6 h, 12 h and 24 h post-SE), and identified five overexpressed master regulator genes implicated in glial responses, i.e., Cyba, Cybb, Lgals1, Nedd9 and Vim. Interestingly, four of these genes were also present in our upregulated meta-DEGs associated with the acute stage (Supplementary Table S2). In addition, both this study and another transcriptomic study of chronic epilepsy (10 months after SE) [16] point out that glial oxidative stress holds a key role during epileptogenesis. This is in consistent with our observation that the glia-enriched module M1, which was involved in multiple immune response processes, including “response to oxidative stress” (Fig. 3B), was constantly upregulated in all three epileptogenesis stages (Fig. 3D). Although most studies had focused on the hyper-expressed genes or upregulated pathways in epileptogenesis [13–16], our analysis also revealed several downregulated functional modules (M3, M6 and M7), which were enriched in pyramidal neurons and/or interneurons, suggesting that an impaired synaptic transmission and circuit-level dysfunction in the hippocampus may also underline epileptogenesis [3]. Complex disease phenotypes are typically drove by defects in multiple gene regulators which exhibit concurrent and aberrant activities [50]. For identifying key gene regulators involved in epileptogenesis, we inferred the differential activity of regulators in each the three stages. The regulator types include not only TFs, which are commonly regarded as the direct regulators controlling transition between different biological conditions [51], but also proteins on the synaptic membrane and intracellular signaling proteins. Given that various proteins located at the membrane of pre- or post-synapse are under the most directly impact when seizure activity occurs, the inclusion of synaptic and signaling proteins can help better depict the abnormalities of the synapse-to-nuclear signaling underlying epileptogenesis. Among the key regulators present in the integrated synaptic signaling pathway in Fig. 5A, several have already been validated in animal epilepsy models for their activity changes. For instance, increased immunolabeling of the γ2 subunit of the GABAA receptor (GABRG2) had been observed throughout the hippocampal formation from 7 to 60 d after SE in a mouse epilepsy model, and the density of labeling increased most significantly in 60 d [52]. Another study of mouse epilepsy model has revealed that downstream PLCγ1 activation is the dominant signaling effector by which excessive activation of TrkB promotes epilepsy [45]. For MAPK signaling, the JNK2/3 isoform phosphorylation was reported to be decreased in the chronic period (50 days after SE) in the hippocampus [53]. Some other key regulators' activity changes can be backed up by their genetic associations with epilepsy. The loss-of-function

patients with TLE, we utilized an RNA-Seq dataset of the hippocampal tissue resected from 12 medically intractable TLE patients with seizure frequencies ranging from 0.33 to 120 seizures per month. Genes with very low counts across samples were filtered out based on their CPM values, and a final set of 18,057 genes were obtained for further analysis. To quantitively assess the association between genes and seizure frequency, we calculated the Spearman's rank correlation coefficient between the log2-transformed counts of genes and the seizure frequency of patients. A corresponding P-value for the coefficient < 0.05 was considered as a significant correlation, resulting in 2125 genes whose expression was associated with SF (Supplementary Table S5). 50 of these genes were key regulators, including above-mentioned PLCB3 (r = 0.81, P = 1.4E-3), NTRK2 (r = 0.58, P = 0.047), SHC1 (r = 0.61, P = 0.034), HTR1A (r = −0.62, P = 0.03) and CHRM3 (r = −0.63, P = 0.026). It is noteworthy that all the five key regulators exhibited consistent changing directions in epileptogenesis and seizure frequency, that is, upregulated activity in epileptogenesis corresponds to a positive correlation with seizure frequency, and downregulated activity points to negative correlation. Indeed, we found that the differential activity of 94% of these key regulators (47 out of 50) in epileptogenesis were in accord with their sign of correlation coefficient (Supplementary Fig. S7A). These results demonstrate the validity of the identified key regulators, and also their conservation in human TLE. We further investigated whether the key regulators were also associated with another common epilepsy pathological condition, hippocampal sclerosis. HS is featured with severe neuronal cell loss and gliosis in the hippocampus, and can be both the cause and outcomes of the epileptogenesis [47]. Utilizing the RNA-seq profiles of dentate granule of MTLE patients with and without HS, we studied how the expression of epileptogenesis-associated key regulators were regulated by HS. After filtering low-expressed genes with CPM, 22,020 genes were retained. Before conducting differential expression analysis between the no-HS and HS groups, we applied a data normalization method, RUVSeq. [38], to remove the unwanted variation from the RNA-Seq data (Supplementary Fig. S8). Then, by taking into account both the HS status and the estimated factors of unwanted variation, differential expression analysis was performed with DESeq2 [39]. This resulted in 1470 HS-related DEGs (absolute log2FC > 0.5 and FDR < 0.05), which contain 32 key regulators (Supplementary Table S6). 84% (27 out of 32) of the key regulators showed consistent expression changes in HS compared with their epileptogenesis activities, including KCNIP3 (log2FC = 1.28), NTRK2 (log2FC = 0.90), RRAS (log2FC = 1.28) and MAPK8 (log2FC = −0.56) (Supplementary Fig. S7B). By checking those TLE phenotype-associated key regulators, we discovered 10 regulators that were associated with both seizure frequency and hippocampal sclerosis, namely ASIC2, C1QA, CNN3, FAM107A, HSPB1, LAMB2, LTBP1, NTRK2, SPARC and TGIF1 (Fig. 6). Neural related functions of these regulators were listed in Table 2. Six of them have been previously reported associating with epilepsy or epileptogenesis, while the others may need further investigation. Only ASIC2 (Acid-sensing ion channel 2) exhibited opposite changing trends in SF and HS conditions (Fig. 6). ASIC2 displayed decreased activity during epileptogenesis and had negative correlation with seizure frequency, while it was significantly upregulated in the HS group, consistent with a recent study reported elevated ASIC2a expression in both rat epilepsy models and human TLE patients [48]. Altogether, 72 key regulators of epileptogenesis were associated with seizure frequency and/or hippocampal sclerosis, indicating that our systems-level analysis has provided a valuable set of genes that may serve as potential therapeutic targets for epilepsy. 4. Discussion Epilepsy is a heterogeneous disorder with multiple origins and many different mechanisms of pathogenesis among patients [49]. To better 10

Genomics xxx (xxxx) xxx–xxx

Y. Fu, et al.

Fig. 6. Key regulators associated with seizure frequency and hippocampal sclerosis in human TLE. A. Scatter plots showing the correlations between KRs' expression and seizure frequency of patients. r = Spearman's rank correlation coefficient between the log2(counts) of genes and seizure frequencies. The blue lines are the regression lines fitted by linear modeling. B. Boxplots showing the differential expression of KRs between the no-HS and HS groups in TLE patients. Bottom of each plot showing the adjusted P-values of genes derived from DESeq2 analysis. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

extracellular protein that covalently interact with latent TGF-β and regulate its function, and may contribute to the malignant phenotype of glioma cells [58]. TGIF1 is expressed in the fetal and adult nervous system, and its mutation and deletion have been related to neurological diseases like holoprosencephaly [59]. Though a recent case report of a child with atypical absence seizures suggests that TGIF1 deletion may have a link to the development of epilepsy, further studies are still needed to confirm this result [60]. Another study using DNA methylation analysis mentioned that TGIF1 was upregulated in a rat model of chronic epilepsy [61]. In summary, our work provides a landscape of the gene network dynamics underlying epilepsy progression and highlighted candidate regulators controlling epileptogenesis that may warrant further investigation as potential anti-epileptogenic targets. Supplementary data to this article can be found online at https:// doi.org/10.1016/j.ygeno.2019.09.020.

mutation of the GAT-1 (SLC6A1), for example, can cause the genetic disorder “epileptic encephalopathy” [54], which is in accord with the decreased activity of SLC6A1 inferred by our analysis. The same applies to AKT3, whose deletion can also cause epileptiform phenotypes [55]. Finally, the association of key regulators with seizure frequency and hippocampal sclerosis in human TLE further proved their validity and clinical relevancy. We identified 10 epileptogenesis-related regulators whose modulation may contribute to both seizure frequency and hippocampal sclerosis (Table 2). Three of them, i.e., FAM107A, LAMB2, and LTBP1, have not been related to epilepsy previously, and one (TGIF1) with limited evidence for epilepsy association. FAM107A, also known as down-regulated in renal cell carcinoma 1 (DRR1), is expressed in the nervous system during embryogenesis, and is downregulated during neuroblastoma carcinogenesis [56]. LAMB2, the β2 chain of laminins, is a critical cortical basement membrane component. Ablation of LAMB2 gene can disrupt cortical lamination and produce dysplasia [57]. LTBP1 (Latent TGF-β binding protein 1) is an 11

Genomics xxx (xxxx) xxx–xxx

Acknowledgments [48] [63,64] [65] [56] [66,67] [57] [58] [46] [68–70] [60,61]

References

Y. Fu, et al.

This work was supported by the by the National Natural Science Foundation of China (NO. U1603285 and NO. 81803960).

glucose hypometabolism; enhance the intrinsic excitability of neurons present in the microglia; neuroinflammation axonal and dendritic outgrowth; synaptic rearrangement nervous system development neuroprotective; related to many neurological diseases critical for cortical laminar organization extracellular regulator of TGF-β activity regulation of neuron survival, proliferation, and differentiation astrocyte-secreted protein; a negative regulator of synapse formation a transcription factor; regulating neuronal development

Yes Yes Yes No Yes No No Yes Yes Uncertain

[1] S.L. Moshe, E. Perucca, P. Ryvlin, T. Tomson, Epilepsy: new advances, Lancet 385 (2015) 884–898. [2] A. Pitkanen, K. Lukasiuk, Mechanisms of epileptogenesis and potential treatment targets, Lancet Neurol. 10 (2011) 173–186. [3] E.M. Goldberg, D.A. Coulter, Mechanisms of epileptogenesis: a convergence on neural circuit dysfunction, Nat. Rev. Neurosci. 14 (2013) 337–349. [4] S.T. Herman, Epilepsy after brain insult: targeting epileptogenesis, Neurology 59 (2002) S21–S26. [5] J. Maguire, Epileptogenesis: more than just the latent period, Epilepsy Curr. 16 (2016) 31–33. [6] R. Citraro, A. Leo, A. Constanti, E. Russo, G. De Sarro, mTOR pathway inhibition as a new therapeutic strategy in epilepsy and epileptogenesis, Pharmacol. Res. 107 (2016) 333–343. [7] C. Heinrich, S. Lahteinen, F. Suzuki, L. Anne-Marie, S. Huber, U. Haussler, C. Haas, Y. Larmet, E. Castren, A. Depaulis, Increase in BDNF-mediated TrkB signaling promotes epileptogenesis in a mouse model of mesial temporal lobe epilepsy, Neurobiol. Dis. 42 (2011) 35–47. [8] S. McClelland, G.P. Brennan, C. Dube, S. Rajpara, S. Iyer, C. Richichi, C. Bernard, T.Z. Baram, The transcription factor NRSF contributes to epileptogenesis by selective repression of a subset of target genes, Elife 3 (2014) e01267. [9] M.R. Johnson, J. Behmoaras, L. Bottolo, M.L. Krishnan, K. Pernhorst, P.L.M. Santoscoy, T. Rossetti, D. Speed, P.K. Srivastava, M. Chadeau-Hyam, et al., Systems genetics identifies Sestrin 3 as a regulator of a proconvulsant gene network in human epileptic hippocampus, Nat. Commun. 6 (2015) 6031. [10] P.K. Srivastava, J. van Eyll, P. Godard, M. Mazzuferi, A. Delahaye-Duriez, J.V. Steenwinckel, P. Gressens, B. Danis, C. Vandenplas, P. Foerch, et al., A systemslevel framework for drug discovery identifies Csf1R as an anti-epileptic drug target, Nat. Commun. 9 (2018) 3561. [11] M. Thom, Review: hippocampal sclerosis in epilepsy: a neuropathology review, Neuropathol. Appl. Neurobiol. 40 (2014) 520–543. [12] M. Mazzuferi, G. Kumar, C. Rospo, R.M. Kaminski, Rapid epileptogenesis in the mouse pilocarpine model: video-EEG, pharmacokinetic and histopathological characterization, Exp. Neurol. 238 (2012) 156–167. [13] O.K. Okamoto, L. Janjoppi, F.M. Bonone, A.P. Pansani, A.V. da Silva, F.A. Scorza, E.A. Cavalheiro, Whole transcriptome analysis of the hippocampus: toward a molecular portrait of epileptogenesis, BMC Genomics 11 (2010) 230. [14] A.M. Bot, K.J. Debski, K. Lukasiuk, Alterations in miRNA levels in the dentate gyrus in epileptic rats, PLoS One 8 (2013) e76051. [15] G. Kalozoumi, O. Kel-Margoulis, E. Vafiadaki, D. Greenberg, H. Bernard, H. Soreq, A. Depaulis, D. Sanoudou, Glial responses during epileptogenesis in Mus musculus point to potential therapeutic targets, PLoS One 13 (2018) e0201742. [16] K.D. Winden, S.L. Karsten, A. Bragin, L.C. Kudo, L. Gehman, J. Ruidera, D.H. Geschwind, J. Engel Jr., A systems level, functional genomics analysis of chronic epilepsy, PLoS One 6 (2011) e20763. [17] E. Ravasz, A.L. Somera, D.A. Mongru, Z.N. Oltvai, A.L. Barabasi, Hierarchical organization of modularity in metabolic networks, Science 297 (2002) 1551–1555. [18] B. Zhang, S. Horvath, A general framework for weighted gene co-expression network analysis, Stat. Appl. Genet. Mol. Biol. 4 (2005) (Article 17). [19] K. Basso, A.A. Margolin, G. Stolovitzky, U. Klein, R. Dalla-Favera, A. Califano, Reverse engineering of regulatory networks in human B cells, Nat. Genet. 37 (2005) 382–390. [20] P. Langfelder, S. Horvath, WGCNA: an R package for weighted correlation network analysis, BMC Bioinforma. 9 (2008) 559. [21] S. van Dam, U. Vosa, A. van der Graaf, L. Franke, J.P. de Magalhaes, Gene coexpression analysis for functional classification and gene-disease predictions, Brief. Bioinform. 19 (2018) 575–592. [22] A.A. Margolin, I. Nemenman, K. Basso, C. Wiggins, G. Stolovitzky, R. Dalla Favera, A. Califano, ARACNE: an algorithm for the reconstruction of gene regulatory networks in a mammalian cellular context, BMC Bioinforma. 7 (Suppl. 1) (2006) S7. [23] M.J. Alvarez, Y. Shen, F.M. Giorgi, A. Lachmann, B.B. Ding, B.H. Ye, A. Califano, Functional characterization of somatic mutations in cancer using network-based inference of protein activity, Nat. Genet. 48 (2016) 838–847. [24] A. Kauffmann, R. Gentleman, W. Huber, arrayQualityMetrics–a bioconductor package for quality assessment of microarray data, Bioinformatics 25 (2009) 415–416. [25] K.Y. Yeung, W.L. Ruzzo, Principal component analysis for clustering gene expression data, Bioinformatics 17 (2001) 763–774. [26] M.E. Ritchie, B. Phipson, D. Wu, Y. Hu, C.W. Law, W. Shi, G.K. Smyth, Limma powers differential expression analyses for RNA-sequencing and microarray studies, Nucleic Acids Res. e43 (2015) e47. [27] F. Del Carratore, A. Jankevics, R. Eisinga, T. Heskes, F. Hong, R. Breitling, RankProd 2.0: a refactored bioconductor package for detecting differentially expressed features in molecular profiling datasets, Bioinformatics 33 (2017) 2774–2775. [28] P. Bardou, J. Mariette, F. Escudie, C. Djemiel, C. Klopp, jvenn: an interactive Venn diagram viewer, BMC Bioinforma. 15 (2014) 293. [29] J. Hardin, A. Mitani, L. Hicks, B. VanKoten, A robust measure of correlation between two genes on a microarray, BMC Bioinforma. 8 (2007) 220.

Acid-sensing ion channel 2 Complement C1q A chain Calponin 3 Family with sequence similarity 107 member A Heat shock protein family B (small) member 1 Laminin subunit beta 2 Latent transforming growth factor beta binding protein 1 Neurotrophic receptor tyrosine kinase 2 Secreted protein acidic and cysteine rich TGFB induced factor homeobox 1 ASIC2 C1QA CNN3 FAM107A HSPB1 LAMB2 LTBP1 NTRK2 SPARC TGIF1

M3 M1 M1 M1 M1 M1 M1 M1 M1 M1

Protein names Key regulators

Table 2 Key regulators associated with both seizure frequency and hippocampal sclerosis.

Modules

Related functions

Known epilepsy association

References

12

Genomics xxx (xxxx) xxx–xxx

Y. Fu, et al.

[30] A. Zeisel, A.B. Munoz-Manchado, S. Codeluppi, P. Lonnerberg, G. La Manno, A. Jureus, S. Marques, H. Munguba, L. He, C. Betsholtz, et al., Brain structure. Cell types in the mouse cortex and hippocampus revealed by single-cell RNA-seq, Science 347 (2015) 1138–1142. [31] Y. Zhou, B. Zhou, L. Pache, M. Chang, A.H. Khodabakhshi, O. Tanaseichuk, C. Benner, S.K. Chanda, Metascape provides a biologist-oriented resource for the analysis of systems-level datasets, Nat. Commun. 10 (2019) 1523. [32] A. Sergushichev, An algorithm for fast preranked gene set enrichment analysis using cumulative statistic calculation, Bio. Rxiv (2016) 60012. [33] S. Durinck, P.T. Spellman, E. Birney, W. Huber, Mapping identifiers for the integration of genomic datasets with the R/Bioconductor package biomaRt, Nat. Protoc. 4 (2009) 1184–1191. [34] A. Lachmann, F.M. Giorgi, G. Lopez, A. Califano, ARACNe-AP: gene network reverse engineering through adaptive partitioning inference of mutual information, Bioinformatics 32 (2016) 2233–2235. [35] da W. Huang, B.T. Sherman, R.A. Lempicki, Systematic and integrative analysis of large gene lists using DAVID bioinformatics resources, Nat. Protoc. 4 (2009) 44–57. [36] M. Kanehisa, S. Goto, Y. Sato, M. Furumichi, M. Tanabe, KEGG for integration and interpretation of large-scale molecular data sets, Nucleic Acids Res. 40 (2012) D109–D114. [37] M.D. Robinson, D.J. McCarthy, G.K. Smyth, edgeR: a Bioconductor package for differential expression analysis of digital gene expression data, Bioinformatics 26 (2010) 139–140. [38] D. Risso, J. Ngai, T.P. Speed, S. Dudoit, Normalization of RNA-seq data using factor analysis of control genes or samples, Nat. Biotechnol. 32 (2014) 896–902. [39] M.I. Love, W. Huber, S. Anders, Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2, Genome Biol. 15 (2014) 550. [40] F. Hong, R. Breitling, C.W. McEntee, B.S. Wittner, J.L. Nemhauser, J. Chory, RankProd: a bioconductor package for detecting differentially expressed genes in meta-analysis, Bioinformatics 22 (2006) 2825–2827. [41] Y. David, L.P. Cacheaux, S. Ivens, E. Lapilover, U. Heinemann, D. Kaufer, A. Friedman, Astrocytic dysfunction in epileptogenesis: consequence of altered potassium and glutamate homeostasis? J. Neurosci. 29 (2009) 10588–10599. [42] Y. Bozzi, M. Dunleavy, D.C. Henshall, Cell signaling underlying epileptic behavior, Front. Behav. Neurosci. 5 (2011) 45. [43] K. Monory, F. Massa, M. Egertova, M. Eder, H. Blaudzun, R. Westenbroek, W. Kelsch, W. Jacob, R. Marsch, M. Ekker, et al., The endocannabinoid system controls key epileptogenic circuits in the hippocampus, Neuron 51 (2006) 455–466. [44] Y. Fukata, M. Fukata, Epilepsy and synaptic proteins, Curr. Opin. Neurobiol. 45 (2017) 1–8. [45] B. Gu, Y.Z. Huang, X.P. He, R.B. Joshi, W. Jang, J.O. McNamara, A peptide uncoupling BDNF receptor TrkB from phospholipase Cgamma1 prevents epilepsy induced by status epilepticus, Neuron 88 (2015) 484–491. [46] G. Liu, B. Gu, X.P. He, R.B. Joshi, H.D. Wackerle, R.M. Rodriguiz, W.C. Wetsel, J.O. McNamara, Transient inhibition of TrkB kinase after status epilepticus prevents development of temporal lobe epilepsy, Neuron 79 (2013) 31–38. [47] M.C. Walker, Hippocampal sclerosis: causes and prevention, Semin. Neurol. 35 (2015) 193–200. [48] H. Zhang, G. Gao, Y. Zhang, Y. Sun, H. Li, S. Dong, W. Ma, B. Liu, W. Wang, H. Wu, H. Zhang, Glucose deficiency elevates acid-sensing Ion Channel 2a expression and increases seizure susceptibility in temporal lobe epilepsy, Sci. Rep. 7 (2017) 5870. [49] S. Watanabe, S. Yamamori, S. Otsuka, M. Saito, E. Suzuki, M. Kataoka, H. Miyaoka, M. Takahashi, Epileptogenesis and epileptic maturation in phosphorylation sitespecific SNAP-25 mutant mice, Epilepsy Res. 115 (2015) 30–44. [50] E. Guney, J. Menche, M. Vidal, A.L. Barabasi, Network-based in silico drug efficacy screening, Nat. Commun. 7 (2016) 10331. [51] C. Lefebvre, P. Rajbhandari, M.J. Alvarez, P. Bandaru, W.K. Lim, M. Sato, K. Wang, P. Sumazin, M. Kustagi, B.C. Bisikirska, et al., A human B-cell interactome identifies MYB and FOXM1 as master regulators of proliferation in germinal centers, Mol. Syst. Biol. 6 (2010) 377. [52] Z. Peng, C.S. Huang, B.M. Stell, I. Mody, C.R. Houser, Altered expression of the delta

[53]

[54]

[55]

[56]

[57]

[58]

[59]

[60]

[61]

[62]

[63]

[64]

[65]

[66]

[67] [68]

[69]

[70]

13

subunit of the GABAA receptor in a mouse model of temporal lobe epilepsy, J. Neurosci. 24 (2004) 8629–8639. M.W. Lopes, F.M. Soares, N. de Mello, J.C. Nunes, F.M. de Cordova, R. Walz, R.B. Leal, Time-dependent modulation of mitogen activated protein kinases and AKT in rat hippocampus and cortex in the pilocarpine model of epilepsy, Neurochem. Res. 37 (2012) 1868–1878. G.L. Carvill, J.M. McMahon, A. Schneider, M. Zemel, C.T. Myers, J. Saykally, J. Nguyen, A. Robbiano, F. Zara, N. Specchio, et al., Mutations in the GABA transporter SLC6A1 cause epilepsy with myoclonic-atonic seizures, Am. J. Hum. Genet. 96 (2015) 808–815. A. Coppola, E. Cellini, H. Stamberger, E. Saarentaus, V. Cetica, D. Lal, T. Djemie, M. Bartnik-Glaska, B. Ceulemans, J. Helen Cross, et al., Diagnostic implications of genetic copy number variation in epilepsy plus, Epilepsia 60 (2019) 689–706. Y. Asano, S. Kishida, P. Mu, K. Sakamoto, T. Murohara, K. Kadomatsu, DRR1 is expressed in the developing nervous system and downregulated during neuroblastoma carcinogenesis, Biochem. Biophys. Res. Commun. 394 (2010) 829–835. S. Radner, C. Banos, G. Bachay, Y.N. Li, D.D. Hunter, W.J. Brunken, K.T. Yee, Beta2 and gamma3 laminins are critical cortical basement membrane components: ablation of Lamb2 and Lamc3 genes disrupts cortical lamination and produces dysplasia, Dev. Neurobiol. 73 (2013) 209–229. I. Tritschler, D. Gramatzki, D. Capper, M. Mittelbronn, R. Meyermann, J. Saharinen, W. Wick, J. Keski-Oja, M. Weller, Modulation of TGF-beta activity by latent TGFbeta-binding protein 1 in human malignant glioma cells, Int. J. Cancer 125 (2009) 530–540. K. Taniguchi, A.E. Anderson, A.E. Sutherland, D. Wotton, Loss of Tgif function causes holoprosencephaly by disrupting the SHH signaling pathway, PLoS Genet. 8 (2012) e1002524. A. Verrotti, C. Palka, G. Prezioso, M. Alfonsi, G. Calabrese, G. Palka, F. Chiarelli, Deletion 18p11.32p11.31 in a child with global developmental delay and atypical, drug-resistant absence seizures, Cytogenet. Genome. Res. 146 (2015) 115–119. K. Kobow, A. Kaspi, K.N. Harikrishnan, K. Kiese, M. Ziemann, I. Khurana, I. Fritzsche, J. Hauke, E. Hahnen, R. Coras, et al., Deep sequencing reveals increased DNA methylation in chronic rat epilepsy, Acta Neuropathol. 126 (2013) 741–756. N.G. Griffin, Y. Wang, C.M. Hulette, M. Halvorsen, K.D. Cronin, N.M. Walley, M.M. Haglund, R.A. Radtke, J.H. Skene, S.R. Sinha, E.L. Heinzen, Differential gene expression in dentate granule cells in mesial temporal lobe epilepsy with and without hippocampal sclerosis, Epilepsia 57 (2016) 376–385. M.I. Fonseca, S.H. Chu, M.X. Hernandez, M.J. Fang, L. Modarresi, P. Selvan, G.R. MacGregor, A.J. Tenner, Cell-specific deletion of C1qa identifies microglia as the dominant source of C1q in mouse brain, J. Neuroinflammation 14 (2017) 48. Y. Chu, X. Jin, I. Parada, A. Pesic, B. Stevens, B. Barres, D.A. Prince, Enhanced synaptic connectivity and epilepsy in C1q knockout mice, Proc. Natl. Acad. Sci. U. S. A. 107 (2010) 7975–7980. Y. Han, H. Yin, Y. Xu, Q. Zhu, J. Luo, X. Wang, G. Chen, Increased expression of calponin-3 in epileptic patients and experimental rats, Exp. Neurol. 233 (2012) 430–437. H.J. Bidmon, B. Gorg, N. Palomero-Gallagher, F. Behne, R. Lahl, H.W. Pannek, E.J. Speckmann, K. Zilles, Heat shock protein-27 is upregulated in the temporal cortex of patients with epilepsy, Epilepsia 45 (2004) 1549–1559. S.E. Brownell, R.A. Becker, L. Steinman, The protective and therapeutic function of small heat shock proteins in neurological diseases, Front. Immunol. 3 (2012) 74. E.V. Jones, Y. Bernardinelli, J.G. Zarruk, S. Chierzi, K.K. Murai, SPARC and GluA1containing AMPA receptors promote neuronal health following CNS injury, Front. Cell. Neurosci. 12 (2018) 22. H. Kucukdereli, N.J. Allen, A.T. Lee, A. Feng, M.I. Ozlu, L.M. Conatser, C. Chakraborty, G. Workman, M. Weaver, E.H. Sage, et al., Control of excitatory CNS synaptogenesis by astrocyte-secreted proteins Hevin and SPARC, Proc. Natl. Acad. Sci. U. S. A. 108 (2011) E440–E449. K.L. van Gassen, M. de Wit, M.J. Koerkamp, M.G. Rensen, P.C. van Rijen, F.C. Holstege, D. Lindhout, P.N. de Graan, Possible role of the innate immunity in temporal lobe epilepsy, Epilepsia 49 (2008) 1055–1065.