Neuronal and glial DNA methylation and gene expression changes in early epileptogenesis

Background and aims Mesial Temporal Lobe Epilepsy is characterized by progressive changes of both neurons and glia, also referred to as epileptogenesis. No curative treatment options, apart from surgery, are available. DNA methylation (DNAm) is a potential upstream mechanism in epileptogenesis and may serve as a novel therapeutic target. To our knowledge, this is the first study to investigate epilepsy-related DNAm, gene expression (GE) and their relationship, in neurons and glia. Methods We used the intracortical kainic acid injection model to elicit status epilepticus. At 24 hours post injection, hippocampi from eight kainic acid- (KA) and eight saline-injected (SH) mice were extracted and shock frozen. Separation into neurons and glial nuclei was performed by flow cytometry. Changes in DNAm and gene expression were measured with reduced representation bisulfite sequencing (RRBS) and mRNA-sequencing (mRNAseq). Statistical analyses were performed in R with the edgeR package. Results We observed fulminant DNAm- and GE changes in both neurons and glia at 24 hours after initiation of status epilepticus. The vast majority of these changes were specific for either neurons or glia. At several epilepsy-related genes, like HDAC11, SPP1, GAL, DRD1 and SV2C, significant differential methylation and differential gene expression coincided. Conclusion We found neuron- and glia-specific changes in DNAm and gene expression in early epileptogenesis. We detected single genetic loci in several epilepsy-related genes, where DNAm and GE changes coincide, worth further investigation. Further, our results may serve as an information source for neuronal and glial alterations in both DNAm and GE in early epileptogenesis.


Introduction
Epilepsy is defined as an inherent predisposition of the brain to recurrently generate epileptic seizures [1] and affects an estimated 65 million people world-wide [2]. Temporal lobe epilepsy (TLE) is the most common subtype amongst the focal epilepsies [3], with hippocampal sclerosis being detected in 70% of drug resistant TLE patients [4,5]. The typical clinical course of the sub-entity, mesial temporal lobe epilepsy with hippocampal sclerosis (mTLE-HS), is characterized by an initial precipitating event (e.g. cerebral trauma, inflammation, prolonged febrile seizure), a seizure-free latency period and finally the onset of spontaneous and progressive seizures [6]. This metamorphosis into a brain prone to spontaneous recurrent seizures of progressive nature, also referred to as epileptogenesis [7,8], is characterized by a plethora of cellular and molecular changes in both neurons and glia [5,[9][10][11][12][13][14][15][16][17][18].
One third of people with epilepsy respond inadequately to treatment with the primarily symptom-alleviating antiepileptic drugs of today [19], rendering the identification of potential upstream effectors of epileptogenesis and the development of disease modifying antiepileptic drugs a task of upmost importance [20,21].
DNAm, in the context of this paper, the methylation of CpG nucleotides in the DNA [22], plays a primordial role in brain development, cell fate, tissue specific gene expression [22][23][24][25]. It has further been shown to be modified by neuronal activity [26]. Alterations of DNAm in epileptogenesis encompass upregulation of DNA-methyl-transferases, enzymes methylating the DNA base cytosine, in human TLE patients [27], genome wide changes in DNAm during epileptogenesis [28,29], and a later onset of spontaneous seizures in murine epilepsy models under treatment with a DNA-methyl-transferase inhibitor [30].
With the dawn of new site-and cell specific epigenetic modulators like modified CRISPR, zinc finger proteins and transcription-activator-like-effectors [31][32][33], the identification of potential genomic sites for antiepileptogenic intervention is of utmost importance.
DNAm in neurons and glial cells is mostly cell specific [34,35]. Sorting of brain tissue into specific cell types prior to downstream analysis has been applied in previous studies of gene expression [36][37][38] and DNA methylation [39][40][41]. This approach provides information about the cellular origin of observed effects on the epigenome and transcriptome level and an elevated detection sensitivity of more subtle changes in DNAm and GE.
Hypotheses for this study are i) that epileptogenesis affects DNAm and GE in a cell specific manner and ii) that differential methylation (DM) in neurons and glial cells correlates with differential gene expression (DGE). To our knowledge, this is the first study to investigate DNAm and GE changes as well as their possible association in neurons and glia separately, and the first one of its kind conducted in the widely used intracortical mouse model of mTLE [15].

Tissue collection and pooling
At 24 hours after status epilepticus, cervical dislocation was performed under anesthesia and right hippocampi were extracted. After extraction, each hemisphere was placed in a 2 mL polypropylene tube, instantly shock frozen in liquid nitrogen, and stored at -80˚C. Right hippocampi were pooled in 2 mL tubes from 4 (KA n = 4, S n = 4) or 2 (KA n = 4, S n = 4) mice prior to further processing. The number of mice per group (KA, SH) amounted to 8 mice per group, the number of biological samples to 3 per group (KA, SH). Tissue was kept on dry ice during pooling. See also S1 Fig.

Fluorescent Activated Nuclear Sorting (FANS)
A modified version of a nuclear sorting protocol by Jiang et al. [42] was used to sort into NeuN+ (refered to as neurons) and NeuN-(refered to as glia) nuclei. Immediately after pooling, hippocampi were placed on ice, and 1 mL homogenization buffer was added to each pool. Tissue was homogenized using a GentleMACS dissociator (Miltenyi), and homogenate filtered through a 70 μm filter. A NeuN-negative control sample from adult mouse liver was processed in parallel with the hippocampal samples. Debris was removed by density gradient centrifugation, using Debris Removal Solution (Miltenyi), and nuclear pellets were resuspended in 100 μL incubation buffer per one million nuclei. Anti-NeuN Alexa Fluor488 (Merck Millipore) was added to each sample to a final concentration of 0.1 μg/mL, and samples incubated for 1 h on ice in a light protected environment. Sorting of nuclei was performed on a FACSAria (BD Biosciences). Propidium iodide (PI) was added prior to sorting, and the following strategy was used for gating (S2 Fig): 1) A nuclear gate was defined by PI-positive events. 2) Aggregated nuclei were excluded in a dot plot using the pulse width of side scatter (SSC-w) versus the pulse area of forward scatter (FSC-a). 3) NeuN-negative gate was drawn based on signal from anti-NeuN stained liver sample. NeuN-positive and NeuN-negative hippocampi nuclei were sorted into tubes, and nuclei pelleted by centrifugation. Pellets were resuspended in lysis buffer for downstream DNA and RNA isolation. A full description of the FANS procedure is given in S1 Supporting Information.

Isolation of DNA and total RNA from sorted nuclei
DNA was extracted from sorted nuclei with MasterPure Complete DNA and RNA Purification Kit (Epicentre), DNA purity was assessed on NanoDrop, and DNA concentration measured on Qubit (DNA HS assay). Further details are available in S1 Supporting Information. Lysates were thawed on ice and total RNA extracted with mirVana miRNA Isolation Kit (Ambion). Up-concentration was performed using RNA Clean & Concentrator-5 kit (Zymo Research). RNA concentration and integrity were assessed on Bioanalyzer with the RNA Pico Kit (Agilent Technologies). Further details are found in S1 Supporting Information.

RRBS
A modified version of the gel-free protocol by Boyle et al. [43] was used for RRBS library preparation. The main changes to the protocol were inclusion of a two-sided size selection prior to bisulfite conversion, and sample pooling performed after completion of single libraries. A full description of the RRBS library prep and sequencing is given in S1 Supporting Information. Libraries were subjected to either 75 bp single read sequencing on NextSeq500 (Illumina), or 150 bp single read sequencing on HiSeq2500 (Illumina). For sequencing on NextSeq500, a pool of 14 libraries were run twice, with 50% PhiX spike-in at each run. On HiSeq2500, pools of 15 libraries were sequenced over two lanes, using 10% PhiX spike-in.

High throughput mRNAseq
SMART-Seqv4 Ultra Low InputRNA Kit for Sequencing (Takara Bio) was used to amplify mRNA from total RNA, and the resulting cDNA was used as input in library preparation with ThruPlex DNAseq Kit (Rubicon Genomics). See S1 Supporting Information for details regarding cDNA synthesis and library preparation. Libraries were sequenced on NextSeq500 (75 bp single read), or HiSeq3000 (150 bp paired end). On NextSeq500, 12 libraries were pooled for one sequencing-run. The remaining 27 libraries were sequenced in one pool over three lanes on HiSeq3000.

RRBS-and mRNAseq-post processing.
Post-processing included trimming of reads using Trim Galore! v0.4.3 with parameters "-rrbs-illumina" and quality control with FastQC. Alignment to the reference mouse genome mm10 was performed with Bismark v0.20, powered by bowtie2. Quality metrics were collected from the resulting BAM files using the Picard tool CollectRrbsMetrics v2. 18.15. Alignment of the mRNAseq reads was accomplished with the Subread package through its R interface Rsubread, after trimming with Trim Galore! v0.4.3. Quality control of the resulting BAM files was undertaken with CollectRnaSeqMetrics from Picard v2.18.15. Uniquely mapped reads were assigned to genes and counted by the featureCounts function of Rsubread, using default parameters. As reference for the gene assignment we used release M16 of the comprehensive gene annotation of mm10 available from GENCODE. Only RNA aligning to mRNA regions was used for further analysis.
Annotation. The mouse genome build mm10 was used as reference in all analyses. Only autosomal data were analyzed. Coordinates of genes, exons and introns were taken from GENCODE's comprehensive annotation (www.gencodegenes.org/mouse/release_M16.html). The R package annotatr [44] was used to bioinformatically link CpG sites to the gene annotations.
Using predefined genomic features within annotatr [44], the promoter region for any gene was defined as the segment from -1kb upstream to the transcription start site and the upstream region from -5 kb to -1 kb, where negative numbers indicate positions upstream of transcription start site.
we also performed an alternative estimate of the conversion rates directly from the untrimmed fastq files, by checking the methylation status specifically on the (unmethylated) cytosines added in the end-repair step of the RRBS preparation (private bash script). See S1 Supporting Information for further information of bisulfite conversion rate.
Multidimensional scaling. In order to validate our cell sorting procedures, and look for outliers among the samples, multidimensional scaling (MDS) plots were produced for the mRNAseq and RRBS data sets. The MDS computations were done by the plotMDS function of edgeR, selecting the top 100 most variable loci. The actual plots were created with ggplot2 [47].
Expression of neuronal and glial genes in NeuN+ and NeuN-fraction. Normalized counts for expression of neuronal (RBFOX3), glial (ALDHL1L1, CX3CR1, MBP), pericytal (PDGFRB) and endothelial (PECAM1) genes were used to visualize enrichment of neurons in the NeuN+ fraction and glia in the NeuN-fraction.

Selection of relevant GO and KEGG terms
Relevant (Tables 1 and 2) and epilepsy-relevant (Fig 2) GO and KEGG terms in neurons and glia were selected manually based on reviews on the subject [9] and personal knowledge. The list of GO and KEGG terms derived from our DGE analysis (for a full list see S1 Table) was manually filtered for specific terms and relevant up-/downregulated genes within these terms in neurons and glia selected for presentation.

Quality control
Bisulfite conversion rates were above 98% (S1 Supplementary Information and S3 Fig)

Differential methylation
A statistical analysis of Differentially Methylated CpGs (DM CpGs) compared right hippocampi of KA to SH mice at 24 hours post injection. After filtering, 928 430 CpG sites remained and were used in subsequent analyses. On average, across all CpG sites and all samples, each CpG was covered by 29.8 reads. In individual samples the mean read depth varied from 20.0 to 35.2 (median = 30.6, interquartile range = [27.4-32.7]).
Differentially methylated sites. The analysis of significantly altered DM CpGs revealed 1060 hyper-and 899 hypomethylated (ratio 1.2:1) CpG sites in neurons and 464 hyper-and 274-hypomethylated (ratio: 1.7:1) CpG sites in glia (Fig 1). Most of the DM CpGs localized to either gene bodies or intergenic regions (for full list see supplementary 2) and were distributed evenly across chromosomes (sex chromosomes excluded), apart from a possible higher ratio of hyper-/hypomethylated CpGs at chromosome 13 to 15 in glia. The ratio of hypermethylated to hypomethylated CpG sites was highest at upstream (1.3:1) and intergenic (1.2:1) regions for neurons and upstream (2.7:1) and promoter (2.0:1) for glia. For detailed information including GO and KEGG annotation of DM CpGs see S1 Table. Differentially methylated CpG sites common to both neurons and glia. Neurons and glia shared four commonly hypermethylated (fraction: 0.3%) and zero commonly hypomethylated DM CpGs. One DM CpG was hypermethylated in neurons and hypomethylated in glia (fraction: 0.0%) and one DM CpGs hypomethylated in neurons and hypermethylated in glia (fraction: 0.0%). Three of the four commonly hypermethylated DM CpGs localized to gene bodies and one to an upstream region of the associated gene. The other two CpGs were associated with intergenic regions. Differentially methylated regions. In order to obtain information about genomic features with significantly altered differential methylation, an aggregated analysis was conducted for each genomic feature (upstream, promoter, UTR5, exon, intron, gene body, UTR3) separately. Most DMR were found in gene bodies (e.g. exonic and intronic areas) and promoters in both neurons and glia. The greatest hyper-/hypomethylated regions ratio was found at 5'UTRs (ratio 4.5:1), gene bodies (ratio 3.7:1) and promoters (ratio 3.1:1) in neurons and at promoters (ratio 8.2:1), 3'UTRs (ratio 7.0:1) and exons (ratio 5.3:1) for glia. For a full list of significantly altered DMR and their annotated GO and KEGG terms see S1 Table. Differential gene expression DGE analysis compared right hippocampi of KA to SH at 24 hours post injection. After filtering, 23369 genes were used for downstream analysis. After processing, alignment and filtering, the mRNASeq samples yielded on average 7.5 giga bases (Gb) data aligned to the mouse genome (median = 9.0 Gb; range = [2.4 Gb-13.0 Gb]). The fraction aligning specifically to mRNA regions varied from 16% to 41.2% (median = 28.5), resulting in an average of 2.2 Gb per sample informative for DGE analysis (median = 2.2 Gb; range = [0.4 Gb-4.3 Gb]). In neurons, 135 genes were up-and 15 downregulated (Table 1), while in glia 147 genes were up-and 85 downregulated (Table 2). A relevant selection of broader GO / KEGG terms is presented in Tables 3 (neurons) and 4 (glia). For neuronal and glial contribution to epileptogenesis in terms of number of differentially expressed genes as part of epilepsy-relevant GO / KEGG terms, see Fig 2. For a detailed list of up-and downregulated genes and associated GO / KEGG terms see S1 Table. Genes up-or downregulated in both neurons and glia. A total number of 45 genes were upregulated in both neurons and glia. GO terms of these included "positive regulation of transcription" and "cytokine mediated signaling" whilst KEGG terms contained "ECM receptor interaction", "VEGF-signaling" and "TNF-signaling". Three genes were downregulated in both neurons and glia.

Association between differential methylation and differential gene expression
An association between DM and DGE was calculated by alignment of significantly altered DMR and differentially expressed genes (Combined FDR 0.25, for details see S1     (upstream, promoter, UTR5, exon, intron, UTR3) was found (S1 Table and S7-S20 Figs). However, significant DMR coincided with DGE at 41 loci for neurons and 12 loci for glia (Tables 5 and 6).

Discussion
Neurons and glial cells execute specific complementary tasks in normal brain functioning as well as in the pathological processes precipitating neurological diseases [9, 48,49]. This diversity is represented on the transcriptome-[36, 50] and epigenome level [34,35]. In this study we investigated neuronal and glial alterations of DNAm and gene expression and their possible association in a mouse model of epilepsy. In order to explore the epigenetic and transcriptomic signature of both cell types in early epileptogenesis, we separated neurons and glia by FANS [51]. This approach was recently applied in epigenetic studies [39-41]. We identified specific neuronal and glial DNAm and DGE changes at particular genomic loci, potentially including important upstream mechanisms worth further investigation.

Differential methylation
With an overlap of only 0.22% of DM CpGs between neurons and glia, differential methylation occurs primarily in a cell specific manner during early epileptogenesis. Accordingly, the attributed GO terms to DM CpGs and DMR reveal mostly neuronal and glia specific terms (S1 Table). Apart from a near even ratio at neuronal promoters, we observe an overweight of hypermethylated to hypomethylated CpG sites at 24 hours post injection. In early stages of epileptogenesis, previous studies suggest no general alteration of DNAm [30] or a tendency towards hypomethylation [52]. In chronic phases of epilepsy, both hyper-and hypomethylation have been reported [29,53]. Differences in pro-convulsant agents, mode and region at which they are applicated, stage of epileptogenesis, anatomical regions investigated, as well as the number of CpG sites covered, may account for varying results. A comparison of genes associated with significant DM CpGs and DMR from our study with results from hippocampal tissue from chronic TLE patients [54] unvails marginal but interesting overlaps. This may be due to different methylation analysis methods and stage specific character of hippocampal DNAm during epileptogenesis. A comparison of our results to those from peripheral blood of TLE patients [55], reveals 15 overlapping DM CpGs and 11 overlapping DMRs.

Differential gene expression
Most differentially expressed genes detected in this study, were specific to either neurons or glia (S1 Table). Only 45 genes were commonly up-and three downregulated in neurons and glia. At this early stage of epileptogenesis, glial cells apparently contribute a higher number of altered gene transcripts than neurons within epilepsy-related GO and KEGG terms (Fig 2). Many significantly differentially expressed genes found in this epilepsy model overlap with data from other experimental and human TLE studies [29,[56][57][58]. Comparing our differentially expressed genes in neurons and glia with differential gene expression from various epilepsy models and stages of TLE [58], glia contributes with more upregulated gene transcripts (13 glia and 11 in neurons), supporting the notion that glia contributes essentially to epileptogenesis [59][60][61][62][63][64][65][66]. Neurons contribute with slightly higher number of altered gene transcripts when comparing to a amygdala stimulation model of epilepsy while glia contributes with a higher number of differentially expressed genes when comparing to a traumatic brain injury model of epilepsy [29]. Within up-and downregulated gene transcripts from a study on refractory human TLE, neurons and glia exhibit an almost even number of gene transcripts [57]. For an overview of essential up-and downregulated pathways see Tables 1 and 2. A summary of neuronal and glial contribution (number of genes within a GO / KEGG term) to epileptogenesis is presented in Fig 2.

Altered epilepsy-relevant pathways based on differentially expressed genes
Upregulation of growth arrest and DNA-damage-inducible beta/gamma (GADD45B/G). One of the genes upregulated in both neurons and glia is GADD45G, a member of the environmental stress inducible GADD45-like genes that mediate activation of various pathways, including c-Jun N-terminal protein kinase family of mitogen-activated protein kinases [67]. It has previously been shown to be elevated after KA induced status epilepticus [68] and to possess DNA demethylation qualities [69], thus potentially linking epileptic activity to changes in DNA methylation. We also find elevated mRNA levels of GADD45B, which is another member of the GADD45-like genes. GADD45B has recently been shown to promote neuronal activity induced neurogenesis via demethylation of the BDNF and fibroblast growth factor promoter, linking neuronal activity to DNA methylation alterations [69].

Upregulation of sphingosine-kinase 1 (SPHK1) and sphingosine 1 receptor 3 (S1R3)
Another interesting finding is the upregulation of SPHK1 mRNA in neurons and glia. SPHK1 phosphorylates sphingosine to sphingosine-1-phosphate (S1P) [70]. S1P in turn is involved in neural development, signaling, autophagy and neuroinflammation as well as a plethora of pathological central nervous conditions [71] and has been shown to modulate histone deacetylase activity [72]. A recent study revealed antiepileptogenic effects of fingolimod [71], a SP1-receptor modulator and FDA approved drug for the treatment of multiple sclerosis [73], possibly via attenuation of astro- [74] and microglial [75] reactions. Further, we find elevated expression of S1R3, a S1P receptor, in glia. This receptor has been shown to be elevated in hippocampi of kainic acid and pilocarpine epilepsy models as well as in humans with TLE, and is mainly expressed in astrocytes [76].

Upregulated mitogen-activated protein kinase (MAPK) pathways in neurons and glia
MAPK is a type of protein kinase specific to the amino acids serine and threonine. MAPKs are involved in directing cellular responses to a variety of different stimuli, such as proinflammatory cytokines, mitogens, osmotic stress and heat shock. They regulate cell functions including proliferation, gene expression and differentiation, mitosis, cell survival and apoptosis [77]. We find several genes within MAPK pathways (GO / KEGG) to be expressed more in the KA than SH group with glia contributing a higher number of differentially expressed genes within the pathways than neurons (S1 Table). In the context of epilepsy, MAPK are thought to play a role in Cx43 phosphorylation, involving TNF-α, interleukin (IL)-1b [78] and VEGF [79]. We find elevated expression levels of genes within KEGG pathways for both TNF-α (mostly in glia) and VEGF (slightly more in neurons, see S1 Table). Phosphorylation of Cx43 in turn has been associated with its elevated internalization and degradation [80], possible contributing to astrocyte uncoupling in both mTLE mice and humans [15]. Astrocytic calcium signaling pathways altered in glia. Neuronal activity induced elevations in astrocytic intracellular calcium levels may in turn facilitate astrocytic release of neuroactive substances including glutamate, aggravating epileptic activity [81]. In acute stages of epileptogenesis, calcium transients in astrocytes are increased, possibly contributing to elevated extracellular potassium levels via Calcium-dependent protease induced cleavage of the dystrophin associated protein complex [82,83]. Elevated extracellular potassium levels in turn may lead to increased excitability of neurons and thereby generate epileptiform activity [84]. We also find elevated gene expression of inositol 1,4,5-trisphosphate 3-kinase A (ITPKA), a protein kinase inactivating inositol triphosphate dependent calcium release from the astrocyte endoplasmic reticulum [85,86], in glia. Further, we find 1-Phosphatidylinositol-4,5-bisphosphate phosphodiesterase epsilon-1 (PLCE1), a member of the phosphatidylinositol-specific phospholipase C family that via G-protein coupled receptors are involved in Inositol-triphosphate and diacylglycerol generation and as such mediate intracellular Calcium elevation [87], elevated in glia. Thirdly, we find CACNG5, a calcium permissive AMPA receptor subunit [88], elevated in glia. Possibly all three genes mediate pro-epileptogenic effects via astrocytic calcium signaling. For a complete summary of differentially expressed genes see S1 Table.

Relationship between differential methylation and differential gene expression
We did not find a general correlation of DM and DGE at specific genomic regions (S1 Table, S7-S20 Figs). However, DM coincided with DGE at 41 genes for neurons and 10 genes for glia. Our results are in line with previous studies that did not find a general correlation between DNA methylation and gene expression in epilepsy, but rather a number of singular genes where significant DM with DGE coincided [30,89,90]. Other studies did report a certain degree of general association between DM at specific genomic regions and DGE [29]. For a full list of coinciding alterations in DM and DGE see Table 5. The identified coinciding alterations of DM and DGE in this study are relatively few and their role in epileptogenesis remains uncertain. Interestingly, they point to genes and pathways previously implicated in epilepsy, TLE and epileptogenesis. In the following we present a selection of these in depth. Coinciding alterations of differential methylation and differential gene expression in neurons. Osteopontin (SPP1) promoter hypomethylation associated with elevated gene expression. SPP1 mediates diverse aspects of cellular functioning in the central nervous system, e.g. the recruitment and activation of microglia and astrocytes, the cumulative effect possibly being neuroprotective [91]. In multiple sclerosis, SPP1 mediates pro-inflammatory pathways contributing to the relapse remission phenotype via e.g. NF-κB [92]. In our study, SPP1 mRNA is significantly upregulated in both neurons and glia, and in neurons this elevated gene expression is associated with significant hypomethylation at the associated upstream region and promoter. These results confirm previous findings of elevated SPP1 mRNA levels in epilepsy [93]. We further found elevated mRNA levels of CD44, a Osteopontin receptor [94] involved in epileptogenesis [95], in glia.
Hypermethylation at a Galanin (GAL) exon associated with elevated gene expression. GAL encodes the neuropeptide Galanin which previously has been shown to possess seizure attenuating properties and discussed as a possible antiepileptogenic target [96]. Further, in a recent study, a de novo mutation in GAL has been unveiled as a possible cause for TLE [97].
In our study, GAL mRNA is significantly upregulated in both glia and neurons (slightly higher log fold change (logFC) and lower FDR in neurons) and we find a significant association of hypermethylation of an exonic region of GAL with elevated gene expression in neurons. Thus, the hypermethylation at the exonic region of GAL with possible consecutive elevated levels of GAL might represent a crucial endogenic seizure attenuating mechanism.

Hypomethylation at synaptic vesicle protein 2 c (SV2C) exon associated with upregulated gene expression
SV2C is, together with SV2A and SV2B, part of the family of synaptic vesicle proteins that are involved in Ca 2+ dependent synaptic vesicle exocytosis and neurotransmission [98]. The most prominent epilepsy-related member, SV2A, is the main target through which levetiracetam and brivaracetam exert their antiepileptic and possibly antiepileptogenic effects [99]. In hippocampi of patients with chronic TLE, SV2C was the only of three synaptic vesicle proteins found to be significantly elevated. It was associated with mossy fiber sprouting and glutamatergic synapses and was proposed as a potential antiepileptogenic target [100]. Recent findings suggest a role of SV2C in the disruption of dopamine signaling in Parkinson's Disease [101]. At 24 hours post injection, we find elevated levels of SV2C mRNA in glia and neurons. In neurons this elevated gene expression is associated with significant hypomethylation of its exonic regions. Hypomethylation of the SV2C exon may thus exert upstream pro-epileptogenic effects.

Coinciding alterations of differential methylation and differential gene expression in glia
Promoter hypermethylation at HDAC11 associated with reduced gene expression levels. In line with previous results [102], we find HDAC11 mRNA levels decreased after SE. This reduced gene expression coincides with a significantly increased methylation at its associated promoter in glia. Reduced levels of HDAC11 may cause an increased acetylation at H4 [103], previously shown to correlate with elevated levels of c-fos, c-jun and BDNF [104]. BDNF in turn has been associated with seizure-aggravating effects in acute phases of epileptogenesis [105] and higher levels of the microRNA miR-132, possibly via ERK and MAPK pathways [106]. Further, miR-132 has recently been associated with seizure induced neuronal apoptosis [107]. We find elevated expression levels of BDNF (only glia), miR-132 (both neurons and glia) and MAPK-(glia more transcripts than neurons) and ERK-(glia more than neurons) pathways in early epileptogenesis. The hypermethylation of the HDAC11 promoter with its possible downstream effects ultimately leading to elevated levels of BDNF and miR-132 might represent a possible antiepileptogenic target for site specific alteration of DNAm.
Hypermethylation at the intron and promoter at dopamine receptor D1 (DRD1) associated with elevated DRD1 mRNA levels. Dopamine exerts its seizure inducing effects via DRD1 mediated ERK1/2 pathways [108]. We find elevated levels of DRD1 mRNA in both neurons and glia. In glia, this augmentation in gene expression is associated with significant hypermethylation in the intronic region and hypermethylation at the promoter of DRD1. Hypermethylation at glial DRD1 (intronic region) may facilitate epileptogenesis.

Technical limitations
This study features several technical limitations worth mentioning. Firstly, we cannot rule out that the NeuN-fraction contains a minor number of non-glial cells (pericytes, endothelial cells) [109,110]. RRBS associated technical limitations involve loss of information associated with e.g. msp1 enzyme cleavage, library pooling/ fragment size selection, bisulfite conversion, sequencing (depth/coverage) [43]. The use of nuclear mRNA results in enrichment of mRNA coding for proteins with nuclear functions [111] and a potentially lower level of immediate early genes [112]compared to when using cytosolic mRNA [113].

Conclusion
In this study, we found DNAm and DGE in early epileptogenesis to occur primarily in a cellspecific manner. We identified several potential neuronal and glial upstream targets worth further investigation. Information on the cellular origin of epigenomic and transcriptomic effects increases our understanding of involved pathological processes and provides a basis for possible future cell specific therapeutic approaches.
Supporting information S1 Table. DM and DGE in neurons and glia. Table of  (TIF) S11 Fig. DM and DGE (neurons, UTR5). Visualization of DM and DGE (FDR 0.25) for neurons (UTR5). Genes associated with significantly altered DMR and DGE are indicated in the figure.
(TIF) S13 Fig. DM and DGE (neurons, exon). Visualization of DM and DGE (FDR 0.25) for neurons (exon). Genes associated with significantly altered DMR and DGE are indicated in the figure.
(TIF) S14 Fig. DM and DGE (glia, exon). Visualization of DM and DGE (FDR 0.25) for glia (exon). Genes associated with significantly altered DMR and DGE are indicated in the figure.
(TIF) S15 Fig. DM and DGE (neurons, intron). Visualization of DM and DGE (FDR 0.25) for neurons (intron). Genes associated with significantly altered DMR and DGE are indicated in the figure.
(TIF) S16 Fig. DM and DGE (glia, intron). Visualization of DM and DGE (FDR 0.25) for glia (intron). Genes associated with significantly altered DMR and DGE are indicated in the figure.