New differentially expressed genes and differential DNA methylation underlying refractory epilepsy

Epigenetics underlying refractory epilepsy is poorly understood, especially in patients without distinctive genetic alterations. DNA methylation may affect gene expression in epilepsy without affecting DNA sequences. Herein, we analyzed genome-wide DNA methylation and gene expression in brain tissues of 10 patients with refractory epilepsy using methylated DNA immunoprecipitation linked with sequencing and mRNA Sequencing. Diverse distribution of differentially methylated genes was found in X chromosome, while differentially methylated genes appeared rarely in Y chromosome. 62 differentially expressed genes, such as MMP19, AZGP1, DES, and LGR6 were correlated with refractory epilepsy for the first time. Although general trends of differentially enriched gene ontology terms and Kyoto Encyclopedia of Genes and Genome pathways in this study are consistent with previous researches, differences also exist in many specific gene ontology terms and Kyoto Encyclopedia of Genes and Genome pathways. These findings provide a new genome-wide profiling of DNA methylation and gene expression in brain tissues of patients with refractory epilepsy, which may provide a basis for further study on the etiology and mechanisms of refractory epilepsy.


INTRODUCTION
65 millions of people were affected by epilepsy in the world according to International League Against Epilepsy (ILAE) [1], and approximately 36% of epilepsy patients were drug-resistant [2]. Varieties of genes encoding channels, receptors, transporters, synaptic transmission, etc. have been associated with different types of epilepsy, among which, some were associated with refractory epilepsy [3,4,5,6]. Many environmental factors, such as economic situation, diet, trauma, stroke, etc. were also associated with epilepsy or seizure [7,8,9,10].
Epigenetic modifications, including DNA methylation, histone modification and aberrant microRNA expression, can affect genomic reprogramming, tissue-specific gene expression and global gene silencing without affecting DNA sequence itself [11,12]. The most common form of DNA methylation occurs at the 5'carbon of cytosine in CpG dinucleotides which often locates in CpG islands within the promoters [13]. More recently, DNA methylation is raised as one of the main epigenetic mechanisms in epilepsy [14]. Previous genome-wide DNA methylation profiling in epileptic animal models presented altered DNA methylation in promoters of genes, and identified many genes that were associated with epilepsy [15].
DNA methylation in promoter may decrease gene expression, for example, increased methylation of reelin promoter resulted in the decrease of reelin expression in epilepsy model [16]. Katja Kobow and colleagues found a genome-wide distinctive DNA methylation pattern in rat Research Paper www.impactjournals.com/oncotarget models of pilocarpine-induced epilepsy in addition to an inverse relationship between gene expression and DNA methylation in promoter, exon and intron [17]. Meanwhile, ketogenic diet could attenuate seizure progression via ameliorating DNA methylation [17]. Selective changes in genome-wide DNA methylation and increased DNAmethyltransferase were also discovered in patients with temporal lobe epilepsy (TLE) [18,19]. However, the profiling of genome-wide DNA methylation and gene expression in patients with refractory epilepsy remains unclear.
In this study, methylated DNA immunoprecipitation linked with sequencing (MeDIP-Seq) and mRNA sequencing (mRNA-Seq) were used to analyze the pattern of genome-wide DNA methylation and gene expression, as well as the relationship between DNA methylation and gene expression. Our findings identified a new distribution pattern of DNA methylation and gene expression in refractory epilepsy patients. Most of the differentially methylated genes (DMG) were methylated in gene element of coding sequences (CDS) and introns. More importantly, some new refractory epilepsy-related genes that have not been documented previously were found in this study.

Demographic and clinical characteristics of subjects
The mean age (mean±SD) of the 10 epileptic samples (5 male/5 female) was 17.10±5.84, and the age of epilepsy onset was 6.49±6. 16. The mean age of the 10 controls (7 male/3 female) was 39.00±17. 40. The detailed data of demographic and clinical characteristics of epileptic samples were presented in Table 1 &  Supplementary Table S1.

No significant difference in distribution of DNA methylation reads
In each sample, 81632654 methylation reads (49 bp) were sequenced. In epileptic samples, an average 71.20% of the reads were uniquely mapped to the reference genome, and in controls, 70.77% of the reads were uniquely mapped. There was no significant difference of uniquely mapped reads between the two groups by T test (Supplementary Table S2).
Moreover, no significant difference was identified between epileptic samples and controls in 1) genome coverage distribution across sequencing depth, 2) distribution of CpG, CHG, and CHH sites across sequencing depth, 3) reads distribution in genome regions with different CpG density, 4) distribution of reads in different gene elements and repetitive elements, 5) distribution of reads around CpG island and gene body.  Table S2). No significant difference between the two groups was found. In addition, no significant difference in the number of peaks with different length, the distribution of peaks with different CpG density and distribution of peaks in gene elements (including peak number and peak coverage) was observed between epileptic samples and control.

Analysis of differentially methylated regions (DMR) and DMGs
The median of DMRs identified by pairwise comparison was 7604.50, covering a median of 8580791.00 bp. The distribution of DMGs in paired samples was mapped to human genome using Circos [20] (Figure 1). DMGs appeared on all of chromosomes extensively except for Y chromosome in refractory epilepsy patients.
The significant enriched gene ontology (GO) terms of DMGs were presented in Table 2. Most of the DMGs were differentially methylated in gene element of CDS and introns. Significant enrichment of DMGs was observed in GO terms of binding of various molecules, such as ATP binding, ion binding, cation binding, and nucleoside binding. In addition, DMG enrichment was also identified in the GO terms involved in receptor activity, transporter activity, kinase activity, transducer activity and channel activity. Those DMGs participate mainly in the biological processes of adhesion and ion transport.

Gene expression profiling
In each epileptic sample, an average of 78044550.80±8806016.67 reads covering 7024009572.00±792541500.40 bp were sequenced, among which an average 83.71% of reads were uniquely mapped to reference genome, and 66.95% of reads were uniquely mapped to reference genes. In controls, an average of 79185195.80±9170582.72 reads covering 7126667622.00±825352445.10 bp were sequenced, among which an average 83.80% of reads were uniquely mapped to reference genome, and 65.99% of reads were uniquely mapped to reference genes. No significant difference was found between the two groups by T test (p<0.05 was considered statistically significant). (Supplementary Table S3). www.impactjournals.com/oncotarget   glutamate receptor activity; GTPase regulator activity; guanyl-nucleotide exchange factor activity; ion binding; ion channel activity; ion transmembrane transporter activity; ionotropic glutamate receptor activity; kinase activity; ligand-gated channel activity; ligand-gated ion channel activity; metal ion binding; metal ion transmembrane transporter activity; molecular transducer activity; nucleoside binding; nucleoside-triphosphatase regulator activity; passive transmembrane transporter activity; phosphoric diester hydrolase activity; phosphotransferase activity, alcohol group as acceptor; protein kinase activity; protein tyrosine kinase activity; protein tyrosine phosphatase activity; purine nucleoside binding; Rho GTPase activator activity; signal transducer activity; transferase activity, transferring phosphorus-containing groups; transmembrane receptor protein kinase activity; transmembrane receptor protein tyrosine kinase activity P cell adhesion; biological adhesion; homophilic cell adhesion; cell-cell adhesion; ion transport; metal ion transport Promoter C synaptosome F transmembrane receptor activity; receptor activity; G-protein coupled receptor activity; signal transducer activity; molecular transducer activity; olfactory; voltage-gated ion channel activity; voltage-gated channel activity P G-protein coupled receptor protein signaling pathway; cell surface receptor linked signal transduction; sensory perception of chemical stimulus; inorganic; anion transport www.impactjournals.com/oncotarget A total of 21353 genes were sequenced, among which 17665 genes were expressed by all samples. The sequencing coverage of 65.01% and 65.45% genes was over 90% in epileptic samples and controls, respectively, indicating a good sequencing quality.
Pairwise comparison presented a total of 395 GO terms significantly enriched in epileptic samples compared to controls (Component 68, Function 70, Process 257). 77 GO terms were significantly enriched in ≥5 epileptic samples compared to controls, including 12 component terms, 23 function terms and 42 process terms ( Table  5). The GO enrichment analysis of gene expression revealed that the biological functions of DEGs are mainly correlated with channel activity, transporter activity, and receptor activity. Furthermore, the biological processes mainly targeted by DEGs were immune system, biological regulation, response to stimulus, signaling, development, and behavior.
In Kyoto Encyclopedia of Genes and Genome (KEGG) pathway enrichment analysis, 75 pathways were significantly enriched in epileptic samples compared to controls, while 27 pathways were significantly enriched in ≥5 epileptic samples compared to controls (Table 5). DEGs mainly participate in calcium signaling pathway, neuroactive ligand-receptor interaction, and pathways involved in inflammation, immune response, and autoimmune diseases.

DNA methylation and gene expression
The distribution of hyper-, hypo-and unmethylated gene expression levels in different gene elements were presented in Figure 3. The trend of gene expression of the three groups were similar, and in all the 5 elements, the percentage of hyper-methylated genes with log 2 (RPKM ratio) (RPKM ratio = RPKM of epileptic sample/RPKM of control) approximately ranging from 1 to 3 were higher than that of hypo-methylated and unmethylated genes. Generally, no significant relationship in modulation was found between DNA methylation and gene expression.

DISCUSSION
It is the first genome-wide report on DNA methylation and gene expression in refractory epilepsy patients. 62 differentially expressed genes such as MMP19, AZGP1, DES, and LGR6 were discovered to be correlated with refractory epilepsy, and diverse distribution of differentially methylated genes was found in X chromosome instead of in Y chromosome.
Although many genes were observed differentially methylated or expressed, the general distribution of DNA methylation reads, DNA methylation peaks and mRNA sequencing reads were similar between refractory epileptic samples and controls, which indicate no significant   In pairwise comparison of gene expression analysis, we identified distinct gene expression signatures. 34 DEGs are correlated with epilepsy or seizure and 14 of these genes are associated with refractory epilepsy, such as AQP1 [21], CCR5 [22], EMP1 [23], CXCL8 [24], ITGA2 [25], and CCL2 [26]. For the first time, 62 DEGs differentially expressed in ≥8 pairs of samples were found related to epilepsy/seizure in our study. These newlyidentified refractory epilepsy-related genes may possibly reveal new mechanisms of refractory epilepsy. In all the DEGs, only MMP19 and AZGP1 were differentially expressed in all the 10 pairs. Compared to the controls, 9 of the 10 epileptic samples showed increased expression   of MMP19, while 1 epileptic sample showed decreased expression. MMP19 and other matrix metalloproteinases can cleave and remodel the extracellular matrix, including tenascin and laminin, and thus influence synapse formation and remodeling, N-methyl-D-aspartate receptor activity, learning and memory, and hippocampal long-term potentiation [27]. Inhibition of MMP19 and other matrix metalloproteinases may prevent development of epilepsy at the early stage of epileptogenesis [28]. Meanwhile, we found the expression of AZGP1 were decreased in 8 of the 10 epileptic samples, and increased in 2 of the 10 epileptic samples when compared to controls. AZGP1 encodes Zinc-a2-glycoprotein, which is an adipokine participates in lipid mobilization, lipolytic effect, and immune response [29,30]. Moreover, both DES in 9 pairs and LGR6 in 8 pairs showed consistently increased expression in epileptic samples. To the best of our knowledge, it is the first time that MMP19, AZGP1, DES, and LGR6 were reported to be correlated with refractory epilepsy. Significant enrichment of DMGs in GO terms of binding, transport, and enzymatic activity were found, which is consistent with previous studies [18]. Interestingly, most of the DMGs were differentially methylated in CDS and intron, while previous research showed differential methylation in all the gene elements in rat models of chronic epilepsy induced by pilocarpine [17]. These findings indicate DNA methylation in CDS and intron may play critical roles in refractory epilepsy besides promoter methylation which has been a very popular target in research on epilepsy [18,31]. The GO enrichment analysis of gene expression revealed a trend similar to a previous report [17] that the DEGs are mainly correlated with biological functions such as protein binding, receptor binding, channel activity, transporter activity, and receptor activity, as well as being involved in biological processes such as immune system, biological regulation, response to stimulus, signaling, development, and behavior.
The change of DNA methylation in this study is not exactly corresponded with alteration of gene expression. Kobow and his colleagues found that DNA methylation in promoter, exon and intron were inversely correlated with gene expression in rat models of chronic epilepsy induced by pilocarpine [17], but our study found that most of the hyper-and hypo-methylated genes were not differentially expressed in epilepsy patients, and the percentage of hyper-methylated genes with log2(RPKM ratio) ranging from 1 to 3 were higher than that of hypomethylated and unmethylated genes, which indicate a complicated modulation between DNA methylation and gene expression in refractory epilepsy in human beings.
In KEGG analysis of gene expression, DEGs significantly enriched in calcium signaling pathway, neuroactive ligand-receptor interaction, and pathways participating in inflammation, immune response, autoimmune diseases. Calcium signaling pathway has been increasingly recognized as a vital factor in epileptogenesis and the excess synchronization, and hyperexcitability of neurons for seizures can be linked to various calcium signaling pathways [32]. The aberrantly neuroactive ligand-receptor interaction can enhance the susceptibility to epileptic seizures [33,34]. It may explain partially that the drugs regulating the function of calcium signaling pathway and neuroactive ligand-receptor interaction are able to alleviate the seizure frequency [35]. The roles of immune response and inflammation in epilepsy have been recognized in previous studies [36]. Autoimmune epilepsy frequently present drug-resistance which can be controlled by immunosuppressive and immunomodulatory therapies [36,37]. Consistent with Lukic and his colleagues' study, we also identified that prion disease was significantly targeted by DEGs which indicate both refractory epilepsy and prion diseases may share some common pathway [38].

Study approval
The research protocol was approved by the Ethics Committees of the Second Affiliated Hospital of Chongqing Medical University. All procedures performed in studies involving human participants were in accordance with the ethical standards of the institutional and/or national research committee and with the 1964 Helsinki declaration. Written informed consent was obtained from all individual participants included in the study or their proxies.

Patients and tissues preparation
Resected brain tissues were retrospectively but consecutively collected from 10 patients with refractory epilepsy and 10 patients with post-traumatic intracranial hypertension who underwent surgical treatment since 2008 to 2014. Patients with refractory epilepsy were diagnosed following the definition of ILAE [39]. Briefly, all patients were resistant to maximum doses of at least three antiepileptic drugs (AED), and evaluated by detailed history, neurological examination, neuropsychological test and neuroimaging data. For presurgical evaluation and epileptogenic zones identification, a combined assessment of ictal simiology, brain magnetic resonance imaging, video-electroencephalography, sphenoidal electrode monitoring and intracerebral electroencephalography and intraoperative electrocorticography were applied. After evaluation, standard en bloc resection was performed. No refractory epilepsy patient received adjustment of AEDs during the 2 months before surgery. Brain tissues as control from the 10 post-trauma intracranial hypertension patients were identified normal by neuropathologist. These patients had no history of epilepsy or exposure to AEDs.
All the resected brain tissues were immediately frozen in liquid nitrogen and then stored at -80°C.
The 10 epileptic samples and 10 controls were paired, the difference in genome-wide DNA methylation and gene expression between the paired samples were analyzed using MeDIP-seq and mRNA-seq.

DNA methylation profiling
Genomic DNA was extracted using QIAamp DNA Micro Kit (Qiagen, Hilden, Germany) according to the manufacturer's instruction. Extracted DNA was fragmented to a size of 100-500 bp by sonication (Bioruptor NGS, Digenode, Liege, Belgium), and subjected to DNA-end repair, 3'-dA overhang and ligation of sequencing adaptors according to manufacturer's instruction (Paired-End DNA Sample Prep kit, Illumina, San Diego, USA) and denatured to single-stranded. Then the methylated DNA was immunoprecipitated by 5mc antibody (Magnetic Methylated DNA Immunoprecipitation kit, Diagenod, Liege, Belgium). After Real-time Quantitative polymerase chain reaction (PCR) (TaqMan Probe, Applied Biosystems, Thermo Fisher Scientific, Waltham, USA) validation and quality control of sample library (Agilent 2100 BioAnalyzer, Agilent, Santa Clara, USA), electrophoretically selected DNA fragments sizing from 200-300 bp were subjected to high-throughput sequencing (Illumina HiSeqTM 2000, Illumina, San Diego, USA). Sequencing strategy was Single-end 50 bp, and reads size was 49 bp.
To describe the distribution of MeDIP-Seq data on genome, the following items were calculated: 1) Genome coverage distribution across sequencing depth; 2) Distribution of CpG, CHG and CHH sites varies with sequencing depth; 3) Reads distribution in genome regions with different CpG density; 4) Distribution of reads in different gene elements, including CpG islands, promoters, 5'-Untranslated regions (UTR), CDS, introns, 3'-UTR, repeat regions and each class of repetitive elements (Repeat dataset is obtained from RepeatMasker (Transposons) and Tandem Repeats Finder (Tandom repeats), and is available at: http://hgdownload.cse. ucsc.edu/goldenPath/hg19/bigZips/chromOut.tar.gz); 5) Distribution of MeDIP-Seq reads around CpG island and gene body.
Whole genome scanning of enrichment region of methylation / Peak was based on a defined analysis model, MACS 1.4.0 (Website: http://liulab.dfci.harvard. edu/MACS/) with default parameters [41]. The following items were calculated: 1) Distribution of peaks with different length; 2) Distribution of peaks with different CpG density; 3) Number and coverage of peaks in gene elements (promoter, 5'-UTR, CDS, intron, and 3'-UTR).
DMGs based on peak of all the paired samples were analyzed. Briefly, peaks of the paired two samples were merged as candidate DMRs. For each candidate DMR, the number of reads of each sample was calculated and tested to get true DMRs. The downtrend DMRs indicate that the number of reads of the control sample was larger than the epileptic sample and the uptrend DMRs indicate the opposite. DMGs were defined as the genes overlapping DMRs. Genes with an element merely overlap uptrend DMRs were considered hyper-methylated with such an element, genes with an element merely overlap downtrend DMRs were considered hypo-methylated with such an element.
To clarify the biological functions of DMGs, GO enrichment analysis was performed. Briefly, DMGs were mapped to GO terms in reference database (http://www. geneontology.org) and gene numbers for every term were calculated and tested to identify the significantly enriched GO terms.

Gene expression profiling
RNA was extracted using Trizol method. After DNase I treatment, mRNA was isolated with magnetic beads with Oligo (dT) and fragmented. Then cDNA was synthesized using the mRNA fragments as templates. The synthesized cDNA fragments were purified and subjected to end reparation, single nucleotide adenine addition and adapter connection. cDNA fragments suitable for PCR amplification were selected with electrophoresis. Quality control of sample library was performed using Agilent 2100 Bioanaylzer (Agilent, Santa Clara, USA) and Applied Biosystems StepOnePlus Real-Time PCR System (Applied Biosystems, Thermo Fisher Scientific, Waltham, USA). The library was sequenced using Illumina HiSeqTM 2000 (Illumina, San Diego, USA).
After sequencing quality control, the mRNA-Seq data was mapped to reference genome and reference genes using SOAP software, version 2.21 (Website: http:// soap.genomics.org.cn) [40]. Then the distribution of reads on reference genome and genes was calculated and coverage analysis was performed. After alignment quality control, the DEGs were selected. And for further analysis, expression pattern analysis of DEGs were also performed.
To clarify the biological functions of DEGs, GO enrichment analysis of DEGs was performed as described above. To further clarify the biological functions of DEGs, KEGG pathway analysis was performed using the same calculating formula as GO enrichment analysis with database available at http://www.kegg.jp/kegg/.

Correlation analysis of DNA methylation and gene expression
The distribution of hyper-, hypo-and unmethylated gene expression levels in different gene elements were calculated to analyze the relationship between DNA methylation and gene expression as previously described [42].

Statistics and analysis
To identify true DMRs, the numbers of reads were calculated and assessed using chi-square statistics and False discovery rate (FDR) statistics (p≤0.01, and the difference of read numbers should be more than twice). To identify significantly enriched GO terms and KEGG pathways, gene numbers for every term or pathway were calculated and then assessed using hypergeometric test, p-value of hypergeometric test was corrected using Bonferroni Correction [43]. GO terms with corrected p-value ≤0.01, and KEGG pathways with corrected p-value ≤0.05 were considered significantly enriched. In selection of DEGs, the gene expression level was calculated using RPKM method [44], and DEGs were selected as previously described [45]. The adjusted p-value was calculated using Benjamini, Yekutieli. 2001 FDR method [46] and DEGs was defined as genes with FDR≤0.001 and the RPKM difference between the paired samples should be more than twice. Hierarchical cluster was performed to analyze the expression pattern of DEGs using Cluster [47] and presented using Java Treeview [48]. The DEGs were clustered by Euclidean distance.