Abstract
The identification of mechanisms transforming normal to seizure-generating tissue after brain injury is key to developing new antiepileptogenic treatments. MicroRNAs (miRNAs) may act as regulators and potential treatment targets for epileptogenesis. Here, we undertook a meta-analysis of changes in miRNA expression in the hippocampal dentate gyrus (DG) following an epileptogenic insult in three epilepsy models. We identified 26 miRNAs significantly differentially expressed during epileptogenesis, and five differentially expressed in chronic epilepsy. Of these, 13 were not identified in any of the individual studies. To assess the role of these miRNAs, we predicted their mRNA targets and then filtered the list to include only target genes expressed in DG and negatively correlated with miRNA expression. Functional enrichment analysis of mRNA targets of miRNAs dysregulated during epileptogenesis suggested a role for molecular processes related to inflammation and synaptic function. Our results identify new miRNAs associated with epileptogenesis from existing data, highlighting the utility of meta-analysis in maximizing value from preclinical data.
Significance Statement
Meta-analyses of data from human research studies are an invaluable tool, and the methods to conduct these investigations are well established. However, meta-analyses of preclinical data are rarely undertaken, due to the typically small sample sizes and the substantial heterogeneity between studies. We implemented a meta-analysis of microRNA (miRNA) expression changes in animal studies of epilepsy. This is the first study of its kind in the field of epilepsy and one of the first in preclinical research. Our analyses identify new miRNAs associated with epileptogenesis and epilepsy, highlighting common mechanisms across different animal models. These miRNAs and their predicted effects on gene expression generate new hypotheses about the causes of epilepsy that will prompt new studies in the field.
Introduction
Epilepsy is a serious, common neurologic disorder primarily characterized by the occurrence of spontaneous seizures. The treatment of epilepsy remains as one of the major unmet medical needs in neurology because, despite of over 20 antiepileptic drugs on the market, seizures are not controlled in about one third of the patients. The most common form of epilepsy in adults originates in temporal structures of the brain (temporal lobe epilepsy, TLE; Hauser et al., 1993). Epilepsy (TLE in particular) frequently arises as a consequence of brain injury (“acquired epilepsy”). While acquired epilepsies are in principle preventable by the therapeutic targeting of molecular processes underpinning their development (i.e., antiepileptogenic therapies), there are currently no established treatment options for halting the transformation of normal brain tissue to epileptic (Simonato et al., 2014). Identification of the mechanisms underlying epileptogenesis would therefore facilitate the identification of therapies for preventing the development of epilepsy and may inform new strategies for overcoming drug resistance in epilepsy more generally (Simonato et al., 2012; 2013; 2014).
Recent studies have suggested that microRNAs (miRNAs) play an important role in the pathogenesis of acquired epilepsy and may represent novel therapeutic targets (Jimenez-Mateos et al., 2012; Brennan et al., 2016; for review, see Cattani et al., 2016; Henshall et al., 2016). miRNAs are a family of small (21-25 nucleotides) noncoding RNAs, which can modulate various cellular and biological processes by degrading or repressing translation of specific mRNAs (Bartel, 2004, Guo et al., 2010). In systems analysis, miRNAs and their gene targets are described as following a “many-to-many” data model, such that each miRNA may regulate many transcripts and a single transcript may be regulated by many miRNAs (Ebert and Sharp, 2012). miRNAs have been implicated in various neuronal functions that are relevant in the pathogenesis of neurologic diseases, including epilepsy (Tan et al., 2013; Rajman et al., 2017; for review, see Karnati et al., 2015).
Interpretation of data investigating the dysregulation of miRNAs in the brains of patients with epilepsy is challenged by the absence of appropriate human control brain tissue (Roncon et al., 2016). Research on the role of miRNAs in epilepsy has therefore focused on the use of experimental models of epilepsy, revealing changes in the expression of hippocampal miRNAs at different stages of the epileptogenic process (Henshall et al., 2016). However, methodological differences between the various preclinical animal models of epilepsy have made comparisons between studies difficult and the identification of common pathways dysregulated in epileptogenesis and epilepsy problematic. The current preclinical miRNA studies vary in multiple parameters, including brain region analyzed, animal model, sample size, microarray platform and analysis technique. Moreover, these studies are generally substantially underpowered to reliably detect modest changes in miRNA expression.
One particular concern is that analysis of large brain regions (like the hippocampus or the cortex) across different studies may confound interpretation and comparison because of variable cellular composition (e.g., relating to variable neuronal loss, astrocytosis, microgliosis, etc.). One way to address this issue could be to focus on a specific cell population. In this respect, dentate gyrus (DG) granule cells (GCs) seem particularly attractive as a target of analysis as the DG GC layer is a compact layer of (almost) identical cells, facilitating the dissection of a nearly pure cell population (GCs). Moreover, the DG has been traditionally described as a “gate” to inhibit hippocampal overexcitation (Chevaleyre and Siegelbaum, 2010) and recently, this hypothesis found support from new technologies; optogenetic GC hyperpolarization was found to stop spontaneous seizures, whereas optogenetic GC activation exacerbated spontaneous seizures, and activating GCs in nonepileptic animals evoked acute seizures (Krook-Magnuson et al., 2015). Finally, the DG is known to undergo important functional changes during epileptogenesis (neurogenesis, mossy fiber sprouting, increased excitation; Pitkänen and Lukasiuk, 2011).
To date, three studies have investigated differential expression of miRNAs in the DG during the epileptogenesis and the chronic phase of epilepsy in rats (Bot et al., 2013; Gorter et al., 2014; Roncon et al., 2015). Each of these studies used a different method to trigger epileptogenesis: focal electrical stimulation of the lateral nucleus of the amygdala (Bot et al., 2013); focal electrical stimulation of the angular bundle, a major afferent pathway from the entorhinal cortex to the hippocampus (Gorter et al., 2014); and the systemically administered chemoconvulsant pilocarpine (Roncon et al., 2015). Interestingly, all these models imply a key involvement of the DG in the development of epilepsy but via a different epileptogenic insult: direct activation in the case of angular bundle stimulation, indirect in amygdala stimulation, and a widespread brain activation in the case of pilocarpine (Peng and Houser, 2005).
Here, we aimed to overcome some of the limitations related to individual studies by combining these three studies in a meta-analysis, the aim being to increase the statistical power for detecting differentially expressed miRNAs while accounting for study heterogeneity, ultimately leading to more robust and accurate predictions of dysregulated downstream pathways (Ramasamy et al., 2008; Yang et al., 2014). Moreover, this approach offers the opportunity to identify miRNA changes that are independent of the model of epilepsy, i.e., more likely to be disease rather than model related. Our meta-analysis was performed at two time points in the “natural history” of the experimental disease: epileptogenesis and the chronic phase of epilepsy.
Materials and Methods
Inclusion criteria and study design
We collected datasets for meta-analysis based on available genome-wide expression profiles of miRNAs from the DG from epileptic and control hippocampi during epileptogenesis and chronic epilepsy. To assist the functional inference of differentially expressed miRNAs we analyzed publicly available gene expression data obtained from the DG of epileptic rodents.
To identify relevant studies, we first undertook a systematic search to identify all published studies of miRNA expression levels and/or gene expression between cases (epileptic) and controls, in the DG of animal models of epilepsy. We conducted a PubMed search based on the string: “(microRNA OR miRNA) AND (dentate gyrus OR dentate cells OR granule cells) AND epilepsy.” miRNAs or genes expression profiles data obtained from the whole hippocampi or different brain regions to the DG were not included in the meta-analysis. This search and inclusion criteria identified only three relevant articles (Bot et al., 2013; Gorter et al., 2014; Roncon et al., 2015).
Time points for each stage of epileptogenesis and for the chronic phase of epilepsy have not been standardized. For the purposes of this study, we considered 7-8 d after SE as the “epileptogenesis phase” and more than two months after SE as the “chronic phase.”
The following models were used in the three relevant papers and considered for the meta-analysis. (1) Pilocarpine model: a microarray study based on the investigation of miRNAs differentially expressed in the laser-microdissected GC layer of the DG of the hippocampi of pilocarpine-induced epileptic rats and matched controls (n = 4), killed during the late phase of latency, 8 d after SE (n = 5) and in the chronic stage, 50 d after the first spontaneous seizure (n = 5; Roncon et al., 2015). (2) Amygdala stimulation: a microarray study focused on the differential expression of miRNAs and genes in hand-dissected DG in the amygdala stimulation rat model during the phase of epileptogenesis, 7 d after the stimulation (n = 5), in the chronic stage, 60 d after the stimulation (n = 5) and controls (n = 5; Bot et al., 2013). (3) Perforant path stimulation: a microarray miRNA study based on the perforant path stimulation rat model of epilepsy. We analyzed DG samples obtained from stimulated and control rat hippocampi, during latency (8 d after SE; n = 8) and three to four months after the stimulation for the chronic stage (n = 6) and controls (n = 10; Gorter et al., 2014). The datasets considered for the meta-analysis are summarized in Table 1.
Datasets included in the meta-analysis
Power calculation
The statistical power of cases and controls for each individual model was calculated using pooled SD of each expressed miRNA. Power to detect miRNAs differentially expressed at multiple fold changes (1, 1.5, 2, 2.5, 3, 3.5, and 4) were calculated, considering miRNA expression variability ranging from 20th, 40th, 60th, and 80th to 100th percentile of the respective SD profiles per model (Fig. 1). Power calculations were performed using R bioconductor package ‘pwr’ version 1.2.
Data processing
From each identified study the following information was extracted: platform, number of cases and controls, and miRNA expression data at different time points of the disease. When available, GEO accession number and gene expression data were extracted (Table 1).
Data transformation
Since different platforms had been used to generate miRNA expression values, a linear transformation approach was applied to each miRNA using a Z-score transformation according to the formula:
Where Xi is the normalized intensity data for each miRNA, X is the average normalized miRNA intensity within a single study, and δ is the SD of cases and controls within respective studies.
Effect size estimation
A meta-analysis is “a technique for quantitatively combining and integrating the results of multiple studies on a given topic” (Polit and Beck, 2004). Thus, a key aspect of meta-analysis is to measure differences and direction of change from quantitative research studies (Polit and Beck, 2004; Berben et al., 2012). A common metric used to provide this important information is the effect size calculation. Accordingly, to give a statistical expression of the magnitude of the difference between groups (i.e., epileptogenesis vs controls and chronic stage vs controls) in regard of miRNAs expression, we estimated the effect size of each individual study defined as the standardized mean difference (SMD) between cases and controls. The SMD has been calculated using the Hedge’s method with the following formula:
Where X1 is the mean of cases, X2 the mean of the control group, and δ is the SD.
Statistical heterogeneity
Different animal models, tissue collection methods, and platforms were used to generate the datasets. This makes it difficult to directly compare the data, the risk being of skewing comparison results, reducing the reliability of measurements of individual miRNA expression changes (Yang et al., 2014). Statistical heterogeneity was assessed using Cochrane meta-regression approach calculated by Q test, I2 statistics, and Tau2 statistics. These measures were applied at each dataset to assess the overall heterogeneity (Higgins et al., 2003; Ioannidis et al., 2007). To test the total variance of each miRNA within the studies, the Cochran Q test have been run, according to the formula:
Where k is the number of the studies included in the meta-analysis, T is the number of variables observed, b is the number of miRNAs included in the test, and N is the total number. A Benjamini-Hochberg (BH) adjusted p < 0.05 was considered statistically significant for the Cochrane Q test.
Furthermore, I2 statistic has been employed to describe the percentage of the variability in the effect size estimates, following the formula:
Where Q is the value derived from the Cochran Q test and df are the degrees of freedom. Tentatively, I2 statistic can be considered as an indicator of heterogeneity, where low, moderate, and high heterogeneity corresponds to I2 values of 0–0.3, 0.3–0.7, and 0.7–1, respectively (Higgins et al., 2003).
Finally, to estimate the variance across studies the Tau2 has been applied:
Where Q is the value derived from the Cochran Q test, df are the degrees of freedom and C is a scaling factor which takes into consideration that the Q value is the weighted sum of squares.
Meta-analysis
The presence of statistical heterogeneity among the studies led us to use a random-effect model for the meta-analysis rather than a fixed-effect model. The pooled effect size (PES) for each miRNA was obtained applying the random effect size model based on the DerSimonian and Laird method (DerSimonian and Laird, 1986; DerSimonian and Kacker, 2007). We generated one forest plot for each miRNA in both the epileptogenesis and the chronic phase to depict the SMD along with its 95% confidence interval (95% CI) for individual studies as well as the pooled MD from the meta-analysis.
miRNAs:mRNAs inverse-fold change
miRNAs predominantly act by repressing their target genes by decreasing target mRNA levels (Bartel, 2004; Guo et al., 2010). Therefore, we investigated the correlation of mRNAs predicted targets expression in the available transcriptomic datasets.
To predict miRNA target genes, the miRNA-target interactions were analyzed with the web-based tool miRwalk (Dweep et al., 2011; Dweep et al., 2014). The rule that the 5’ region of miRNA from nucleotides 2-8 (“seed region”) has importance in targeting, is commonly accepted as the canonical mechanism by which miRNAs complementary convey functional binding to mRNA targets (Jinek and Doudna, 2009). However, despite the importance of the seed region, the 3’ end of a mRNA also contributes to the binding in ∼2% of all preferentially conserved sites (Grimson et al., 2007; Shkumatava et al., 2009). In addition, some validated miRNAs have a binding site that exhibits few continuous base pairs in the control region (Shin et al., 2010). Thus, to figure out the complete mechanism of miRNA regulation, we expanded the miRNAs-binding site prediction within the 3’, 5’ untranslated regions (UTRs), and the seed sequence, with a minimum seed length of seven nucleotides. Furthermore, to exclude overprediction we applied a comparative analysis by six prediction programs: miRMap, RNA22, miRanda, RNAhybrid, PICTAR2, and Targetscan. Following this approach, a candidate mRNA target has to be identified by all these programs. As we did for the miRNAs, we conducted PubMed search based on the keywords: “gene expression, epilepsy, dentate gyrus.” Two studies were identified that were then included in our analysis (Table 1): Dingledine et al. (2017; GEO repository accession number: GSE47752) and Bot et al. (2013; GEO repository accession number: GSE49850). The first includes gene expression data obtained from laser-microdissected DG of rats that received different SE models of epilepsy: pilocarpine, self-sustained SE (SSSE), and kainite, it also included kindling, not considered in this study not being a post-SE model. We considered in our analysis only those rats killed 10 d after SE that did not develop spontaneous seizures, as the best time point matching the epileptogenesis phase used for the miRNAs meta-analysis (7-8 d after SE). In Dingledine et al. (2017), the kainate and pilocarpine models were performed in two independent labs, while the SSSE model that was performed in only one lab. For studies undertaken in multiple labs, we combined the p values (obtained from differential expression analysis) from the independent labs using Fisher’s method. The second dataset (GSE49850) was obtained from the same amygdala-stimulated rats used for the miRNA analysis.
To investigate miRNA-mRNA interactions, we included only those mRNAs that had inverse relationship to miRNA changes in at least three datasets in the epileptogenesis phase; while for the chronic phase, we considered all predicted mRNAs that presented an inverse relationship to miRNA changes in the amygdala stimulation dataset only.
miRNA functional enrichment analysis
Functional enrichment analysis using gene ontology (GO) and pathways enrichment analysis based on the Kyoto encyclopedia of genes and genomes (KEGG) database were performed using Webgestalt webserver (Zhang et al., 2005; Wang et al., 2013). The enrichment was performed with a hypergeometric test separately for the list of predicted targets based on those miRNAs dysregulated in epileptogenesis and in the chronic stage. Significant canonical pathways maps were selected according to a false discovery rate (FDR) < 5%.
To infer functional relationships between miRNAs identified using meta-analysis, a network of miRNAs based on their ability to target common pathways has been generated. A connection was made between a pair of miRNA, if respective mRNA targets belonged to the same pathways or GO terms that were significantly (FDR < 0.05) enriched for combined miRNA targets. The network was visualizes using Cytoscape version 2.8.2.
Results
Study design
We included in the meta-analysis all published miRNA expression datasets from dissected DG of the hippocampus that compared control (baseline) tissue with tissue from epileptogenesis and chronically epileptic rats (Fig. 2A). Based on these inclusion criteria, we identified three datasets that used different epilepsy models (Table 1): (1) pilocarpine (Roncon et al., 2015), (2) amygdala stimulation (Bot et al., 2013), (3) angular bundle stimulation (Gorter et al., 2014). We first calculated the power of each individual study to detect significant changes in miRNA expression and found that all three individual studies were substantially underpowered to detect modest fold changes (<2.0) in miRNA expression (Fig. 1). Of the total number of miRNAs expressed at any time point in any of the three models, the expression levels of 176 miRNAs were detectable across all three studies (Fig. 2B). The expression values of these 176 miRNAs were then Z-score transformed and meta-analyzed across the three studies.
Power calculation. Power calculation is plotted as the power (y-axis) to detect a miRNA with fold change (x-axis) according to the percentile of the ranked SDs for miRNAs for each study. Across all three models, the power to detect miRNA with fold change 2 or less falls below 80% for at least 40% of the miRNAs.
Study design and data preprocessing. A, Study design. B, Venn diagram showing miRNAs commonly expressed between the three studies included in the meta-analysis. C, Statistical heterogeneity estimation. I2 scores of commonly expressed miRNAs (n = 176) in epileptogenesis and in chronic stage. miRNAs are ordered based on the adjusted p value after meta-analysis. I2 < 0.3, low heterogeneity; 0.3 < I2 > 0.7, moderate heterogeneity; I2 > 0.7, high heterogeneity. SRS, spontaneous recurrent seizures; pilo, pilocarpine model; amy stim, amygdala stimulation model; AB stim, angular bundle stimulation model.
Estimation of statistical heterogeneity
Meta-regression analyses were performed separately for epileptogenesis and chronic stages of epilepsy for the 176 miRNAs that were expressed in all three datasets (Fig. 2B). There are two models that are commonly used to perform meta-analysis, the fixed effect and the random effects models. The fixed effect model assumes that the effect size is the same in all studies, while the random effects analysis assumes that the effect can vary from one study to another. To determine the correct model for this study, we first estimated statistical heterogeneity using the Cochrane’s Q test, separately for the epileptogenesis and the chronic phase. This revealed significant heterogeneity between studies at both stages of the disease (BH, adjusted p < 0.05). Next, to assess the proportion of miRNAs that were differentially expressed between studies, during epileptogenesis and the chronic stage separately, we calculated the I2 statistics (Higgins et al., 2003). Of the 176 miRNAs measured across the three studies, 26.14% revealed a high level of heterogeneity: 31.25% a moderate level and 42.61% a low rate of heterogeneity during epileptogenesis (Fig. 2C). In the chronic stage, 18.75% displayed high, 30.11% moderate, and 51.14% low level heterogeneity (Fig. 2C). Collectively, these observations favor the use of the random effects model.
Differential expression of miRNAs in the epileptic DG
Using a random effects meta-analysis of miRNA changes in the three models of epileptogenesis and adopting a stringent correction for multiple testing to minimize false positives (Bonferroni adjusted p < 0.05), we identified 26 and 5 differentially expressed miRNAs between control and latency and between control and chronic epilepsy, respectively. Full results of all these miRNAs including PES estimations, I2, Tau2 and p values are shown in Tables 2, 3. Forest plots for selected miRNA are shown in Figure 3. Comparing these results with those presented in each original studies that have been meta-analyzed here, our meta-analysis identified 11 miRNAs differentially expressed in epileptogenesis compared to control and two miRNAs (i.e., miR-324-3p and miR-130a-3p) in the chronic stage of epilepsy that were not identified as significantly differentially expressed in any of the individual studies (Tables 2, 3, miRNAs highlighted in bold). The datasets employed in this meta-analysis do not allow a precise evaluation of the abundance of expression of these 26 plus five miRNAs under control conditions, but a relative abundance estimate based on the internal standards employed in each study suggests that almost all are expressed at medium to high abundance in the control DG (only miR-212-5p was expressed at relatively low levels but upregulated during epileptogenesis).
Differentially expressed microRNAs in the epileptogenesis period compared with controls
Differentially expressed microRNAs in the chronic stage compared with controls
Forest plots of selected miRNA. Forest plots for miR-7a-5p, miR-92b-3p, miR-101a-3p, miR-138-5p, miR-150-3p, miR-153-3p, miR-335, miR-383-3p, and miR-3573-3p are shown for the phase of epileptogenesis, and miR-130a-3p and miR-148b-3p for the chronic stage. For each miRNA, the effect size of the individual studies is reported as MD and 95% CI. The % weight refers to random effects analysis. Individual effect sizes are represented by colored boxes (green for epileptogenesis and blue for the chronic period) and 95% CI are denoted by black lines. The combined effect sizes are represented by diamonds, where diamond width correspond to the 95% CI bounds; boxes and diamonds size is proportional to effect size estimation precision. For each miRNA, the weight of the dataset in the combined analysis has been reported in percentage.
Relationship between miRNAs and mRNAs expression changes
Previous transcriptomic studies in epilepsy models have revealed dysregulation of many genes in the different phases of the experimental disease (Lukasiuk and Pitkänen, 2004). It can be hypothesized that differentially expressed miRNAs may contribute to these alterations (upregulated miRNAs may downregulate their mRNA targets whereas downregulation of miRNAs may allow upregulation of their mRNA targets). To assess the role of differentially expressed miRNAs during the course of epilepsy, that is, to infer their mRNA regulatory targets, we explored the potential regulatory relationship between miRNA and mRNA changes.
First, to predict miRNA target transcripts we used the miRWalk database (Dweep et al., 2011; Dweep et al., 2014), which combines information across six miRNA target prediction programs (miRMap, RNA22, miRanda, RNAhybrid, PITCAR2 and Targetscan). As expected, this analysis identified a very large number of predicted targets but, obviously, the large majority of these may not be expressed in the DG or may not be expressed in a negatively correlated fashion relative to the miRNAs. Therefore, we asked whether the miRNAs that were identified by the meta-analysis as dysregulated in epileptogenesis and in the chronic stage were negatively correlated with changes in DG expression in their predicted mRNA targets. To this end, we took advantage of publicly available gene expression data generated in three separate datasets that investigated mRNA expression changes in the DG of rats from the epilepsy models used in our meta-analysis: pilocarpine (GEO accession: GSE47752; Dingledine et al., 2017), angular bundle stimulation (called SSSE in this database; GEO accession: GSE47752), and amygdala stimulation (Bot et al., 2013; PMID: 24146813). Only the amygdala stimulation mRNA dataset (Bot et al., 2013) was obtained from the same animals employed for obtaining the miRNA dataset and included data on the chronic phase. The other datasets were generated by separate research groups on a separate group of animals, the experimental procedures were slightly different from those used in the corresponding miRNAs studies (Gorter et al., 2014; Roncon et al., 2015) and included mRNA data for these models related to the epileptogenesis phase only. In addition, Dingledine et al. (2017) also included datasets for another SE model (kainate) that we also considered in our analysis.
Analysis of the epileptogenesis data identified inverse relationship, based on significant fold changes (gene FDR < 0.1), between 22 (of 26) miRNAs and 122 unique predicted gene targets in at least three of the four epilepsy models of the Dingledine dataset (Dingledine et al., 2017) that we considered in this study (Table 4). The mRNAs encoding the mitogen-activated protein kinase kinase kinase 4 (Map3k4) and the enhancer of zeste 2 polycomb repressive complex 2 subunit (Ezh2) were predicted targets and had inverse expression relative to four miRNAs, i.e., miR-92b-3p, miR-101a-3p, miR-153-3p, and miR-3575-3p for Map3k4 and miR-92b-3p, miR-101a-3p, miR-138-5p, and miR-153-3p for Ezh2; synapsin type 2 was inversely correlated to three miRNAs (miR-101a-3p, miR-139-5p, and miR-551b-3p); the mitogen-activated protein kinase kinase kinase 14 (Map3k14) and the protein tyrosine phosphatase, nonreceptor type 5 (Ptpn5) were inversely correlated with two miRNAs, i.e., miR-7a-5p and miR-138-5p for Map3k14; miR-150-5p and miR-383-5p for Ptpn5. All the above transcripts were upregulated during epileptogenesis and predicted targets of miRNAs that were downregulated. In contrast, miR-132-3p, miR-146a-5p, miR-212-3p, and miR-212-5p were upregulated during epileptogenesis. The 5-hydroxytryptamine receptor 5 (Htr5b) and the β-1,3-galactosyltransferase 5 (B3galt5) were predicted targets of miR-146a-5p and were downregulated. The γ-aminobutiric acid receptor subunit δ (Gabrd) was a predicted target of and inversely correlated with miR-212-5p. Finally, we observed that miR-344b-2-3p, let-7d-3p, miR-21-5p, miR-29c-5p, and miR-324-5p were not anticorrelated with any of their predicted mRNA targets. Representative graphs for the anticorrelations between miRNAs and predicted mRNA targets are shown in Figure 4.
miRNA-mRNA fold changes inverse correlation in epileptogenesis
Relationship between selected miRNA and their predicted targets in different model of TLE. All panels show selected miRNAs-mRNA anticorrelation based on miRNAs and mRNAs fold changes in epileptogenesis. A, Inverse relationship between four downregulated miRNAs (miR-92b-3p, miR-101a-3p, miR-153-3p, and miR-3573-3p) and the commonly predicted target Map3k4. B, Inverse relationship between the downregulated miR-138-5p, miR-7a-5p, and the upregulated Map3k14. C, Inverse relationship between miR-101a-3p, miR-139-3p, miR-551b-3p, and Syn2. D, Inverse relationship between miR-150-5p, miR-383-5p, and Ptpn. E, F, Examples of the opposite anticorrelation, the upregulated miR-146a-5p with the downregulated Htr5b transcript, and the upregulated miR-212-5p and the downregulated Gabrd transcript.
The relationship between the changes in expression of miRNAs and their mRNA targets in the chronic stage of epilepsy was analyzed using only the amygdala stimulation dataset (Bot et al., 2013). We observed negative correlations (based on fold changes) between all five miRNAs that emerged as significantly downregulated from the meta-analysis and 29 unique predicted mRNA targets in the dataset. Five of these 29 anticorrelated mRNAs were predicted targets and had inverse expression relative to two miRNAs and one, the glutamate ionotropic receptor δ-type subunit 2 (Grid2), was a predicted target and had inverse expression relative to three miRNAs, namely, miR-130a-3p, miR-148b-3p, and miR-551b-3p (Table 5). Furthermore, interestingly, three mRNA targets, the transmembrane protein 176B (Tmem176b), the EFR3 homolog A (Efr3a), and the zinc finger, MIZ-type containing 1 (Zmiz1) were downregulated in both epileptogenesis and the chronic stage.
miRNA-mRNA fold changes inverse correlation at the chronic stage
Recent evidence supports the notion that miRNAs not only decrease levels of their mRNA targets (Guo et al., 2010), but additionally may have nuclear functions capable of influencing gene expression, and which may be reflected by a correlation between a miRNA and its target gene mRNA levels (Catalanotto et al., 2016). Analysis of the epileptogenesis data revealed significant correlation (gene FDR < 0.1), for 21 (of 26) miRNAs and 77 unique predicted gene targets in at least three of the four epilepsy models of the Dingledine dataset (Dingledine et al., 2017; Table 6). In addition, we found positive correlations between five of the five miRNAs that were downregulated in the chronic period and 39 predicted mRNA targets in the amygdala stimulation dataset (Bot et al., 2013; Table 7). Interestingly, 29 of the mRNAs identified as potential targets in epileptogenesis were inversely correlated to some miRNAs and directly correlated to others (e.g., map3k14 is inversely correlated to miR-7a-5p and miR-138-5p and directly correlated to miR-212-5p, while bdnf is inversely correlated to let-7b-3p and directly correlated to miR-212-5p). This observation prompts the hypothesis that some mRNAs may be subject to a dual control by different miRNAs at cytosolic and nuclear level. This hypothesis should be challenged and investigated.
miRNA-mRNA fold changes positive correlation in epileptogenesis
miRNA-mRNA fold changes positive correlation in the chronic period
Genes that are anticorrelated with miRNAs are enriched for “epileptogenic” ontology categories
To further investigate the functional role of miRNAs significantly differentially expressed and anticorrelated with their predicted mRNA targets, we examined the functional enrichment of the mRNA targets identified in epileptogenesis and chronic phases of epilepsy.
Target genes that inversely correlated with differentially expressed miRNAs during the epileptogenesis period were enriched for GO terms related to synaptic function [like “response to stimulus” (p = 0.0013), “signaling” (p = 2.68−05), “signal transduction” (p = 0.0047), and others] and immunity [like “humoral immune response” (p = 0.0009), “regulation of immune system process” (p = 0.0013), and others]. In addition, terms related to complement activation [like “complement activation” (p = 0.0002) and “complement activation, classical pathway” (p = 0.0009)] are in prominent position (Fig. 5A; Table 8). Proteins of the classical complement pathway not only play a role in the innate immune system, but have been also shown to be released from neurons, and serve as a new class of synaptic organizers (Yuzaki, 2017). These GO terms are potentially relevant to changes occurring at the level of the DG in epileptogenesis (Dudek and Sutula, 2007; Vezzani et al., 2015). In addition, at the level of cell signaling pathways, analysis of KEGG pathways enriched among the mRNA targets suggested a key role for the MAPK cascade (p = 0.354; Fig. 5B; Table 8). Notably, changes in the activation state of kinase pathways and altered kinase expression patterns have been reported in the hippocampus by previous studies (Xi et al., 2009). In the chronic period, the predicted and anticorrelated mRNA targets revealed enrichment in biological processes that have been previously implicated in chronic epilepsy (Ludewig et al., 2016; Robel and Sontheimer, 2016) such as “regulation of dendritic cell differentiation” (p = 4.8 × 10−6), “glial cell development” (p = 4.0 × 10−4), “proliferation” (p = 6.0 × 10−4), and “cell proliferation” (p = 4.0 × 10−4; Fig. 5C; Table 9). Notably, we performed a permutation test to check for false positive enrichment in GO term and KEGG pathway analysis, but this did not change any of the results.
Functional enrichment of dysregulated miRNA-mRNAs targets modules. A, Horizontal bar plots (on the left) show the GO enrichment status (top 20 terms) for 112 predicted mRNAs targets that anticorrelate with 22 miRNAs expression level in epileptogenesis (FDR < 5%, hypergeometric test). The miRNA-mRNA module is represented by a network graph (on the right) showing the connections between miRNAs based on the function of their mRNAs predictive targets revealed by the GO enrichment. B, Horizontal bar plots (on the left) show KEGG enrichment analysis for predicted mRNAs targets that anticorrelate with miRNAs expression level in epileptogenesis (FDR < 5%, hypergeometric test). miRNA-mRNA modules are represented with network plot (on the right) showing the connection between miRNAs based on the pathways in which are involved their predicted targets revealed by KEGG analysis. C, D, GO and KEGG enrichment status for 29 predicted miRNA targets that anticorrelate with five miRNAs differentially expressed in the chronic stage (FDR < 5%, hypergeometric test).
GO and KEGG enrichment of 122 predicted and anticorrelated mRNAs targets of 22 miRNAs differentially expressed in epileptogenesis
GO and KEGG enrichment results for the 29 predicted and anticorrelated mRNAs targets of five miRNAs differentially expressed in the chronic period
Target genes that directly correlated with differentially expressed miRNAs during the epileptogenesis period were enriched for GO terms related to glia proliferation [“regulation of glial cell proliferation” (p = 0.0001) and “glial cell proliferation” (p = 0.0001)]. In addition, and as in the inverse correlation analysis, terms related to complement activation [like “complement activation, classical pathway” (p = 0.0004)] were significantly enriched. In the chronic period, the positively correlated mRNA targets revealed significant (p < 0.00001) enrichment in GO terms related to receptor function like “receptor activity,” “signaling receptor activity,” “G protein-coupled receptor activity,” “signal transducer activity,” “molecular transducer activity,” and “transmembrane signaling receptor activity.”
To infer the functional relationships between miRNAs identified as differentially expressed in the meta-analysis, we created a network of miRNAs based on their predicted anticorrelated target pathways (Fig. 5). These results highlight that several distinct miRNAs may contribute to the regulation of functionally related processes and pathways, and so prioritizing individual miRNAs as potential therapeutic targets will require downstream experimental analysis.
Discussion
Main findings
The present meta-analysis provides a miRNA differential expression signature in the DG of rats during epileptogenesis and in the chronic phase of epilepsy. We identified 26 miRNAs significantly differentially expressed during epileptogenesis, and five miRNAs significantly differentially expressed in the chronic phase of epilepsy. We also identified 11 miRNAs in epileptogenesis and two in chronic epilepsy that were identified as significantly differentially expressed by the meta-analysis but not in any of the individual studies. Further, we explored the negative correlation between the significantly differentially expressed miRNAs and their predicted mRNA targets in the same models of epilepsy. We identified 122 predicted mRNAs targets with an anticorrelated expression relationship to 22 of the 26 miRNAs significantly differentially expressed in epileptogenesis. Below, we discuss these findings and their possible implications in the development and maintenance of epilepsy. Together, we also discuss the intrinsic limitations of this study that must be taken into account.
Epileptogenesis
Functional annotations of the target genes of miRNAs significantly differentially expressed during the latent interval between brain injury and the development of spontaneous seizures (epilepsy) support a relationship between dysregulated miRNAs and molecular and cellular reorganizations that are known to occur during epileptogenesis. First, the GO enrichment analysis identifies many terms that suggest a role for modulation of synaptic transmission during epileptogenesis. This is not surprising, given the critical role of the DG in the temporal lobe seizure network (Krook-Magnuson et al., 2015) and previous experimental evidence for changes in synaptic efficacy and connections during epileptogenesis (Dudek and Sutula, 2007). Another set of terms broadly refers to immunity and inflammation, events that are deeply associated with epileptogenesis (Vezzani et al., 2015). This is also supported by the identification of individual differentially expressed miRNAs (e.g., miR-146a-5p; Aronica et al., 2010) and mRNA targets (e.g., CD74 and C1r; Teo and Wong, 2010; Zeis et al., 2016) involved in inflammation.
Worthy of note is the enrichment for genes in the MAPK signaling pathway. Enrichment within the MAPK cascade has been reported during latency in the pilocarpine model in a hippocampal RNA expression study based on high-throughput RNA sequencing (Hansen et al., 2014). In particular, we found robust upregulation of two MAPKs, Map3k14, also called NIK, and Map3k4, also called MEKK4. Map3k14, a target of miR-7a-5p and miR-138-5p, mediates the neuron specific suppression of the nuclear factor κ-B (NF-kB; Mao et al, 2016) that is upregulated in epilepsy patients (Teocchi et al., 2013) and in an experimental model of traumatic brain injury (Lipponen et al., 2016). NF-kB has been linked to traumatic brain injury relevant outcomes, including epileptogenesis and tissue repair, the hypothesis being that it plays an antiepileptogenic role (Lipponen et al., 2016). Thus, Map3k14 activation may favor epileptogenesis and damage. The transcript of Map3k4, the other upregulated MAPK, is a target of four significantly downregulated miRNAs (namely, miR-92b-3p, miR-153-3p, miR-101a-3p, and miR-3573-3p). This enzyme activates the p38 and JNK pathways that are known to contribute to the apoptotic and inflammatory responses after kainate injection in mice (Yang et al., 1997; Jeon et al., 2000).
Activation of each of these suggested proepileptogenic kinase pathways might be counterbalanced by upregulation of inhibitors such as phosphatases. In our analysis, the protein tyrosine phosphatase, nonreceptor type 5 (Ptpn5), also called STEP, was found upregulated, and anticorrelated with the downregulated miR-150-5p and miR-383-5p. Contrary to Map3k4, Ptpn5 has been shown to inhibit p38 by selectively dephosphorylating its activation loop tyrosines and by sequestering it in the cytosol (Francis et al., 2014). Ptpn5 can also target the glutamate receptor subunits GluN2b and GluA2 leading to receptor internalization and decreased synaptic efficiency (Snyder et al., 2005; Xu et al., 2009). In addition, Ptpn5 inhibits the ERK2 pathway. Whereas p38 downstream molecules lead to the activation of neuroinflammation and apoptotic processes, the ERK2 cascade triggers neuronal differentiation and survival through the activation of the antiapoptotic gene bcl-2 (Cruz and Cruz, 2007).
In addition to above, miRNAs such as miR-101a-3p, miR-551b-3p, and miR-139-5p may act together to modulate the expression of the predicted target syn2, a gene that is mutated in epileptic patients (Cavalleri et al., 2007). Synapsin 2 is a member of the synapsin family composed by synaptic vesicle phosphoproteins that modulate synaptic transmission and plasticity. Notably, Syn2-knock-out mice show a decreased vesicle density at inhibitory synapses of DG GCs and are prone to epileptic seizures (Medrihan et al., 2013). Here, we found an inverse correlation of three downregulated miRNAs (miR-101a-3p, miR-551b-3p, and miR-139-5p) with the Syn2 mRNA, which levels are slightly increased suggesting that Syn2 phophoproteins, at this stage (i.e., epileptogenesis), are still able to control neuronal transmission at the DG synapses. Notably, Syn2 interacts with presynaptic Ca2+ channels to promote GABA asynchronous release (Medrihan et al., 2013), maintaining the tonic inhibition of excitatory neurons and contrasting the aberrant network synchronization that lead to seizures development in the chronic phase. These findings suggest an antiepileptogenic role of this inverse-correlation. The dentate cells, in this case, may slow down the miRNA levels to contrast the upcoming epileptogenic process.
High levels of miR-212-5p may favor the epileptogenic process through reduced expression of Gabrd. We found that the expression of the subunit δ of GABAA receptors was decreased in DG and inverse-correlated with the upregulated miR-212-5p. δ-Subunit-containing receptors are found in extrasynaptic and perisynaptic locations in hippocampal DG GCs (Wei et al., 2003). Because of their high affinity for GABA, they mediate tonic GABAA inhibition (Stell et al., 2003). Therefore, a decrease in GABAA receptor δ-subunits may impair tonic GABA inhibition, contributing to GCs hyperexcitation and seizures onset.
Chronic epilepsy
Our exploration of the chronic stage of epilepsy was more limited than that for epileptogenesis as anticorrelations (inferred via statistically significant gene and miRNA fold changes in response to the disease) with miRNA targets could be evaluated only for the amygdala stimulation model. This analysis highlighted the downregulation of five miRNAs and the upregulation of several mRNAs targets, and identified genes enriched in GO terms related to glial cells and dendritic cells (DCs). Whereas the proliferation of glia cells and their contribution to neuroinflammation and hyperexcitability in chronic epilepsy are well recognized (Devinsky et al., 2013; Robel and Sontheimer, 2016), the role of DCs in the context of epilepsy remains elusive. It can be hypothesized that DCs might be involved in epilepsy by maintaining a chronic inflammatory response (Ludewig et al., 2016).
Among the list of mRNAs predicted targets anticorrelated with significantly differentially expressed miRNAs in the chronic stage of epilepsy, of note is the glutamate ionotropic receptor δ type 2 (Grid2), which is anticorrelated with miR-130a-3p, miR-148b-3p, and miR-551b-3p. The involvement of ionotropic glutamate receptors in epileptic hyperexcitability is well established, but not much is known specifically on glutamate ionotropic type δ receptors or the regulation of this process. Further studies are needed to establish a role of these receptors in hippocampus and more specifically in epilepsy.
Limitations
The purpose (and, in our view, the strength) of this work was to maximize information from underpowered individual studies, increasing power and allowing the identification of a set of miRNAs that, being significantly and similarly dys-regulated in multiple experimental models, may be related to the disease rather than specific to a particular model. The study, however, also has limitations that should be taken into account.
A technical limitation is that the comparison of datasets in which tissue was obtained through different methods and that used different microarray platforms may have led to some miRNAs being detected in one experimental model and not in another, due to technical differences related to the assay system. Therefore, we cannot exclude the possibility that additional miRNAs were significantly dys-regulated.
Other limitations refer to biological aspects. First, miRNAs are only one mechanism of regulation of gene expression. Other changes may occur depending on other epigenetic mechanisms (histone modifications, DNA methylations) or changes in transcription factors. Second, this analysis has been conducted on one specific hippocampal subarea, the DG, enriched in a specific cell population, the GCs. Other brain areas and cell populations may be equally or even more important in epileptogenesis. Third, due to limitation in the availability of datasets, we analyzed data from a single time point in all models, but epileptogenesis may develop differently in different models. Finally, given the lack of miRNA and mRNA datasets in epileptic patients matched with valid controls, we could not verify whether the miRNA-mRNA interactions identified in rats may be relevant for the human disease.
In addition, the comparison with mRNA datasets should be viewed as a secondary outcome of the study and considered with caution. In this study, this comparison is primarily based on the assumption that miRNAs decrease levels of their mRNA targets (Guo et al., 2010) and, therefore, that target mRNAs will undergo changes in anticorrelation with those of miRNAs. Although this is the best characterized mechanism of miRNA action, it is becoming evident that miRNAs also have specific nuclear functions, including transcriptional control of gene expression and regulation of alternative splicing (Catalanotto et al., 2016), which may not lead to anticorrelation between miRNA and mRNA levels. Therefore, we also performed a further analysis of direct correlation between miRNAs and mRNAs. These analyses now require verification and further studies to establish the exact patterns of interaction between miRNAs and mRNAs in the epileptic tissue and their functional impact on epileptogenesis and maintenance of an epileptic condition.
Conclusions
The present meta-analysis identified many significantly differentially expressed miRNAs in epileptogenesis and chronic epilepsy, several of which were not uncovered in the individual studies, highlighting the additional information that can be gained by meta-analysis. Our results also highlight the added value of meta-analysis of existing data and so avoid unnecessary animal experimentation to generate new hypotheses on miRNAs involved in epileptogenesis and chronic epilepsy. Our results highlight a possible key role for a few miRNAs that are worthy of further investigation. As it may be expected, however, the number and heterogeneity of mRNAs identified by this meta-analysis suggest that therapies focused on a single miRNA target may be not sufficient to reverse or ameliorate the epileptogenic process.
Acknowledgments
Acknowledgements: We thank Dr. Ray Dingledine (Emory University, Atlanta, GA) for helpful discussions and for critically reading this manuscript.
Footnotes
The authors declare no competing financial interests.
Authors contributions: P.K.S., P.R., M.R.J., and M.S. performed research; E.P., K.L., J.A.G., E.A., and A.P. designed research; P.K.S., P.R., M.R.J., and M.S. wrote the paper.
This work was supported by funding from the European Union’s Seventh Framework Programme (FP7/2007-2013) under Grant Agreement 602102 (EPITARGET; to M.R.J., M.S., E.P., K.L., E.A., A.P.), the Imperial College National Institute for Health Research Biomedical Research Centre Scheme (M.R.J., E.P.), and the Polish Ministry of Science and Education Grant W19/7.PR/2014 (to K.L.).
↵† M.R.J. and M.S. share the position of senior author and Authors who analyzed data are P.K.S., P.R., M.R.J., and M.S.
This is an open-access article distributed under the terms of the Creative Commons Attribution 4.0 International license, which permits unrestricted use, distribution and reproduction in any medium provided that the original work is properly attributed.
References
In this issue
