Integrated Genome-Wide DNA Methylation and RNAseq Analysis of Hippocampal Specimens Identifies Potential Candidate Genes and Aberrant Signalling Pathways in Patients with Hippocampal Sclerosis
Correspondence Address: Source of Support: None, Conflict of Interest: None DOI: 10.4103/0028-3886.280649
Source of Support: None, Conflict of Interest: None
Keywords: DNA methylation, drug resistance, hippocampal sclerosis, gene expressionKey Messages: DNA methylation plays a crucial role in the regulation of the pathogenic mechanisms of epileptogenesis associated with HS.
Mesial temporal lobe epilepsy with hippocampal sclerosis (HS) is the most common type of drug resistant epilepsy (DRE) observed in adults, accounting for 40% of cases. The molecular pathogenesis and the mechanism of drug resistance in DRE pathologies like HS remain poorly understood and requires intensive investigations., DNA methylation mediated changes in gene expression has been implicated in experimental and human TLE. Increased DNA methyltransferase (DNMT) gene expression is reported in human TLE and inhibition of DNMTs in animal models were shown to reduce DNA methylation and epileptogenesis.,,, More recent studies reported genome-wide DNA methylation changes and correlation of differential methylation patterns of selected genes by quantitative PCR in tissue samples and recently on blood samples in human epilepsy.,, To evaluate the methylation regulated gene expression in HS patients, differential methylation pattern was correlated to the differential gene expression (DGE) analysis in this study. Simultaneously we have also analyzed the expression levels of the two DNMTs, DNMT1 and DNMT3a in the resected tissues. We have identified and discussed epigenetically deregulated genes and pathways with potential role in the pathophysiology of HS patients for the first time.
Patients included in the study
The patients who were diagnosed to have DRE due to HS and underwent surgery were included in this study [Table 1]. We have used histopathologically normal hippocampus tissues (data not shown) obtained from the post-mortem cases without any history of seizures or other neurological disorders as non-epileptic controls [Table 1]. All the autopsies were performed within 8 hours of death. The prospective study was reviewed and approved by Institutional Ethics Committee (IEC), All India Institute of Medical Sciences, New Delhi, India. Informed consent was obtained from all the participants. Tissue samples obtained from four HS patients (H1-H4) and two autopsy cases (A1-A2) were processed for DNA methylation. Two autopsy (A1-A2) samples were processed for RNAseq analysis. For HS, we included previously published data of RNA sequencing obtained from HS  patients (H1 and H2) for DGE analysis and correlation with methylation pattern.
An independent set of 10 patients (H5 to H14) and 6 non-epileptic controls (A3-A8) along with (H1-H4) and (A1-A2) included for validation of DNA methylation as well as DGE by qPCR analysis.
Genome wide DNA methylation analysis
Genomic DNA was extracted from brain tissue samples using genomic DNA extraction kit (Qiagen, Hilden, Germany) and methylated DNA was immunoprecipitated (MeDIP) as described previously. The immunoprecipitated DNA and a sample of input DNA were amplified and subsequently labeled using the SureTag DNA labeling kit (Agilent #5190-3400) and hybridized to Agilent Human DNA Methylation Microarray slides (G4495A, AMAIDID 023795) containing 237,227 biological probes covering all annotated 27,627 CpG islands and 5081 UMR regions with median probe spacing of 97 bp, following the manufacturer's instructions. The slides were scanned using Agilent scanner P/N G2565BA and methylation status analysis was carried out using GeneSpring GX software (version 13.0). Differentially methylated genes between the control and HS groups were identified using an adjusted Student's t-test (P < 0.05) corrected for multiple comparisons with the false discovery rate (FDR) and fold-change (FC) cut off of ≥2.
RNAseq methodology and analysis
For HS, previously published RNAseq was used in this analysis. RNAseq analysis for autopsy samples was carried out on Illumina platform as described previously. Briefly, total RNA was extracted from surgically resected tissues using the RNeasy Mini Kit (QIAGEN) according to the manufacturer's protocol, and treated with DNase I. RNA quality was confirmed to be high for all samples on the Agilent 2100 Bio-Analyzer (Agilent Technologies, Santa Clara, CA, USA), with an RNA integrity number (RIN) range of 8.1 to 10, and concentrations were determined using the Nanodrop ND-1000 (NanoDrop Technologies, Wilmington, DE, USA). Transcriptome sequencing libraries were prepared using the TruSeq RNA Access Library Prep Kit (Illumina, CA, USA). Paired-end sequencing (100 bp) was performed at an IlluminaHiSeq 2500 instrument. All the 100 bp paired-end raw reads were quality checked for low quality bases and adapter sequences. Quality check was done using NGS QC Toolkit (version 2.3.1). Low quality bases and reads were filtered out from the datasets prior to read mapping. Although, RNA seq analysis for HS and control were not performed simultaneously, they were performed on the same platform using same kit for transcriptome sequencing library (TruSeq RNA Access Library Prep Kit (Illumina, CA, USA).
The analysis was done using the raw data that was generated for the two groups on Illumina HiSeq 2500. After removal of the reads mapping to ribosomal and mitochondrial genome, the paired-end reads were aligned to the reference human genome, Feb. 2009 release, downloaded from UCSC database (GRCh37/hg19). Alignment was performed using HISAT (a fast spliced aligner with low memory requirements) program (version 0.1.7). Alignment QC which includes Splice junction and aligned reads distribution was performed using RSeQC (2.3.7) and RNA-SeQC (v1.1.8). The alignment percentage of the reads to the human genome was good albeit the differences in the read length of the two groups. The aligned reads were used for estimating expression of genes and transcripts using cufflinks program (version 2.2.1) and are reported as FPKM/RPKM (Fragment per kilo per million/Read per kilo per million) units for each of the genes and transcripts. Differential gene expression analysis was performed using cuffdiff (version 2.2.1). Significant mRNA fold-change was determined by an adjusted P value lower than 0.05 based on the Benjamini and Hochberg multiple testing correction.
Integrative analysis of DNA methylation and gene expression
To identify those DNA methylation changes with concomitant changes in gene expression, integrative analysis of DNA methyaltion and RNAseq data was carried out using gene spring software GX version 13.0 as described previously.
Functional classification/gene network analyses
Sets of genes with concordant DNA methylation and expression patterns were analyzed for gene Ontology (GO) enrichments using the DAVID (version 6.7) tool to cluster differentially-regulated genes by their common functionality [Table S1]. Pathway analysis was done using PANTHER) Classification System and analysis tools and network analysis using Natural language processing-based (NLP) network discovery algorithms in gene spring software version 13.0 as described previously.
Quantitative real time PCR (qPCR)
qPCR was performed to evaluate the expression of DNMT1 and DNMT3 α and validate the expression of differentially methylated and differentially expressed genes using specific primers [Table S2] for selected genes on 14 HS patients and 8 control samples. For MeDIP validation, previously published primers for GAPDH unmethylated and H19 imprinted were used as negative and positive controls. All samples were amplified in triplicates. In order to normalize qPCR reactions for gene expression analysis, HPRT (hypoxanthine phosphoribosyl-transferase) was included as reference gene. Real time PCR amplifications were performed in CFX 96 real time systems (Biorad) with Biorad CFX software manager with the following cycling parameters: an initial hot start of 95°C for 3 min followed by 40 cycles of 95°C for 5 s and 60°C for 30 s. The 2−ΔΔCq method was used to quantify the relative normalized expression of studied genes based on the average Ct values across samples (Schmittgen and Livak, 2008). Test of significance was performed using Student's t-test (two-tailed, unpaired) and P values < 0.05 were considered significant.
Differential DNA methylation and gene expression patterns in HS patients
DNA methylation profiles readily discriminated samples and controls into discrete groups [Figure 1]a. There are a total of 12886 hypermethylated and 3154 hypomethylated sites with adjusted P values < 0.05 and fold change (FC) ≥2 [Table 2]a. Out of these, 2205 genes were hypermethylated and 980 were hypomethylated in the promoter regions [Table 2]a. The distribution and location of the CpG sites with altered methylation are shown in [Figure 1]b. A high percentage (61%) of altered methylation sites are clustered within the gene body [Figure 1]b. Although several recent studies have begun to investigate the potential for how gene body DNA methylation affects gene expression, the function of gene body DNA methylation is yet not clearly understood. Since differential methylation of CpG-rich promoters is well accepted as an indicator of gene expression alterations in physiologic and pathologic conditions,,, we further filtered the gene list and selected genes with differential methylated promoters for further analysis. Enrichment of MeDIP-DNA was confirmed by a qPCR-based assay. The fold enrichment for H19 compared with GAPDH using the comparative Ct method ranged from 9 to 45-fold for each of the samples tested (data not shown). DNMT1 expression levels were significantly upregulated (P < 0.05) as compared to controls whereas DNMT3a mRNA levels were unaltered [Figure 1]c. DGE analysis revealed up regulation of 142 and down regulation of 102 genes in HS patients with log2 FC ≥2 and q-value < 0.05, when compared with the controls [Table 2]b.
Integrative analysis of DGE with the genes of differentially methylated promoters revealed significant number of genes with inverse corelation
102 genes showing down regulation were integrated with the 2205 hypermethylated genes and 142 genes showing up regulation were integrated with the 980 hypomethylated genes using the gene spring software GX version 13.0. A total of 66 DEGs genes were inversely correlating with the promoter methylation patterns. Promoters of 26 downregulated DEGs were hypermethylated and 40 upregulated DEGs were hypomethylated [Figure 2]a and [Figure 2]b suggesting that abnormal DNA methylation might have functional consequences. Fold change expression profiles of the differentially expressed genes, as observed in qPCR correlated with the RNAseq data [Figure 2]c.
Functional gene clustering, pathway and network analysis of inversely correlated genes revealed epigeneticaly regulated canonical pathways associated with HS
We limited further analysis to these 66 genes showing inverse correlation of their promoter methylation and gene expression patterns [Table S3]. The most important signalling pathways identified were axon guidance mediated by semaphorins, ionotropic glutamate receptor pathway and notch signaling pathways [Table S4].
Gene clustering and network analysis also identified various genes from these pathways like TFAP2A, NRP-1, SEMA3B, CACNG2, IGFBP2, SLK, TUBB6, IRF6, SCNN1A, MAP3K11, RAB13 and ADAM17 providing important insight into the understanding mechanisms underlying epigenetic regulation, via DNA methylation [Figure 3], [Figure S1] and [Figure S2]. Further classification of these modulated genes were done based on their molecular functions [Table S5].
Methylation array and RNA-Seq techniques are powerful tools to study global methylation variation and transcription changes, but there are no reports on integrative analysis of the two data in MTLE-HS patients. Here, we are discussing the role of the genes which showed an inverse correlation between methylation and gene expression, and their possible association with epileptogenesis and/or drug resistance in HS by grouping them in the three hubs that were most prominently overrepresented in the network analysis, pathway analysis and cluster analysis [Table S5].
Genes involved in neuronal development and cell-cell interactions
We found that in HS patients, hypermethylation of the IRF6, TJP3, PCDH3 and MAP1S promoter correlated with reduced gene expression and hypomethylation of the SEMA3B, NRP1, IGF2BP2, TFAP2A, TUBB6, CBFB, CAP1, TUBB6 and CNN2 promoters correlated with increased expression. This result suggests that the expression of these genes were regulated by DNA methylation.
Semaphorins have emerged as important regulators of excitatory and inhibitory synaptogenesis and function. Indeed, Sema3F, Sema 3A and Sema4D influence neuronal activity in animal models of epilepsy., However, there are no reports available regarding SEMA3B in epilepsy. Since, SEMA3B-NRP1 mediated immune response and apoptosis is reported, their involvement in neuroinflammation and cell death cannot be ruled out in epileptic condition. NRP-1 mediates the chemorepulsant activity of semaphorins.,,
Mossy fiber sprouting in HS relies largely on proper arrangement of cytoskeletal components, such as microtubule-associated proteins.,, Interestingly, in our study, cytoskeletal-related genes CAP1, CNN2 and TUBB6 expression was found to be upregulated while MAP1S were downregulated. Expression of TJP3 (Tight junction protein 3) and PCDH3 (Protocadherin Beta Cluster) were decreased in our study. These genes are involved in the mechanism regulating cell-cell interaction.
Several transcription factors and other genes involved in transcription/neurodevelopmental processes are found to be differentially expressed, i.e., TFAP2A, CBFB, ARID1A, IRF6 and IRX2. In our network analysis and pathway analysis, transcription factor AP-2 alpha (TFAP2A) emerged as an important molecule found to be upregulated. However, no reports are available in epilepsy. Defects in this gene are a cause of branchiooculofacial syndrome (BOFS). Reduced, or overexpressed TFAP2A expression is often detected in several types of cancers. TFAP2 α induces p21 dependent P53 expression corroborating the observed ability of TFAP2α to induce G1 and G2 cell cycle arrest., Elevated p53 expression in the hippocampus from patients with intractable temporal lobe epilepsy was reported previously. Upregulation of Bcl-2 family proteins in both resected TLE samples and experimental models of seizure-induced neuronal death have been reported earlier. Serum Bcl-2 levels are elevated in children with TLE and is further increased in drug refractory patients. Taking these account into consideration, the role of TFAP2A need to be explored in HS. Mutations in IRF6 gene can cause van der Woude syndrome and popliteal pterygium syndrome. Mutations in this gene are also associated with non-syndromicorofacial cleft type. Children of mothers with epilepsy may suffer from facial clefts more frequently than children of non-epileptic mothers.
Ion channel dysfunction
Ion channel dysfunction represents a fundamental pathophysiological element underlying abnormal neuronal discharge in seizure disorders., CACNG2 and SCNN1A showed hypomethylation and upregulated expression in the HS. CACNG2 constitutes on the major component of ionotropic glutamate receptor pathway regulate both trafficking and channel gating of the AMPA receptors. A mutation affecting affecting CACNG2 (TARP γ-2) has been described in an individual with moderate intellectual disability (ID). Mutations in this gene cause an autosomal dominant form of cognitive disability. Defects in CACNG2 s hown to cause seizures in mice. CACNG2 could be a potential target in HS.,,
Protein kinases and cell signaling
In our study, ADAM17 is an important factor of notch signaling is found to be differentially methylated with increased expression. Sha et al. demonstrated a correlation between aberrant Notch signaling and epileptic seizures. Downregulation of ADAM17, as observed in our study, could represent a homeostatic effect in response to seizure. Apart from Notch 1 cleavage to NICD1, ADAM17 is also important for the release of many cytokines and its receptors, growth factors and some cell adhesion molecules.,
MAP3K11 involved in MAPK signaling and CCK signaling was downregulated. The MAPK pathway has been implicated in epileptic seizures due to hippocampal sclerosis., Apart from that, MAP3K11 plays an important role in CCKR (Cholecystokinin receptor) signaling pathway. CCK controls the perisomatically targeting inhibitory interneurons of the hippocampus., CCK may have anticonvulsant and neuroprotective properties. Pretreatment of hippocampal slices with sulphated CCK blocked the effect of kainic acid on synaptic transmission. Seizures induced by picrotoxin and electroshock are inhibited by intracerebroventricular administration of CCK in rats. Therefore, MAP3K11 which is downregulated in our study is another potential therapeutic agent that should be further studied.
In the present study, expression of EPS8, RAB13, TTBK1 and DIRAS1, were downregulated while expression of SLK and IGF2BP-2 was upregulated. EPS8 controls dendritic spine density and synaptic plasticity through its actin-capping activity. RAB proteins are involved in the maintenance of neuronal vesicular trafficking in the central nervous system. SLK is a microtubule-associated protein inducing actin stress fiber disassembly, it also mediates apoptosis 
TTBK1 transgenic mice showed enhanced activity of CDK5 and GSK3 β, and accumulation of hyperphosphorylated tau in the brain. Hyperphosphorylated tau in patients with refractory epilepsy correlates with cognitive decline. However in our study its expression has been found to be downregulated. Further study on larger sample size is needed to clarify its function in HS.
Our study has several limitations. The first limitation of this study was the small sample size. However, we have validated our MeDIP and RNAseq data on a larger cohort but the sample size for MeDIP array and RNAseq was limited. There are no ideal or acceptable non-seizure controls for such studies involving epilepsy surgery. Owing to several limitations associated with tumor and trauma tissues, we used histologically normal hippocampus autopsy samples as controls.
Furthermore, DNA methylation is dynamically regulated and influenced by a variety of factors, including gender, age, and social environment. Methylation patterns also differ between brain regions and neuronal cell types., We selected the hippocampal tissues from autopsy as well as patients for DNA methylation analysis to reduce interference factors. In addition, all patients were refractory to anti-epileptic drugs. The effects of the drugs on DNA methylation cannot be excluded. Also, since the role of other CpG islands in epigenetic regulation of gene expression is not well-defined, we have only discussed the CpG islands in the upstream promoter regions. Future experimentations addressing these limitations are likely to provide us with a better knowledge of the underlying mechanisms and pharmacological targets of HS.
In summary, the present study presents a preliminary set of data indicating the role of DNA methylation as epigenetic regulator of gene expression in patients with HS. Based on bioinformatics analysis, some of the differentially methylated sites are associated with promoters, genes and molecular pathways important for development, axon guidance, transcription regulation and protein phosphorylation, ion channel activity, cytoskeletal reorganisation and neuroinflammation with most distinct ones included SEMA3B, NRP-1, TFAP2A, ADAM17, CACNG2 and TTBK1.
Declaration of patient consent
The authors certify that they have obtained all appropriate patient consent forms. In the form, the patient(s) has/have given his/her/their consent for his/her/their images and other clinical information to be reported in the journal. The patients understand that their names and initials will not be published and due efforts will be made to conceal their identity, but anonymity cannot be guaranteed.
Authors would like to thank all the patients and their legal guardians for participating in this study. This work was supported by financial grants to the Centre of Excellence for Epilepsy, a collaborative project between National Brain Research Centre (NBRC), Manesar, and All India Institute of Medical Sciences (AIIMS), New Delhi, by Department of Biotechnology, Ministry of Science and Technology, Government of India (#BT/MED/122/SP24580/2018), MG grant and DST-PURSE grant, University of Delhi, Delhi.
Financial support and sponsorship
Conflicts of interest
There are no conflicts of interest.
[Figure 1], [Figure 2], [Figure 3]
[Table 1], [Table 2]