Transcriptomic and genetic analyses reveal potential causal drivers for intractable partial epilepsy

Abstract

Mesial temporal lobe epilepsy with hippocampal sclerosis represents the most common epilepsy syndrome in adult patients with medically intractable partial epilepsy. Mesial temporal lobe epilepsy is usually regarded as a polygenic and complex disorder, still poorly understood but probably caused and perpetuated by dysregulation of numerous biological networks and cellular functions. The study of gene expression changes by single nucleotide polymorphisms in regulatory elements (expression quantitative trait loci, eQTLs) has been shown to be a powerful complementary approach to the detection and understanding of risk loci by genome-wide association studies. We performed a whole (gene and exon-level) transcriptome analysis on cortical tissue samples (Brodmann areas 20 and 21) from 86 patients with mesial temporal lobe epilepsy with hippocampal sclerosis and 75 neurologically healthy controls. Genome-wide genotyping data from the same individuals (patients and controls) were analysed and paired with the transcriptome data. We report potential epilepsy-risk eQTLs, some of which are specific to tissue from patients with mesial temporal lobe epilepsy with hippocampal sclerosis. We also found large transcriptional and splicing deregulation in mesial temporal lobe epilepsy with hippocampal sclerosis tissue as well as gene networks involving neuronal and glial mechanisms that provide new insights into the cause and maintenance of the seizures. These data (available via the ‘Seizubraineac’ web-tool resource, www.seizubraineac.org) will facilitate the identification of new therapeutic targets and biomarkers as well as genetic risk variants that could influence epilepsy and pharmacoresistance.

Introduction

Epilepsy is a common neurological disorder that affects people of all ages. Approximately 0.5–1% of people have epilepsy, and another 1.5–5% of people have a seizure at some point in their lives (Banerjee and Hauser, 2007). Mesial temporal lobe epilepsy with hippocampal sclerosis (mTLE+HS) represents the most common epilepsy syndrome in adult patients with medically intractable partial epilepsy (Wieser and Epilepsy, 2004).

Significant advances have been made over the last few years in the field of epilepsy; however, the pathological mechanisms underlying mTLE+HS are still poorly understood and the cause remains unknown. Mutations in the sodium channels SCN1A and SCN1B have been described in a minority of familial cases with temporal lobe epilepsy (TLE) (Abou-Khalil et al., 2001; Scheffer et al., 2007) or familial febrile seizures (Mantegazza et al., 2005). Our group identified a common variation in and near SCN1A that may increase risk for mTLE+HS with childhood febrile seizures in a large genome-wide association study (GWAS) with 1018 subjects with mTLE+HS (Kasperaviciute et al., 2013). Identification of other risk variants can be expected with the increase of sample sizes and collaborative efforts as already demonstrated in other multifactorial disorders (Lambert et al., 2013; Nalls et al., 2014). With this purpose, at least two large-scale collaborative GWAS of epilepsy and/or pharmacoresistance have been funded, EPIPGX (http://www.epipgx.eu/) and Epi4K (Epi, 2012). Recently Epi4K reported ultra-rare genetic variations for familial non-acquired focal epilepsy in known epilepsy genes such as DEPDC5, LGI1, SCN1A and GRIN2A. However, this study was unable to identify GWA significant genes in sporadic focal epilepsy and only LFI1 achieved an uncorrected significant P-value (Epi and Epilepsy Phenome/Genome, 2017).

Microarray studies in temporal surgical specimens in patients with TLE have provided a novel and sophisticated approach to better understand the molecular mechanisms underlying TLE pathology; however, no studies currently exist addressing genome-wide alternative splicing in human mTLE+HS. Given that defects in splicing are being increasingly implicated in various neurological disorders (Licatalosi and Darnell, 2006; Chabot and Shkreta, 2016; Li et al., 2016) and that alternative splicing has been shown to expand gene function in numerous biological processes including synaptogenesis (Li et al., 2007), progress can be enhanced by including the study of transcript variants.

Furthermore, gene expression can be modified by polymorphisms in regulatory elements. The study of gene expression changes that correlate with genetic variation, namely expression quantitative trait loci (eQTL), has been shown to be a powerful and complementary approach to the detection of novel risk loci (Cookson et al., 2009; Ertekin-Taner, 2011; Piasecka et al., 2018), especially when performed in disease-relevant tissue.

No studies currently exist addressing genome-wide alternative splicing, together with genotyping-expression pairing in a large group of human brains with mTLE+HS. With this background, and through a multifaceted and comprehensive approach, we embarked on a series of experiments with the purpose of: (i) providing novel information regarding genes or transcripts as well as pathways, that are affected in mTLE+HS; (ii) identifying genetic loci that can confer risk for mTLE+HS; and (iii) generating a data resource for researchers working in epilepsy or disorders that include epilepsy as part of the phenotype (>200 disorders).

Consequently, a whole (gene and exon-level) transcriptome analysis was performed on 85 neuropathologically confirmed mTLE+HS cortical tissue samples and 75 neurologically healthy controls. Genome-wide genotyping data from the same individuals (patients and controls), and the same cortical tissue samples [Brodmann areas (BA) 20 and 21] were analysed.

Data and results are publicly available as a web-based, searchable database—www.seizubraineac.org—with interactive access to gene and exon expression patterns and genetic variants that influence transcript levels in brain tissue from patients with mTLS+HS and neurologically healthy controls.

We believe our resource will help to prioritize fine-mapping efforts and provide a deeper understanding of the biology of seizure susceptibility.

Materials and methods

Brain tissue samples from 75 neurologically healthy controls, without history or signs of neurological or neuropsychiatric illness or drug abuse, were collected from the Medical Research Council Sudden Death Brain and Tissue Bank in Edinburgh, UK (Ramasamy et al., 2014). A consultant neuropathologist confirmed the absence of pathology in the control group by histological examination of paraffin-embedded brain tissue sections.

Surgical specimens were obtained from patients with chronic pharmacoresistant mTLE+HS who had previously undergone surgical resection at the National Hospital for Neurology and Neurosurgery (NHNN). Samples were assessed and confirmed by a specialist neuropathologist.

Informed and written consent were obtained from patients and controls for all procedures as approved by the institutional review board.

RNA isolation and processing

RNA was isolated from the middle temporal cortex using the RNeasy® Mini Kit (Qiagen). All samples were randomly allocated and hybridized to Affymetrix Human Exon 1.0 ST. Exon Array data were preprocessed using Affymetrix Power Tools (APT), Partek Genomics Suite (Partek Incorporated, St. Louis, MO, USA) and R software.

Genotyping

Samples with MTLE+HS (63%) were genotyped on Illumina HumanHap550 chip (Kasperaviciute et al., 2013), the Human610-Quadv1 and the Illumina Infinium Omni1-Quad BeadChip. Controls were genotyped on the Illumina Infinium Omni1-Quad BeadChip.

All the analyses were residual-adjusted and unless specified we report only those associations with false discovery rate (FDR) corrected significance (Benjamini and Hochberg method) at 0.05. For further details see the Supplementary material.

Data availability

Processed expression and eQTL data described in this study can be downloaded and queried through www.seizubraineac.org portal. The raw expression data (microarray CEL files) are available on Gene Expression Omnibus (GEO) with the accession code GSE46706. Genotype data are accessible via Ramasamy et al. (2014) and Kasperaviciute et al. (2013). The source code that supports the findings of this study is available from the corresponding author, upon reasonable request.

Results

Genome-wide differential expression analyses

Sample information is listed in Table 1. After correcting for batch, age, gender and RNA integrity number (RIN) differences, 6103 transcripts (5523 genes) showed significant differential expression at a 5% FDR, with 3016 and 3087 transcripts showing significantly higher and lower expression, respectively in the mTLE+HS cortex when compared to controls (Supplementary Table 1). Table 2 shows the top 35 differentially expressed genes (FDR 5%, and fold change ≥1.5).

Gene Ontology (GO) enrichment analyses showed that the top (based on FDR P-values) and most abundant category terms on the list of downregulated genes in the mTLE+HS cortex were mainly related with neuron development, differentiation and synaptic signalling (Supplementary Table 2). Detailed examination of the expression of ion channel genes showed SCN1A, SCN1B and KCNS1 as the genes with the largest differences in expression between patients and controls (Fig. 1). All these genes showed lower expression in patients with mTLE+HS compared to controls (FDR 5%).

Although a large number of mTLE transcriptome studies have been published, all studies with publicly available data have been performed in animal models or on human hippocampal tissue. It has not been possible therefore, to test if our findings could be replicated in a different study. Wang et al. (2010) reviewed 12 papers of genome-wide expression profiling in models of mTLE (animal and human), including one study that compared epileptogenic and non-epileptogenic tissue (based on the intraoperative electrocorticographic recording) in temporal cortex of five patients with mTLE and two studies that compared different brain regions in a rat mTLE model (Wang et al., 2010). The review provided a list of 53 genes showing expression changes in more than two publications. We found that 70% of these genes were contained in our list of significantly differentially expressed genes in mTLE (FDR 5%).

As there were significant differences in age and RIN values between the control and the mTLE+HS group, the analyses were residual-adjusted. We wanted to test if a smaller but more homogeneous sample group would provide the same or similar results. We therefore randomly selected samples with similar ages and RIN values in the patient and control groups. We reduced the sample size to 32 controls and 49 mTLE samples with a mean age and RIN value of 43.7 ± 11.5 and 5.8 ± 1.4, respectively in the control group and 39.41 ± 10.7 and 6.2 ± 1.3, respectively in the mTLE group (P-values > 0.01). After correcting for batch effect, 8738 transcripts showed significant differential expression (FDR 5%). We then checked how many of the 6103 transcripts that showed differential expression in the large mTLE sample dataset were also differentially expressed in the smaller mTLE group when comparing mTLE and controls. We found that 6028 of the 6103 transcripts (99%) showed differential expression with no significant difference in age and RIN values. These analyses provide validation of our data and are supported by the fact that many of the differentially expressed genes identified in this study have previously been reported in the literature.

Weighted gene co-expression networks

To uncover potential regulators underlying mTLE+HS, we next applied Weighted Gene Co-expression Network Analyses (WGCNA; Langfelder and Horvath, 2008). Gene modules attempt to represent groups of genes that act together in a concerted manner and are often enriched in specific pathways or with molecular functions. We first used the ‘module preservation’ method to identify modules that were non-preserved across conditions. Looking for non-preserved or weakly preserved modules potentially allows us to identify pathologically dysregulated pathways or molecular systems that are either acquired or lost in comparison with a healthy network.

First, we constructed separate co-expression networks for mTLE+HS and control groups using the expression values of 16 000 genes with the highest variable expression in each group. This method resulted in 14 modules for the network from healthy samples (‘C’ modules) and 21 modules from mTLE+HS samples (‘TLE’ modules).

We then computed the preservation of modules for the healthy (reference) network as compared to the mTLE+HS (test) network. Using described thresholds (Langfelder et al., 2011), a Z-score <2 indicates no preservation, 2–10 indicates weak to moderate preservation and >10 indicates preservation. All modules had moderate to high preservation. The C-Turquoise module, which was significantly enriched by neuronal pyramidal markers and functionally enriched with genes involved in action potential showed the weakest preservation score (Z-score = 4) (Supplementary Tables 3 and 4).

When the mTLE+HS network was used as the reference network, two modules, TLE-Grey60 and TLE-Purple were not conserved in healthy tissue (Grey60 Z-score = 0.8 and Purple Z-score = 0.9). The module TLE-Grey60 was significantly and exclusively enriched by marker genes for neurons with the strongest enrichment from the neuron marker gene list obtained from a single cell study of epilepsy human temporal cortex (Darmanis et al., 2015) (Fisher’s exact test FDR P = 2.87 × 10⁻¹²) and specifically neuronal community four (FDR P = 4.86 × 10⁻⁵), which is an excitatory type of neurons expressing SLC17A7, NEUROD6 and SATB2. Functionally, the module Grey60 also showed significant over-representation of genes related to synaptic signalling (GO:0099536, FDR P = 5.32 × 10⁻²¹) and pathways related to the neuronal system (Reactome 12687663, FDR P = 4.30 × 10⁻¹⁰).

Functionally, the TLE-Purple module did not show any significant enrichment for any current functional annotation (GO, KEGG or Reactome). Cell type analyses, however, revealed enrichment by an excitatory neuronal subtype identified across six human cortical regions, including temporal cortex and with predominance to be expressed in cortical layer VI (FDR P = 4.53 × 10⁻⁵) (Lake et al., 2016).

Identification of critical network genes

We subsequently attempted to identify hub genes in non-preserved modules, as hub genes are hypothesized to be critical to the proper functioning of the module (Varki et al., 2008) and can thus be potential biomarkers of disease. The identification of hub genes may also allow us to infer the main biological function of the module. Using the 250 gene-gene interactions and defining hub genes as genes that interact with at least 15 other genes, we found a large number of hub genes in the TLE-Grey60 module that have previously been linked to different forms of epilepsy including SCN2A, GABRB2, GABRA1, SCN8A, or ITPR1 (Fig. 2). In addition, upregulation of the ATPase plasma membrane Ca² transporting 1 protein (ATP2B1 or PMCA1) in the inner molecular layer of the dentate gyrus has been reported in early and late stage animal models of TLE (Ketelaars et al., 2004).

Our analyses also revealed other interesting hub genes that have not previously been associated with epilepsy, such as CAP2, MYT1L, CADPS2 or LDB2.

In the TLE-Purple module (Fig. 2), CNNM1 was the gene with the highest number of interactions with other genes within the module. The role of CNNM1 in the brain is not well known, but it has been suggested to be involved in ion transport (Wang et al., 2004). Interestingly, mutations in VPS13A, another hub gene in the TLE-Purple module, cause chorea-acanthocytosis, a rare genetic disorder characterized by hyperkinetic movements, seizures, cognitive decline, neuropsychiatric symptoms and acanthocytosis. Some people with mutations in this gene have seizures as the first and predominant clinical manifestation, with some of the seizures being of temporal lobe origin and difficult to control (Al-Asmi et al., 2005; Benninger et al., 2016).

However, the most remarkable aspect is that the module was significantly enriched by hub genes (MDGA2, POU6F2, SLC8A1, CRTAC1, SEMA3A, FEZF2 and SULF1) involved in the regulation of CNS development, axon guidance, cell movement or neuronal differentiation (Fisher’s exact test P = 2.76 × 10⁻⁵; Fig. 2C). The differential expression and gene-gene relation of this non-preserved module could be associated with abnormal neural circuit development in the temporal cortex with mTLE+HS. Changes in neuronal connectivity as well as abnormal migration and location of hippocampal neurons are known in epilepsy and are not limited to the hippocampus.

To identify co-expressed modules that correlate with the presence of mTLE+HS, and to test if hub genes of the non-preserved modules play an important role in this disorder, we constructed a network of modules of highly co-expressed genes using both mTLE+HS and control samples (‘All’ modules). The expression levels of each module were summarized by the first principal component (the so-called module eigengene) and related to the presence of mTLE+HS. After removing the effect of age, gender and RIN, WGCNA identified 8 of 13 co-expression modules significantly associated with mTLE+HS. All-Turquoise, All-Salmon, and All-Pink modules were under-expressed and enriched for biological processes such as synaptic signalling, heterocycle metabolic processes or transmembrane transport whereas All-Brown, All-Red, All-Green and All-Blue were overexpressed and were enriched for processes such as lipid metabolism, single organism metabolic processes, immune responses and the regulation of cell projection organization (Supplementary Table 3).

The All-Turquoise module was of particular interest as it was mostly composed of genes belonging to the non-preserved mTLE+HS modules TLE-Grey60 and TLE-Purple as well as TLE-Blue module and the less preserved C-Turquoise module from the healthy network. The All-Turquoise module was therefore enriched for GO or KEGG categories such as synaptic signalling (GO:0099536, FDR P = 5.23 × 10⁻¹²), nervous system development (GO:0007399, FDR P = 1.17 × 10⁻¹²) and axon guidance (KEGG:04360, FDR P = 0.0007). Cell type analyses also revealed enrichment of neuronal biomarkers, for pyramidal neurons from corticosensorial areas (FDR P = 5.15 × 10⁻¹⁰), epilepsy human temporal cortex (FDR P = 1.82 × 10⁻⁵) and hippocampus, mainly excitatory (FDR P = 0.001) but also inhibitory subtypes (somatostatin positive interneurons FDR P = 0.01). SCN1A was one of the most highly connected genes together with CADPS2, ATP1A1, CABP1, HIVEP2 and KCNB1 (all of them from the TLE-Blue module and C-Turquoise module) as well as ITPR1, CAP2, ATP2B1, KCNQ5, ROCK2 and MYT1L from the non-conserved TLE-Grey60 module and CNNM1 and MDGA2 from the non-conserved TLE-Purple module.

Overall, these results suggest that these MTLE+HS hub genes and especially ITPR1 and CNNM1, could be driving the disruption of a healthy non-epileptic set of functionally related genes (the C-Turquoise module) and promoting the division and creation of a new group of gene-gene interactions and new gene’s leaders (hub genes; Fig. 2).

As differences in expression do not necessarily cause epilepsy but might result from the epileptic activity per se, we investigated whether any of the modules, and in particular the non-preserved TLE-Grey60 and TLE-Purple modules, were enriched with genes related to epilepsy according to the OMIM database. As a negative control, we also tested another clinical phenotype, which shows abnormalities in the temporal lobe, schizophrenia. The TLE-Grey60 module showed a significant enrichment of OMIM epilepsy genes (P = 6.66 × 10⁻⁶) but no OMIM schizophrenia (SCZ) genes. TLE-Purple module did not show significant enrichment of epilepsy or schizophrenia genes. We also tested for enrichment of de novo mutations from published whole-exome sequencing studies of trio datasets for epileptic encephalopathy, autism disorder, schizophrenia, intellectual disability and developmental disorders (Epi et al., 2013; Euro et al., 2014). Only the TLE-Grey60 module was consistently enriched with de novo mutations and only for epileptic encephalopathy (Fisher’s exact test FDR P = 0.006).

Identification of potential drivers

We first looked for cis-SNPs (single nucleotide polymorphisms) that could affect the expression of individual genes or its exons in the temporal cortex of patients and controls separately. A total of 16 461 genes accounting for 149 184 exons were tested in patients (mTLE+HS-eQTL) and controls (ctrl-eQTL) separately, resulting in 1265 and 1324 cis-eQTLs, respectively, at 5% FDR. Fifty-two per cent of all non-shared eQTLs showed the same direction of effect in both groups and potentially due to insufficient power to detect eQTL signals in one of the groups. Of all non-shared eQTLs, none showed opposite directional effect or were significant eQTLs in both groups, except the SNP rs28673430 (chr7: 73983339) acting on the gene WBSCR27 (probe ID: t3056320). While the G allele increases the expression of the gene in the control group, the same allele decreases the expression in the disease group (Fig. 3). METTL2 (also known as WBSCR27), encodes a protein belonging to the ubiE/CoQ5 methyltransferase family. The biological function of the gene is not known; however, the gene is deleted in Williams syndrome, a multisystem developmental disorder that sometimes presents with seizures and is caused by the deletion of contiguous genes at 7q11.22-q11.23 (Schubert, 2009).

We did not look for enrichment of epilepsy GWAS SNPs within our eQTL results as the existing few focal epilepsy or meta-analysis GWAS have either been unsuccessful or unreproducible. We did, however, look for SNPs that could affect the expression of genes or exons involved in epilepsy according to OMIM or genes containing de novo mutations. Thirty-eight and 41 regulated genes associated with disorders in which epilepsy or seizures are the primary feature, were found to be regulated by common genetic variants on the mTLE+HS and control group, respectively. Thirty-one genes containing de novo mutations were also affected by SNPs (Table 3 and Supplementary Table 5).

We were insufficiently powered to look for trans-eQTLs (SNPs that were at a distance >1 MB from the transcription start site of an associated transcript on the same or different chromosome). However, recent studies have demonstrated that some cis-eQTLs are also trans-eQTLs because the variation in the expression of the cis gene affects the expression of a trans gene or genes (Brynedal et al., 2017). We, therefore, looked for SNPs that could regulate hub genes, as they are potential regulators of the majority of genes in the module. SNPs that can regulate or change the expression of hub genes may also regulate/change the expression of other genes that interact with the hub genes (trans-eQTLs). The only hub-eQTL gene of interest was the SNP rs2241819 that was found to regulate the expression of the hub gene CYP4V2 in the All-Red module. This SNP also correlates with the expression of 12 other genes (trans-eQTLs) within this module (PAX6, SOX2, AGT, HIST2H2BE, TTYH1, SEPT2, MTMR10, SLC44A2, SLC38A3, GATM, FERMT2, FKBP10, and EPHX1) (FDR 5%). The A6-Red module is one of the upregulated co-expression modules significantly associated with mTLE+HS and is enriched mainly by astrocyte markers (Fisher’s exact test FDR P = 6.44 × 10⁻⁴³) but interestingly is also enriched by a cell community identified by Darmanis et al. (2015), which displays characteristics of both neurons and astrocytes, with the suggestion by the authors that it could be a new cell type or bipotent progenitor cell (Fisher’s exact test FDR P = 0.006). Functionally, the A6-Red module is highly enriched by carboxylic acid metabolic processes, drug metabolic process or ion transport as well as CNS, or eye development (FDR 5%).

Mesial temporal lobe epilepsy with hippocampal sclerosis differential exon usage

Given that alternative splicing allows individual genes to produce multiple protein isoforms with distinct properties and functions, differences in exon usage may therefore promote the overexpression or underexpression of certain types of protein isoforms (Blencowe, 2006), resulting in non-desirable effects

After correcting for batch, age, gender, and RIN values, 2618 genes (2627 transcripts) showed differential exon usage (significant differences in expression in at least one exon) when comparing mTLE+HS with control subjects (FDR = 5%). Fifty-seven per cent of these genes also showed significant differential expression (FDR 5%) of the gene in the mTLE+HS cortex when compared to controls. Genes showing differential exon usage in the mTLE+HS cohort showed an enrichment of biological processes such as cell morphogenesis (GO:0000902, FDR P = 1.75 × 10^–10), circulatory system development (GO: 0072358, FDR P = 1.75 × 10^–10) or cell migration (GO:0016477, FDR P = 9.15 × 10⁻⁹).

Using a conservative splicing index algorithm (FDR =5% and a minimum 1.5-fold splice index difference between the mTLE+HS and control group) 124 genes showed differential exon usage with an enrichment of biological processes such as cell adhesion (GO:0007155, P = 9.79 × 10^–10); immune response (GO:0006955, P = 1.40 × 10^–9), or response to drugs (GO:0042493, P = 2.68 × 10^–5).

One of the genes showing differential exon usage is ankyrin 3 (ANK3), which encodes a protein essential for the correct clustering and activation of sodium channels during the development of axonal initial segments, nodes of Ranvier and postsynaptic folds of the neuromuscular junction (Zhou et al., 1998). Mutations and polymorphisms in ANK3 have been associated with a wide spectrum of disorders including a homozygous frameshift mutation affecting the longest isoform of the gene that has been reported to cause epilepsy, in addition to intellectual disability and aggressive behaviour (Iqbal et al., 2013). In our study, there was an enrichment of the transcript variant 1 (NM020987) in controls, which encodes the longer isoform. In the mTLE+HS group, however, there was an enrichment of the transcript variant 2 (NM001149), a transcript considerably shorter than variant 1, which doesn’t include the ankyrin repeats that mediate protein-protein interactions (Fig. 4). The decrease of longer isoforms in mTLE+HS could contribute to a decrease in sodium channel activity and thus predispose to abnormal excitability.

Another gene that displayed different exon use between patients and controls was Neurexin 3 (NRXN3), which encodes a type of trans-synaptic cell adhesion molecule and receptor. The neurexins not only organize and synchronize the assembly of synapses, but also likely determine their properties and diversity. Differential promoter and splicing usage result in the expression of different neurexin isoforms, with transcripts initiating from the upstream promoter encoding alpha isoforms [which contains epidermal growth factor (EGF)-like sequences] and transcripts initiating from the downstream promoter encoding beta isoforms (missing EGF-like sequences). We found that the expression of beta isoforms was higher in patients with mTLE+HS compared with controls (Fig. 4). Although the function of the different neurexin isoforms remains unclear, recent studies have suggested that neurexins perform multiple independent functions that are regulated by distinct events of alternative splicing (Aoto et al., 2013). Supplementary Table 1 and Fig. 4 show genes with differential exon usage in the mTLE+HS group.

Mesial temporal lobe epilepsy with hippocampal sclerosis differential exon usage and cis-eQTLs

Of the 2618 genes showing differential exon usage in the mTLE+HS group compared to control subjects (FDR 5%), the expression of the same exons in 168 genes was regulated by SNPs in the mTLE+HS group (284 exon-level cis-eQTLs, FDR 5%) and the exon expression of 193 genes was impacted by SNPs in the control group (223 cis-eQTLs, FDR 5%) (Supplementary Table 4). One interesting example includes the cis-eQTL involving the gene tubulin tyrosine ligase-like 1 (TTLL1). Although the exact role of the protein is unclear, cellular activities such as neuronal differentiation or ciliary and flagellar motility, as well as synaptic transmission, have been suggested (Ikegami et al., 2007).

Mesial temporal lobe epilepsy with hippocampal sclerosis, febrile seizure and temporal lobe sclerosis

MTLE+HS is frequently preceded by febrile seizures in infancy and early childhood. mTLE+HS is also described in families with epilepsy syndromes in which febrile seizures are a common feature. Furthermore, our mTLE+HS GWAS reported a significant association between a genetic variant in the SCN1A gene and mTLE+HS and with febrile seizures no association was detected for the mTLE+HS group without history of febrile seizures (Kasperaviciute et al., 2013). Looking for potential gene expression differences between both groups, we compared cases with (n = 49) and without (n = 28) known febrile seizure history. To our surprise, there were no differences in gene or exon expression (FDR 5%), even without taking fold changes into account.

This outcome could be the result of studying tissue from the temporal cortex instead of hippocampus. We have focused on surgical tissue from temporal cortex to avoid problems associated with different cellular compositions of the tissue derived from patient or control populations. Both regions play an important role in the pathogenesis of mTLE+HS but the loss of neuronal populations in the hippocampus is a key feature of mTLE+HS, whereas temporal cortex exhibits less dramatic histological changes (Thom, 2009; Thom et al., 2009).

In our cohort, 19% (n = 16) of cases presented with features of temporal lobe sclerosis; this is where neuronal loss from the superficial cortical layer II and gliosis is evident in addition to hippocampal sclerosis (Table 1). A comparison between cases with and without temporal lobe sclerosis produced non-significant differences in gene or exon expression, even on cell-specific markers such as GFAP (P = 0.26), MOG (P = 0.50), or genes involved in neuron death or apoptosis, which support the view that the RNA changes found in our study do not merely reflect differences in tissue composition.

Discussion

Through a whole (gene and exon-level) transcriptome analysis on a large cohort of cortical tissue sample, we have identified a large number of genes that are differently expressed between mTLE tissue and controls.

Downregulated genes in the mTLE+HS cortex were mainly related with neuron development, differentiation and synaptic signalling. SCN1A, SCN1B and KCNS1 were the ion channel genes with the largest differences in expression between patients and controls. SCN1A and SCN1B have been associated with familial cases of mTLE+HS, and SNPs in SCN1A have been shown to increase susceptibility to mTLE+HS with febrile seizure (Kasperaviciute et al., 2013). KCNS1 has not been related with epilepsy until now. Its protein is not functional by itself but modulates the activity of K_v2 alpha subunits (KCNB), acting as an inhibitor of neuronal excitability (Salinas et al., 1997).

Upregulated genes in the mTLE+HS cortex were instead related to immune response and vascular development. Alterations in regional cerebral blood flow and vasoregulation as well as blood–brain barrier breakdown and infiltration of peripheral immune cells have been described in human and experimental epileptic brain tissue (Vezzani et al., 2013; Han et al., 2017). There is increasing evidence of a critical role of immunity and inflammation in the evolution of drug-resistant epilepsy. The disruption of these mechanisms has been shown to occur after the presence of acute seizures but there is also increasing evidence to the contrary (e.g. inflammatory mediators may initiate or trigger early seizures) (Webster et al., 2017) and therefore currently these mechanisms are believed to be involved in the predisposition and the maintenance of epileptogenesis.

We also identified the TLE-Grey60 module that was enriched with genes related to epilepsy according to the OMIM database and de novo mutations for epileptic encephalopathy. Thus this module may be causally involved with the development of epilepsy. Interestingly, two of the genes in the TLE-Grey60 module were LGI1 and SESN3. Mutations in LGI1 result in autosomal dominant lateral temporal epilepsy with auditory features (ADPEAF) (Fanciulli et al., 2012), and ultra-rare variants have been found in familial non-acquired focal epilepsy (Epi and Epilepsy Phenome/Genome, 2017). SESN3, a member of the Sestrin family, has recently been identified as a regulator of a proconvulsive gene co-expression module in the human hippocampus of a group of patients with mTLE (Johnson et al., 2015). Furthermore, our analyses revealed hub genes that haven’t previously been associated with epilepsy such as CAP2, MYT1L, CADPS2 or LDB2. Interestingly, mutations in CAP2 or MYT1L have previously been described in patients with neurological or behavioural abnormalities who present with seizures (Celestino-Soper et al., 2012; Mayo et al., 2015; Blanchet et al., 2017). The CADPS2 gene encodes a member of the calcium-dependent activator of secretion (CAPS) protein family; calcium-binding proteins that regulate the exocytosis of synaptic and dense-core vesicles in neurons and neuroendocrine cells. Mutations on chromosome 7q31 including CADPS2 have been associated with disorders affecting cognitive functions and susceptibility to autism (Okamoto et al., 2011; Sadakata et al., 2013). Finally Lim domain binding (LDB) protein 2 has been suggested to be important for the regulation of cell fate differentiation (Leone et al., 2017) and to play a role in vascular remodelling (Choi et al., 2018).

We also identified a number of genes that are differentially spliced between mTLE tissue and controls. Although the splicing events identified exhibit only small magnitude changes (<2-fold change), even small changes may have dramatic physiological consequences if the protein encoded by the gene plays an important role in a specific cell function. Determination of the role of these splice events requires a more detailed study.

It is currently unclear to what extent misregulation of neuronal exons contributes to epilepsy. Reports of mutations causing abnormal splicing associated with epilepsy are increasing, however, with the development of genomic technologies. For example, splicing mutations in ALDH7A1 or MTHFR have been found in Mendelian diseases characterized by recurrent seizures (Striano et al., 2009; Prasad et al., 2011). Recent studies have also reported how neuronal excitation can be affected by the manipulation of the expression of splicing regulators such as RBFOX1 (Lal et al., 2015).

Aberrant splicing events, however, might also result from the epileptic activity per se. For example, abnormal gephyrin splicing and expression found in the hippocampus of patients with mTLE has been reported to occur in response to the cellular stress arising from seizures (Forstera et al., 2010). We were therefore interested in changes to splicing patterns that were regulated genetically. For example, splicing at exon 5 in human SCN1A is mutually exclusive with the choice of either exons 5A or 5N (for adult and neonatal forms, respectively). Tate et al. (2005) identified a significant association between an SNP that affects alternative splicing of exon 5 and the doses of two anti-epileptic drugs. Patients with the allele that increases transcripts with the neonatal isoform needed less medication than patients with the allele that facilitates the adult isoform (Tate et al., 2005).

We therefore report splicing events that are controlled by SNPs (cis-eQTLS), which exclude those splicing events that might result from the epileptic activity per se. This provides a reduced number of potential genes that contribute to the cause of mTLE or the response to anti-epileptic drugs.

Interestingly we also found that a large proportion of eQTLs were not identifiable in healthy tissue. A recent study looking for cis- and trans-eQTLs in different tissues from patients with coronary artery disease (CAD) reported more CAD-associated risk SNPs identified by GWAS overlapping eQTLs found in tissue from patients with the disease than without (Franzen et al., 2016). These results suggest that the inclusion of tissue from people with and without the disease may increase our chances of identifying the true disease-risk gene. Further studies will be needed to clarify which findings are as a consequence of real biological mechanisms or confounding factors.

We also report an indirect trans-eQTL, the SNP rs2241819, which was found to regulate a module enriched mainly by astrocyte markers and by a cell community identified in a single cell study of epilepsy human temporal cortex (Darmanis et al., 2015) and which displays characteristics of both neurons and astrocytes with the suggestion by the authors that it could be a new cell type or bipotent progenitor cell. Interestingly, some of the genes of the All-Red module are CYP4V2, PAX6, SOX2, AGT, HIST2H2BE, TTYH1, SEPT2, MTMR10, SLC44A2, SLC38A3, GATM, FERMT2, FKBP10, EPHX1 and HEPACAM. Mutations in CYP4V2, as in PAX6 and SOX2, cause eye disorders (Smith et al., 2009; Astuti et al., 2015; Mauri et al., 2015) and some individuals with mutations in PAX6, SOX2, and HEPACAM present with mesial temporal malformation and epilepsy (Sisodiya et al., 2006; Marchese et al., 2014). PAX6 and SOX2 have also been reported to play a crucial role in gliogenic competence (Gomez-Lopez et al., 2011). Microdeletions of chromosome 15q13.3, including MTMR10 and chromosome 2q37, including SEPT2, result in severe developmental delay, epilepsy, and microcephaly (Mulley and Dibbens, 2009; Imitola et al., 2015). Furthermore, SNPs in the EPHX1 gene have been reported to be involved in the adverse effects of anti-epileptic drugs (Hung et al., 2012), and TTYH1 has been implicated in aberrant neuronal structural plasticity in vivo, as increased TTYH1 protein expression was observed in the molecular layer of the dentate gyrus during epileptogenesis (Wiernasz et al., 2014). As brain pathology or seizures themselves can induce astrocytes to become reactive, it is difficult to distinguish whether glial changes found in mTLE+HS tissue could be the cause or consequence of epilepsy. However, we cannot ignore the fact that glial cells play a role in synaptic development/function and recently, several lines of evidence indicate their contribution as a cause of epilepsy (Bedner et al., 2015; Coulter and Steinhauser, 2015) and other neurological disorders (Heneka et al., 2015). Further studies will be required to clarify the specific contribution of this module, but the identification of one SNP regulating the most relevant functional genes of this network, indicates a potential genetic mechanism underlying these glial changes.

One limitation of this study is the different source of tissue between patients and controls. Ideally, one would use tissue from one source only. However, it is difficult to obtain large numbers of autopsy tissue samples from patients with a history of epilepsy, and who died of non-epilepsy related causes. The same can also be said regarding surgical tissue, particularly when trying to source temporal cortex from living patients without a neurological disease. Several studies have compared the transcriptome of human tissues from autopsy and surgery and have showed differential expression only in a minority of genes (1–3%) mostly involved in basal cell metabolism and hypoxia (Sluimer et al., 2007). Furthermore, an increasing number of expression profiling studies have shown how well gene expression measured in post-mortem samples reflects events happening in vivo.

Other related shortcomings have been the difference in age between patients and controls, as surgical tissue is frequently obtained from young people with epilepsy whereas most post-mortem tissue is taken from individuals who died at an older age. In relation to these two factors, the quality of the RNA is normally better in surgical than in post-mortem tissue and hence the differences. To overcome these obstacles we have adjusted and corrected for differences in age and RIN, and we have tested if using a smaller but more homogeneous group, formed by samples without differences in age and RIN, would show the same or similar results. Furthermore, some analyses have been done separately for patients and controls (e.g. eQTLs) and although some of these changes could be a product of the variability between sample groups, the fact that many of the genes showing differential expression have already been reported to be associated with epilepsy in human and animal studies shows the validity of our study.

Finally, despite the advantages in examining transcripts directly from the temporal cortex instead of the hippocampus to avoid problems with changes in cellular composition, it is still expected that some transcripts specific to one cell type will remain below the level of detection. Future transcriptional studies on single cells will better characterize the transcriptional landscape of the different brain regions and their component cells and will provide greater insight into mTLE. However, the study of hundreds of cells in a large number of individuals, as required in eQTLs studies, is a very expensive and time-consuming task. One alternative is to investigate if the common genetic variants of interest influence gene expression in the temporal cortex using our dataset and then test if the expression of the gene shows a higher preference for one region or cell type, using existing single cell studies. (Darmanis et al., 2015; Lake et al., 2016; Krishnaswami et al., 2016). Computational methods like population-specific expression analysis together with some validation using LCM experiments could also be another alternative (Kuhn et al., 2012).

Conclusion

Our study is the first to provide: (i) epilepsy-relevant gene networks that implicate neuronal and glial mechanisms, including one module enriched with genes that have been related to epilepsy according to the OMIM and are therefore likely to have a causal role in mTLE+HS; (ii) mTLE+HS-associated splicing changes in ion and non-ion genes; (iii) different genetic loci that affect the expression of genes and/or transcripts that have been involved directly or indirectly in epilepsy and are therefore strong candidates as genetic risk factors for mTLE+HS; and (iv) a data resource for researchers working on epilepsy or disorders that include epilepsy as part of the phenotype, that could help to speed the discovery of new genetic loci and facilitate the understanding of how this locus contributes to disease. Our findings indicate that information obtained from annotating the functional effects of genetic risk variants in tissue from individuals with the disease are at least complementary to those obtained from corresponding tissues isolated from healthy individuals and therefore essential to fully understand how common genetic risk variants influence complex traits such as epilepsy.

In summary, ‘Seizubraineac’ brings rapid insights into mTLE+HS and will be a complementary resource for exploring GWAS findings moving forward.

Abbreviations

Abbreviations
 
• eQTLs

  expression quantitative trait loci

 
• GWAS

  genome-wide association studies

 
• HS

  hippocampal sclerosis

 
• mTLE

  mesial temporal lobe epilepsy

 
• RIN

  RNA integrity number

 
• SNP

  single nucleotide polymorphism

Acknowledgements

We are grateful to the affected individuals and their families who made this study possible. We acknowledge Epilepsy Research UK (ERUK) for its support. We thank Prof Sanjay Sisodiya for providing us with the mTLE+HS GWAS genotype and Ms Andrea Andreja Avbersek and Dr Claudia Catarino for their help in the mTLE clinical phenotyping.

Funding

Funding was from a Marie Curie International Reintegration Grant FP7-PEOPLE-2009-RG grant No 256545 and an Epilepsy Research UK Fellowship (F1206) to M.M. S.G. was supported by Alzheimer’s Research UK through the award of a PhD Fellowship (ARUK-PhD2014–16). M.P. holds a fellowship from the National Atomic Energy Commission, Argentina.

Competing interests

M.E.W. is an employee of Genomics plc, a company providing genomic analysis services to the pharmaceutical and healthcare sectors. The remaining authors report no competing interests.

References