Skip to main content
Download PDF

Neuronal methylome reveals CREB-associated neuro-axonal impairment in multiple sclerosis

Clinical Epigenetics volume 11, Article number: 86 (2019) Cite this article

Abstract

Background

Due to limited access to brain tissue, the precise mechanisms underlying neuro-axonal dysfunction in neurological disorders such as multiple sclerosis (MS) are largely unknown. In that context, profiling DNA methylation, which is a stable and cell type-specific regulatory epigenetic mark of genome activity, offers a unique opportunity to characterize the molecular mechanisms underpinning brain pathology in situ. We examined DNA methylation patterns of neuronal nuclei isolated from post-mortem brain tissue to infer processes that occur in neurons of MS patients.

Results

We isolated subcortical neuronal nuclei from post-mortem white matter tissue of MS patients and non-neurological controls using flow cytometry. We examined bulk DNA methylation changes (total n = 29) and further disentangled true DNA methylation (5mC) from neuron-specific DNA hydroxymethylation (5hmC) (n = 17), using Illumina Infinium 450K arrays. We performed neuronal sub-type deconvolution using glutamate and GABA methylation profiles to further reduce neuronal sample heterogeneity. In total, we identified 2811 and 1534 significant (genome-wide adjusted P value < 0.05) differentially methylated and hydroxymethylated positions between MS patients and controls. We found striking hypo-5mC and hyper-5hmC changes occurring mainly within gene bodies, which correlated with reduced transcriptional activity, assessed using published RNAseq data from bulk brain tissue of MS patients and controls. Pathway analyses of the two cohorts implicated dysregulation of genes involved in axonal guidance and synaptic plasticity, with meta-analysis confirming CREB signalling as the most highly enriched pathway underlying these processes. We functionally investigated DNA methylation changes of CREB signalling-related genes by immunohistofluoresence of phosphorylated CREB in neurons from brain sections of a subcohort of MS patients and controls (n = 15). Notably, DNA methylation changes associated with a reduction of CREB activity in white matter neurons of MS patients compared to controls.

Conclusions

Our data demonstrate that investigating 5mC and 5hmC modifications separately allows the discovery of a substantial fraction of changes occurring in neurons, which can escape traditional bisulfite-based DNA methylation analysis. Collectively, our findings indicate that neurons of MS patients acquire sustained hypo-5mC and hyper-5hmC, which may impair CREB-mediated neuro-axonal integrity, in turn relating to clinical symptoms.

Background

Multiple sclerosis (MS), a leading cause of neurological disability in young adults, is a chronic inflammatory and neurodegenerative disease of the central nervous system (CNS) characterized by autoimmune destruction of myelin and subsequent neuronal death [1]. Although the cause of MS remains largely unknown, accumulating data support the notion of MS being a complex disease influenced by genetic [2] (primarily affecting immune genes) and environmental factors [3]. While major advances have been made in understanding immune dysfunction in early phases of MS, leading to increasingly more effective disease modulatory treatments, the pathological basis of late-stage brain pathology is still largely unknown. The disease pathology is characterized by episodic inflammatory demyelination and axonal injury, which translates into clinical symptoms such as sensory, motor, visual and coordination deficits, as well as complex cognitive disturbances [4]. Apart from focal lesions, neuroimaging [5] and histopathological investigations have revealed association between seemingly unaffected normal-appearing white matter (NAWM) integrity and cognitive impairment [6,7,8]. Evidence suggests that the degree of axonal degeneration rather than demyelination in corticospinal tracts is the major determinant of clinical motor disability [9]. Neuroradiological measures, such as global and regional brain atrophy, are also correlated to disability status, underscoring the neurodegenerative aspect of MS [10, 11]. Altogether, the MS paradigm proposes that exhaustion of neuro-axonal reserve capacity and compensatory mechanisms caused by energy failure, glutamate excitotoxicity and ionic imbalance, among others, will ultimately result in permanent neurological sequelae [12]. However, due to limited accessibility of the CNS, mechanisms underpinning neuronal dysfunction in MS patients remain largely unresolved, further hindering clinical translation for the care of progressive patients.

One way to tackle this challenge is to exploit the remarkable properties of epigenetic modifications, such as DNA methylation, which lies at the interface between internal (genes) and external (microenvironment) cues. DNA methylation is a stable epigenetic mark that can impact gene regulation and/or reflect genome activity and can be accurately quantified at a single-nucleotide resolution in a genome-wide manner post-mortem. As such, DNA methylation offers the unique possibility to interrogate the genome activity underlying cellular dysfunction in conventionally inaccessible organs from autopsy tissue. DNA methylation is acquired and maintained by the action of DNA methyltransferases (DNMTs) that catalyse the covalent addition of a methyl group to cytosine (5mC). Inversely, active demethylation is catalysed by the Ten-eleven translocation (TET) family of enzymes. Importantly, the brain methylome is highly complex due to tissue- and cell type-specific variation between grey matter (GM) and white matter (WM), functionally distinct regions of the brain [13, 14] and cell types [15]. Moreover, the brain and more specifically neurons are enriched in regulatory non-canonical DNA modifications such as CpG hydroxymethylation (5hmC) [15]. A growing body of evidence suggests that in addition to its crucial role in brain development and function [15], DNA methylation could be implicated in MS brain pathology as well [16]. Nevertheless, the majority of studies of DNA methylation in healthy and diseased brain have been performed on mixed cell populations (bulk tissue) or on total DNA methylation modifications (bulk DNA methylation reflecting 5mC + 5hmC), potentially hindering the capture of complex epigenome architecture in neurons in situ.

Despite clinico-pathological evidence of accumulating neuro-axonal impairment occurring over the MS disease course, the precise mechanisms underpinning neuronal dysfunction in MS are still largely unknown. In the present study, we explored the potential of DNA methylation as a readout of neuronal genome activity in post-mortem brain tissue of MS patients and non-neurological controls. Using a strategy combining isolation of neuronal nuclei and array-based 5mC and 5hmC technology, we here present an analysis of neuronal methylome in MS patients, followed by comparison of observed changes with functional outcomes.

Results

Investigation of bulk DNA methylation changes in neurons from MS patients

We first evaluated potential sources of heterogeneity in brain tissue that could impact DNA methylation variability. To that end, we conducted genome-wide DNA methylation analyses using Illumina Infinium 450K Human Methylation Beadchip arrays (450K) on bisulfite (BS)-treated genomic DNA in NeuN+-sorted neuronal nuclei (> 99% purity) as well as the corresponding unsorted bulk tissue from the same fresh-frozen post-mortem tissue blocks (Additional file 1: Figure S1). The analysed tissue comprised various lesion phenotypes as well as NAWM originating from the same individual (Table 1, Additional file 2). We tested the effect of lesion phenotype, degree of myelination (GM/WM) and cell heterogeneity on DNA methylation values, estimated as β-values for each CpG site. As expected, tissue composition (neurons vs. bulk) had the largest influence on DNA methylation profiles (Additional file 1: Figure S1). Analyses further suggested DNA methylation differences between lesion phenotypes (Additional file 1: Figure S1). To account for neuronal subtypes, we conducted cell type deconvolution stratifying for glutamatergic (GLU) and GABAergic neurons based on previously reported differentially methylated CpG sites between NeuN+/SOX6− (GLU) and NeuN+/SOX6+ (GABA) neuronal nuclei [17] and found that neuronal subtype proportion indeed had an impact on DNA methylation profiles (see the “Methods” section and Additional file 1: Figure S2).

Table 1 Description of pilot samples and study cohorts
Full size table

We next aimed to investigate DNA methylation changes underpinning neuronal vulnerability in MS. Considering the extensive heterogeneity of GM methylome [14], likely due to pervasive phenotypic diversity of GM excitatory neurons [18], and the substantial evidence of WM abnormalities in MS [6,7,8], we focused on sorted neurons from subcortical WM. We performed DNA methylation analysis on BS-treated DNA isolated from WM-neurons in a first cohort of MS patients and non-neurological controls (cohort 1, n = 12, Table 1, Additional file 2) using 450K. After correction for confounders such as sex, age, lesion phenotype, brain region (according to antero-posterior axis) and GLU proportions (the most prominent neuronal subtype within our samples), we detected 13 significant differentially methylated positions (DMPs) between MS and controls (adjusted P value < 0.05, Table 2). Gene ontology (GO) analyses of genes associated with the candidate DMPs (unadjusted P value < 0.001, Additional file 3) suggested modest enrichment, which did not reach significance after correction for multiple testing, of nervous processes, intracellular signalling pathways and immune/xenobiotic processes (Additional file 4). Notably, some candidate DMP-associated genes map to previously identified differentially methylated regions (e.g. CREB5, DLL1, PRKG1, PRKCZ and DYRK1B loci) or dysregulated genes (e.g. SLIT3, NOTCH2, HLA-DR/-DP, GSTT1 or KCNQ1 genes) in bulk MS NAWM tissue compared to control WM tissue [19] (Additional file 3).

Table 2 BS-DMPs (adjusted P value < 0.05) associated with multiple sclerosis (cohort 1)
Full size table

Separation of true 5mC from 5hmC enables identification of predominant gene body hypo-methylation in MS neurons

In addition to the limited sample size and differences in lesion phenotypes, we hypothesized that heterogeneity in DNA modifications could be an additional contributory factor for the lack of functionally significant changes in the sample set. Indeed, 5hmC is a highly abundant stable modification in the CNS, especially in actively transcribed neuronal genes [20, 21], thus antagonizing 5mC action. Given the inability of BS-based technology to differentiate 5mC from 5hmC, it is highly plausible that specific signals might be lost or misinterpreted depending on the overlapping pattern of 5mC and 5hmC. We further tested this in an additional cohort, focusing on low-grade inflammation WM samples to further reduce heterogeneity.

To examine 5mC and 5hmC changes in MS neurons, we performed traditional BS as well as oxidative BS (oxBS) treatments in parallel, followed by 450K arrays in an additional cohort of MS patients and non-neurological controls (cohort 2, n = 17, Table 1, Additional file 2). β-values generated by BS conversion represent the total methylation signal (5mC + 5hmC), while β-values derived from oxBS conversion reflect only 5mC levels (“true” methylation). Unsupervised hierarchical clustering and MDS showed a clear separation of BS and oxBS samples (Additional file 1: Figure S3). As expected, the corresponding density profiles showed a bimodal distribution of BS and oxBS β-values, with oxBS peaking at lower β-values (Fig. 1a). This indicates that 5hmC, which is quantified by subtracting corresponding oxBS and BS β-value (ΔβBS-oxBS), accounts for a substantial fraction of the total methylation levels detected in BS samples, and is prominent in sorted neuronal nuclei. In accordance with previous reports [22], 5hmC β-values mostly ranged between 0 and 0.5 and peaked at β-values ~ 0.25, while a minor fraction of probes displayed slightly negative values (Additional file 1: Figure S4). We therefore applied a stringent stepwise probe-filtering strategy, described in detail in the “Methods” section and summarized in Additional file 1: Figure S5. We could confirm the reliability of our strategy for estimating 5hmC in comparison with the recently published method referred to as OxyBS, which relies on maximum likelihood estimation (oxy-MLE) [23], as 5hmC β-values generated by the two methods strongly correlated (Additional file 1: Figure S6).

Fig. 1
figure 1

Predominant DNA hypo-methylation and hyper-hydroxymethylation in neurons from multiple sclerosis patients. a Density plots of oxBS (purple) and BS (green) β-values after subset-quantile within array normalization (SWAN)-based type I/II normalization used for downstream analysis. Venn diagram (b) and Circos plot (c) illustrating number of significant sites and differences (Δβ) in true DNA methylation (5mC, purple), hydroxymethylation (5hmC, orange) and bulk methylation (BS, green) between patients (n = 10) and non-neurological controls (n = 7) (adjusted P value < 0.05). Red and blue circular ideograms represent hyper- and hypo-methylated sites, respectively. The outer track is an hg19 ideogram illustrating chromosome and cytoband information

Full size image

Of the 419,858 and 272,883 probes used for the subsequent 5mC and 5hmC analyses, respectively, we identified 2811 significant “true” DMPs and 1534 differentially hydroxymethylated positions (DhMPs) between MS and non-neurological controls after correction for confounders (adj. P value < 0.05, Fig. 1b, Additional file 5). It is worth noting that 19 DMPs and 19 DhMPS were classified as non-CpG methylation sites. However, the majority of these, 12 and 16, respectively, were located within non-annotated, intergenic regions, therefore not pursued further. The top DMPs and DhMPs (adj. P value < 0.01, |Δβ| > 0.05) are listed in Table 3. Interestingly, we found a striking predominant hypo-methylation (87%, 2442/2811) and hyper-hydroxymethylation (74%, 1137/1534) in MS patients compared to controls, throughout the genome (Fig. 1c, Additional file 5). Given the high correlation of DNA methylation levels between contiguous CpGs forming co-methylated regions, we sought to identify changes clustering at multiple neighbouring CpG sites, as differentially methylated (DMRs) and hydroxymethylated (DhMRs) regions (Additional file 6). We identified 472 DMRs (5mC-Δβ between − 0.27 and 0.13) and 309 DhMRs (5hmC-Δβ between − 0.10 and 0.22), the large majority of them (87% of DMRs and 88% of DhMRs) encompassing at least one DMP or DhMP, respectively. Consistent with single CpG analyses, most of the DMRs (88%, 416/472) were hypomethylated, contrasting with the predominant (78%, 243/309) hyper-hydroxymethylation of DhMRs (Additional file 6). We could confirm changes at OBSCN locus (chr1: 228503693–228503882) using a BS-free restriction enzyme-based method that distinguishes 5mC from 5hmC (Additional file 1: Figure S7).

Table 3 Top 5mC-DMPs, 5hmC-DhMPs and BS-DMPs (adjusted P value < 0.05) associated with multiple sclerosis (cohort 2).
Full size table

Analyses of BS signals generated 1434 differentially methylated DMPs (BS-DMPs, adjusted P value < 0.05) and 193 DMRs (BS-Δβ between − 0.23 and 0.21) (Additional files 5 and 6). Notably, only a very small fraction of BS-DMPs overlapped with true DMPs (42/1434) and DhMPs (67/1434) (Fig. 1b), implying that substantial differences between MS patients and controls might be missed or diluted using BS-methodology. Indeed, overlapping DMPs-DhMPs displayed anti-correlated changes (Fig. 2a), resulting in unchanged BS-Δβ values for these sites (Fig. 2b, c). This strongly suggests that opposing 5mC and 5hmC levels cancel out differences in the conventional BS signal, and that investigating 5mC and 5hmC alterations separately allows the discovery of a substantial fraction of changes occurring in neurons. Accordingly, we found a strong positive correlation (P value < 2.2.10−16, r = 0.73) between oxBS-Δβ (adjusted P value < 0.05) and BS-Δβ (P value < 0.001) at overlapping sites, and, to a lesser extent (P value < 2.2.10−16, r = 0.33) at all sites present in the array (Additional file 1: Figure S8). Collectively, these findings suggest that many BS changes with P value < 0.001 are likely true and that the sensitivity for detection of significant changes using BS data could be significantly compromised by dilution with 5hmC signals.

Fig. 2
figure 2

True 5mC and 5hmC changes are anti-correlated and enriched in gene bodies. a Spearman’s rank correlation of Δβ-5mC (oxBS) and Δβ-5hmC at 425 overlapping DMPs/DhMPs (adjusted P value < 0.05). Regression line, P value (p) and Spearman’s rank correlation coefficient (r) are given. Scatterplots (b) and boxplots (c) of 5mC (purple), BS (green) and 5hmC (orange) Δβ-values at 425 overlapping DMPs/DhMPs (adjusted P value < 0.05). Boxplots of β-values (left panel), distribution (in comparison to background, bkg) (top right panel) and Δβ-values between multiple sclerosis cases (MS, n = 10) and non-neurological controls (Ctr, n = 7) (bottom right panel) for BS-DMPs, 5mC-DMPs and 5hmC-DhMPs (adjusted P value < 0.05) across d classical gene features (shades of blue) including promoter-like features (TSS1500, TSS200, 1stExon, 5′UTR), gene body and 3′UTR and e CGI-related features (shades of red) including CpG islands, shores, shelves and open sea. In all boxplots, the black line represents the median, the interquartile range is comprised by a box and whiskers are extended to maximum and minimum Δβ values. *p < 0.05, **p < 0.01 and ***p < 0.001 using Fisher’s exact test for enrichment and depletion analyses (full data are shown in Additional file 1: Figure S9)

Full size image

As previously mentioned, 5hmC has been shown to display a unique genomic distribution compared to 5mC and to exert antagonistic roles on regulation of transcription [15]. We stratified 5mC and 5hmC changes based on gene features annotated in 450K and found overall hypo-methylation and hyper-hydroxymethylation in MS patients compared to controls across all gene features (Fig. 2d). The largest 5mC and 5hmC changes occurred within the gene body and 3′UTR, with a smaller contribution from regions harbouring promoter-like features (TSS1500, TSS200, 5′UTRs and 1st exon) (Fig. 2d). Fisher’s exact test revealed that 5mC-DMPs are strongly enriched (P value = 1.16 × 10−39) within gene bodies while depleted (P value = 4.24 × 10−23 for TSS1500) from segments with promoter-like features (Fig. 2d, Additional file 1: Figure S9). Finally, analysis in relation to CpG islands (CGIs) indicated that CGIs displayed the lowest 5mC, 5hmC β-values and |Δβ|-values compared to shores and shelves (Fig. 2e). We found enrichment of 5hmC-DhMPs in CGIs (P value = 1.85 × 10−07), which are depleted (P value = 3.27 × 10−09) of 5mC-DMPs (Fig. 2e, Additional file 1: Figure S9). Of note, the backgrounds used for 5mC and 5hmC analyses are noticeably different.

Thus, all explored gene features generally showed evidence of negatively correlated changes in 5mC and 5hmC, with DMPs and DhMPs predominantly enriched in gene bodies and CpG islands, respectively.

5mC and 5hmC changes of MS neurons associate with genes involved in CREB signalling pathway, synaptic plasticity and axonal guidance

To gain insight into biological functions associated with changes in 5mC and 5hmC, we performed ingenuity pathway analysis (IPA) on 2448 genes associated to all significant sites, i.e. 2811 DMPs + 1534 DhMPs (adj. P value < 0.05). Top enriched pathways associated to D(h) MPs remained when relaxing the significance of CpG sites (unadjusted P value < 0.001), and reflected changes in both 5mC and 5hmC (Fig. 3a, b Additional file 4). The most significant canonical pathways were related to nervous system processes, particularly axonal guidance (Fig. 3c), cAMP response element-binding (CREB) signalling and synaptic plasticity, to other intracellular signalling pathways such as actin cytoskeleton signalling and to oxidative stress/inflammation processes (Fig. 3a, b, Additional file 4). IPA analysis of BS-DMPs showed modest enrichment for intracellular signalling pathways such as xenobiotic metabolism and integrin signalling (Additional file 4). GO findings were confirmed when focusing at region levels where DMRs + DhMRs were enriched in genes implicated in regulation of neuron projections (e.g. PLXN4A, NTN1/5, RHOT1/2, CYFIP1, DOCK1, KIF1A, PACSIN1, RIN1 genes), neuronal development (e.g. CFL1, ISL1, ADAP1, DIABLO, PTK2B, MYBL2 and TP73), as well as ion channels (e.g. KCNQ1, KCNAB2), and GABA/glutamate activity (e.g. GABBR1, GRIN2D, LRP1, ABAT) (Additional file 4). Examples of dendritic and axonal (TBCD, APC2), synaptic (ADORA2, ELFN1, SLC8A2, CACN1H) and neurogenic (DMRTA2, LOXHD1) DMR-associated genes are illustrated in Fig. 3d and e. Findings were further supported by GO analysis of the predicted target genes from differentially methylated miRNAs identified in cohort 2 (exemplified in Fig. 3d, Additional file 4), axonal guidance being the most enriched pathway (Additional file 4).

Fig. 3
figure 3

Functional association of 5mC and 5hmC changes to axonal guidance, CREB signalling and synaptic plasticity. a Top canonical pathways associated with “true” methylation (5mC) and hydroxymethylation (5hmC) between multiple sclerosis cases (MS, n = 10) and non-neurological controls (Ctr, n = 7) using Ingenuity Pathway Analysis (IPA). Results for all adjusted P value < 0.05 5mC-DMPs +5hmC-DhMPs (black) as well as 5mC-DMPs (purple) and 5hmC-DhMPs (orange) separately. Significance is represented as − log10(B-H, P value) after adjustment using Benjamini-Hochberg (B-H) correction. b Multidimensional scaling of over-represented biological function terms associated with 5mC-DMPs (purple) and 5hmC-DhMPs (orange) (adjusted P value < 0.05), separately, between MS cases (n = 10) and Ctr (n = 7), according to semantic similarities, using REVIGO. The circle size represents − log10(B-H, P value) after adjustment using Benjamini-Hochberg (B-H) correction. c Representation of the genes from axonal guidance network using STRING analysis. Grey gradient indicates the strength of data support (darker grey representing stronger evidence) and colours represent different cluster (kmeans clustering set at 5). d Plots illustrating genes associated with top 5mC-DMR and 5hmC-DhMR (left and middle panel) or with BS-DMRs and 5hmC-DhMRs (right panel) (|Δβ| > 0.05). e Plots illustrating genes associated with top 5mC-DMR, 5hmC-DhMR or BS-DMRs (Δβ > 0.05). The hg19 ideogram illustrating chromosome and cytoband information, the complete gene structure of the locus (blue track), CpG Island (CGI) feature (green track), probes and DMR (black) location are shown. Methylation (β-values) of single CpGs within and outside the DMR for individual cases and controls is depicted in red and blue, respectively, with connecting lines indicating mean methylation for each consecutive CpG

Full size image

Considering the limited sample size, we aimed to validate pathways and functions by performing a meta-analysis of the BS data generated for both cohorts. We took advantage of the strong correlation identified between true 5mC and BS-generated DNA methylation changes (Additional file 1: Figure S8). A total of 265,129 common BS-derived probes considered homogeneous between the two studies (I2 < 15%, Additional file 1: Figure S10) were examined. Meta-analysis identified 8281 positions (Benjamini-Hochberg adjusted P value < 0.05, Fig. 4a, Additional file 7), the large majority (97%) exhibiting same direction of changes in both cohorts (Fig. 4b, Additional file 1: Figure S10). Pathway investigation of the genes associated with DMPs confirmed strong enrichment of CREB signalling in neurons and associated pathways, followed by other nervous system processes (e.g. axonal guidance), intracellular signalling pathways and inflammatory/oxidative processes (Fig. 4c, Additional file 4). Accordingly, altered genes are predominantly involved in CREB signalling (Fig. 4d) and comprise several subunits of the ionotropic NMDA (GRIN1, GRIN2A-C genes), AMPA (GRIA4), delta (GRID2) and kainate (GRIK3-4 genes) and metabotropic (e.g. GRIM1, GRM3-4) glutamate receptors as well GABA receptors (e.g. GABBR1, GABRA1, GABRB3, GABRG2). Downstream pathways involve multiple signalling molecules such as PKA and PKC genes, downstream cAMP-dependent kinases (CAMK2/4 genes) and transcriptional regulators (CREBBP, CREB5, EP300), among others. Axonal guidance genes are associated with ephrin/Eph receptors (e.g. EFN3-5, EPHA1, EPHA2, EPHA4, EPHB1, EPHB6) and semaphorin/plexin (e.g. SEMA4A, SEMA6A-C, PLXNA2, PLXND1) with Slit/ROBO (SLIT1/3, ROBO1/2) complexes together with actin-cytoskeleton molecules as well as trophic factors (BDNF, NGF) and downstream TGF-β, Shh and Wnt signalling pathways (e.g. BMP4, WNT5A, GLI2). Finally, inflammatory processes include cytokine signalling (e.g. TNF, IL1RAP, IRAK3, TRAF6, MAP2K4, MAP3K1, MAPK11), oxidative pathways (e.g. NOS3, SOD3, NOSTRIN, HSPA1A/B, HSP90B1) and xenobiotic metabolism such as the aryl hydrocarbon receptor (AHR, AHRR), glutathione- and cytochrome P450-mediated detoxification enzymes (e.g. GSTM1/3/5) (Additional file 4).

Fig. 4
figure 4

Meta-analysis of cohorts and comparison of DNA methylation with transcriptional changes. a Volcano plot illustrating differences in M-values between multiple sclerosis (MS) cases and non-neurological controls (adjusted P value < 0.05) from meta-analysis of cohorts 1 and 2. b Correlation between effect sizes (ΔM) of cohorts 1 and 2 for DMPs identified in meta-analysis (adjusted P value < 0.05) revealed predominant same direction of changes in the two cohorts. Regression line, P value (p) and Spearman’s rank correlation coefficient (r) are given. c Top canonical pathways of DMPs (adjusted P value < 0.05) identified in meta-analysis. Significance is represented as − log10(B-H, p value) with the red line showing threshold of significance after adjustment using Benjamini-Hochberg (B-H) correction. d Schematic representation of CREB-signalling pathway associated with DMPs from meta-analysis, with differentially methylated genes between cases and controls depicted in blue. e Association of 5mC-DMPs (left) and 5hmC-DhMPs (right) with gene expression data (RNA-seq) in bulk NAWM vs control WM [19]. Proportions (percentage) of number of upregulated and downregulated genes are depicted in red and blue, respectively (the dotted line represents the expected proportion). P values (p) generated with the Chi-square test are given for DMPs located in gene body or promoter-associated features (TSS1500, TSS200, 5′UTR and 1stExon) (full data in Additional file 8). f Scatter plots illustrating association between 5mC-DMPs (purple, left) and 5hmC-DhMPs (orange, right) methylation changes (Δβ-value) from cohort 2 (adjusted P value < 0.05) with gene expression data (RNA-seq) in bulk NAWM vs control WM [19]. The number of samples used in cohorts 1 and 2 equals to n = 12 and n = 17, respectively

Full size image

We then investigated the putative functional impact of DNA methylation changes by examining corresponding differentially expressed (DE) genes reported in RNAseq data from bulk MS-NAWM compared to control WM [19]. Of 4669 genes associated with 8281 DMPs from the meta-analysis, 84.4% were detected in the RNAseq analysis, among which 13% (510) were differentially expressed (unadjusted P value < 0.05). Of those, 61% (312/510) were significantly downregulated in bulk MS NAWM compared to control WM [19] (Additional file 8). Interestingly, genes harbouring methylation changes within gene body showed the largest proportion of downregulated genes (72%), which is significantly more than expected by chance (p = 4.5 × 10−13, Chi-square test), as opposed to changes within region of promoter-like features (54%) (Additional file 8). A predominant downregulation was also evident for true 5mC and 5hmC (cohort 2) (Fig. 4e,f), with for example 76% (129/170) of genes affected by gene body 5mC-DMPs displaying a significant decreased expression in bulk NAWM compared to WM (p = 6.0 × 10−08, Chi-square test, Additional file 8). This suggests that gene body hypo-methylation is likely associated with decreased transcriptional activity, as previously suggested [24, 25].

Altogether, these data strongly suggest that WM-neurons from MS patients manifest 5mC and 5hmC differences that might impair neuronal homeostatic functions, possibly through transcriptional regulation of associated genes. Most of the changes converge on alteration of CREB signalling pathway, synaptic plasticity and axonal guidance.

DNA methylation changes associate with reduced CREB activity in NAWM neurons

DNA methylation analysis identified alteration of multiple genes involved in CREB response in MS neurons, some of which were dysregulated in transcriptome or methylome studies of heterogeneous MS brain tissue [19] (Additional file 7). We then asked whether DNA methylation changes in CREB-related genes could associate with alteration of CREB activation in neurons from MS patients. We addressed this by examining the activation status of CREB transcription factor using immunofluorescence (IF), reflected by nuclear phosphorylated CREB (P-CREB+) in MS and non-neurological controls (n = 15, a cohort overlapping cohort 2, Table 1, Additional file 2). Group-representative images are shown in Fig. 5a. Analysis revealed a significantly lower number of neurons with nuclear phosphorylated CREB (P-CREB+) in the NAWM of MS patients (18.3 ± 11.0%) compared to controls (72.5 ± 23.9%) (p = 3.4 × 10−02, Kruskal-Wallis with Dunn’s multiple comparisons test, Fig. 5b). In contrast, no significant change in the activation status of CREB could be detected in normal-appearing GM from MS patients compared to control GM (Fig. 5b), implying that NAWM neurons of the MS brain seem to be more susceptible to dysfunction of the CREB signalling pathway than the GM ones.

Fig. 5
figure 5

Reduction of CREB activity in normal-appearing white matter neurons. a Representative images of multiple sclerosis (MS) and non-neurological control (Ctr) brain tissue sections co-targeted with phosphorylated CREB (P-CREB, green, left panel) and neuronal marker NeuN (red, middle panel); merged with DAPI (blue, right panel) assessed by immunofluorescence (IF) (magnification 630x, scale bar 15 μm). b Quantification of IF images from MS cases (n = 7) and Ctr (n = 8). The amount of NeuN+/CREB (phospho S133)+-neurons is presented as a percentage of the total amount of all detectable NeuN+ cells. In contrast to NAWM, CREB activation status in the NAGM neurons was less affected and comparable between the groups. *p < 0.05, *p < 0.01 using Kruskal-Wallis test and Dunn’s test for multiple comparisons. NAWM, normal-appearing white matter, NAGM, normal-appearing grey matter, WM, White matter, GM, Grey Matter

Full size image

Discussion

We here utilized DNA modifications to investigate processes that occur in neurons from MS patients. To our knowledge, this is the first report of 5mC and 5hmC abnormalities in neurons, especially WM-neurons, in the context of neurological disease. We found striking hypo-methylation and hyper-hydroxymethylation at specific CpGs, some of which clustered into regions and occurring mainly within gene bodies. Our findings implicate alterations of genes involved in axonal guidance, synaptic plasticity and CREB signalling, as putative key processes that could contribute to impaired neuro-axonal integrity and inability to repair after immune insult, and perhaps also be more prone to degenerate. Importantly, we could associate DNA methylation changes with functional outcomes by reporting for the first time a reduction of CREB activity in NAWM neurons compared to control WM neurons.

Our data show opposing hypo-5mC and hyper-5hmC changes, which seem to co-localize mostly within gene bodies. Additionally, we observed that a large majority of genes acquiring gene body hypo-5mC/hyper-5hmC was significantly downregulated in bulk MS-NAWM vs. control WM [19], linking these DNA modifications with transcriptional activity. Several studies have shown that hypo-5mC in post-mitotic neurons associates with impaired survival and excitability [26,27,28]. Likewise, hyper-5hmC can lead to neuronal hypersensitivity [29] and defects in neurite outgrowth and synaptic formation[30]. Moreover, a growing body of evidence suggests that global hypo-5mC and hyper-5hmC changes result in subsequent TET-mediated activation of specific endogenous retroviral (ERV) elements [31, 32], that can ultimately result in severe postnatal neurodegeneration in rodents [31]. Importantly, strong transcriptional activity of human endogenous retroviruses, including in CNS cells, has been detected in MS patients [33]. Thus, altered 5mC and 5hmC patterns might contribute to neuronal vulnerability by impairing chromatin regulatory architecture in MS neurons, leading to dysregulation of neuron-specific genes and enhanced genome instability.

Whereas specific mechanisms affecting levels of 5mC and 5hmC in MS neurons remain to be explored, several processes such as oxidative stress, inflammation or hypoxia together with putative upstream regulators identified in our study (e.g. TGF-β1, POU5F1, CREB1, Additional file 4), may impact locus-specific 5mC and 5hmC profiles. Moreover, an emerging body of evidence suggests a putative role of dysregulated epigenetic enzymes such as DNMT, TET and MBD genes in MS, as observed in blood cells [34] and brain [35] from MS patients compared to controls. This is further supported by a recent study showing reduced methionine levels in plasma of MS patients and the regulatory effect of peripheral methionine on DNMT3A in the mouse brain [36]. Consistent with this, we found that several key DNA methylation enzymes, namely DNMT1, DNMT3A, DNMT3B or MBD3, among others, displayed significant 5mC and/or 5hmC (cohort 2) and BS (meta-analysis) DNA methylation changes in MS neurons compared to controls. Finally, since DNA methylation levels can be affected by variations in DNA sequence, it may be speculated whether MS risk alleles such as that on chromosome 4 encompassing the TET2 gene [37], might contribute to the observed changes in 5mC and 5hmC levels. In line with this, TET2 promoter exhibited reduced 5mC levels in MS neurons (cohort 2). Thus, the combined influence of neuroinflammatory processes occurring as a result of MS as well as putative upstream MS disease risk factors may dually impact on 5mC and 5hmC profiles in MS neurons.

While axonal injury is believed to be a major contributor to disability such as cognitive decline in late MS [6,7,8], comprehensive characterization of the underlying mechanisms is still lacking. Moreover, accumulating evidence points to a role of glutamate excitotoxicity [38, 39] and synaptopathy [40, 41] in neuro-axonal dysfunction and degeneration in MS. Importantly, glutamate levels measured in NAWM of MS patients are predictive of neuro-axonal integrity, brain atrophy and cognitive impairment [42, 43]. In our study, neurons from MS patients displayed epigenetic alterations affecting several genes of the glutamate/GABA signalling, including multiple subunits of glutamate and GABA receptors. Additionally, our data suggest that compromised axonal architecture likely results from changes in interconnected cellular networks ranging from surface complexes (semaphorin/plexin, Ephrin, Slit/ROBO) to genes from the cytoskeleton and Shh/Wnt-signalling pathways. Comparison of samples from the same individual containing different lesion phenotypes (pilot samples, Additional file 1: Figure S1) suggests lesion-associated changes in genes implicated in neuronal projections (e.g. EPHA10, NIN, ALCAM) and synaptic processes (e.g. GABRA5, PRKG1, DLGAP3/SAPAP3) as well. Thus, neurons of MS patients exhibit epigenetic alterations that likely reflect neuro-axonal damage and/or failure in compensatory mechanisms.

Interestingly, DNA methylation changes converged on genes in the CREB signalling pathway, namely PKA, PKCs, CAMKs kinases and CREB genes. Transcription factors of the CREB family play pivotal roles in axonal regeneration, plasticity, cell survival, oxidative stress and neuroprotection [44, 45]. Prior studies have reported dysregulated CREB signalling in lesions from bulk/mixed MS brain tissue [40, 46]. Importantly, we could associate DNA methylation differences with reduction of CREB phosphorylation in NAWM neurons from MS patients compared to controls, which is to our knowledge the first report of altered CREB activity in NAWM in MS. Undeniably, one cannot exclude additional mechanisms affecting CREB activity independently of DNA methylation changes. Nevertheless, this is highly relevant considering the CREB-mediated neuroprotective effect exerted by the MS drug fingolimod (Gilenya™; Novartis, Basel, Switzerland) in vitro and in vivo in a model of neurodegenerative motor disease [47]. Interestingly, CREB activation differed between WM and GM neurons, likely reflecting differences in neuronal subtypes and/or distinct susceptibility towards demyelination pathology [48], such as hypoxia-like metabolic injury [49]. Altogether, these findings strongly suggest that epigenetic and functional alteration of CREB signalling can occur without or prior to focal tissue damage.

Undoubtedly, epigenetic marks are critical features of cellular differentiation and phenotype but multiple confounders could bias proper interpretation of DNA methylation data. Therefore, we aimed to minimize tissue- and DNA modification-driven sources of heterogeneity by studying neuronal nuclei isolated from WM, correcting for neuronal subtype proportions and separating true 5mC from 5hmC. To our knowledge, this strategy has not been previously reported in a neuropathological context. Our results suggest that BS-based DNA methylation studies might preferentially capture strong difference in one epigenetic mark and that examining 5mC and 5hmC separately could greatly aid in detecting novel changes that might have escaped detection by conventional BS-based DNA methylation studies. Nevertheless, in our study, BS-DMP changes appear strongly correlated with 5mC differences. Moreover, some of our findings could replicate previously reported DNA methylation results and could further associate with differential gene expression in bulk NAWM compared to control WM [19]. The aberrations in WM-neurons of MS patients identified here are interesting in light of neuronal circuit disconnection [50] and subsequent disability in progressive MS. Additional cell type-specific studies in larger cohorts can shed further light on mechanisms underlying neuronal dysfunction and associated disability. In the long term, a better understanding of functional implications of epigenetic changes in CNS cells in MS may lead to novel therapeutic strategies aiming at restoring neuronal integrity and ameliorating cognitive deficit in progressive MS [51, 52].

Conclusions

Our study demonstrates functionally relevant DNA methylation alterations in WM-neurons from MS patients in comparison to non-neurological controls, thereby providing new insights into mechanisms underlying neuro-axonal pathology of the WM and clinical symptoms in MS patients. Furthermore, our findings open new perspectives for similar approaches based on deciphering true 5mC and 5hmC changes in neurons specifically in other neurological disorders.

Methods

Subjects and cohorts

Brain tissues of all the cohorts was obtained from the Multiple Sclerosis and Parkinson’s Tissue Bank (Imperial College London). Briefly, the pilot samples were used to assess the impact of sample heterogeneity on DNA methylation profiles, cohorts 1 and 2 were used to analyse neuron-specific DNA methylation changes between MS cases and controls and cohort-IF overlapping with cohort 2 was utilized for functional validation of CREB activity in brain sections. Further sample details are given in Additional file 2. The material comprises snap-frozen brain tissue blocks collected within 33 h post-mortem and divided by specialist into histopathologically characterized lesion categories: active lesion (AL), chronic active lesion (CAL), chronic inactive lesion (CL) and normal-appearing white matter (NAWM). Control subjects were selected based on a non-neurological cause of death. The samples were further annotated for cerebral location (antero-posterior axis), characteristic of the tissue (mixed, white or grey matter) using the human brain atlas sectional anatomy database (http://www.thehumanbrain.info/) and were further dissected accordingly. The anatomical localization of the samples across brain regions is illustrated in Additional file 1: Figure S11. Of note, different brain samples from same individuals were used between cohort 1 and 2, and, four identical specimens derived from eight individuals were used for DNA methylation study (cohort 2) and functional validation. No formal sample size calculation was conducted, samples reaching the following inclusion criteria have been included in DNA methylation analyses: (1) all available samples with sufficient DNA amount, (2) samples that passed DNA methylation quality control and (3) cases with confirmed MS diagnosis and non-neurological controls without any signs of inflammation in the CNS. The pilot samples were selected based on the availability of different lesion phenotypes within the brain of one MS patient. Cohort 1 was selected based on available amount of DNA from WM neuronal nuclei. Cohort 2 was selected based on sufficient DNA amount from WM neuronal nuclei to perform both 5mC and 5hmC analyses and based on the low-grade inflammation of the lesion phenotype with only data from WM-neuronal nuclei being used for the comparison MS versus non-neurological control. Cohort used for functional validation was selected based on the availability and quality of the material.

Sample preparation

Fluorescence-activated cell sorting (FACS)-based neuronal nuclei isolation from dissected brain tissue was performed according to previously published protocol [53] (see representative experiment in Additional file 1: Figure S1). Briefly, following resuspension of the homogenized brain tissue in hypotonic lysis buffer (0.32 M sucrose, 5 mM CaCl2, 3 mM MgAc2, 0.1 mM EDTA, 10 mM Tris pH.8, 1 mM DTT, 0.1 % Triton), nuclei were extracted by ultracentrifugation in sucrose gradient (1.8 M sucrose, 3 mM MgAc2, 1 mM DTT, 10 mM Tris pH.8) for 2.5 hours at 4 °C. Nuclei were further labelled with Alexa Fluor 488 (Invitrogen #A11029)-conjugated anti-NeuN antibodies (1:700, Millipore #MAB377) and separated into neuronal and non-neuronal nuclei by flow cytometry (MoFloTM high-speed cell sorter). First, FSC [Par] × SSC gating was used to separate larger particles from smaller debris. Next, FSC [Par] × Trigger Pulse Width plot was used to remove aggregated nuclei such as duplicates. The fluorescence event plot showed two clear populations including the NeuN-positive fraction (representing 4–40% of the nuclei). We confirmed that positive fractions were negligible in negative controls, i.e. no antibody or unconjugated Alexa Fluor 488-labelled nuclei. Neuronal nuclei were pelleted and stored at − 80 °C until DNA isolation. Genomic DNA was isolated using QIAmp DNA micro kit (QIAGEN), resuspended in water and stored at − 80 °C.

Illumina Human Methylation 450K

We used Illumina Infinium Human Methylation 450K BeadChip (Illumina, Inc., San Diego, CA, USA; 450K) for quantitative and genome-wide DNA (hydroxy) methylation profiling. Genomic DNA was subjected to either conventional BS-treatment for the pilot samples and cohort 1 or BS and oxidative BS (oxBS)-conversion using TrueMethylTM 96 kit of CEGXTM (Cambridge Epigenetix Limited) for cohort 2. BS-DNA from the pilot samples and cohort 1 was hybridized to 450K arrays at BEA core facility (Karolinska Institutet), oxBS/BS-DNA from cohort 2 were processed at GenomeScan (GenomeScan B.V., Leiden, The Netherlands), according to manufacturer’s instructions and the BeadChip images were scanned on the iScan system. Samples were randomized ensuring that disease group, gender and age were balanced to control for potential confounding effects. Technicians performing 450K arrays were blinded to the MS disease status during the experiments. Persons performing statistical analysis were not blinded to disease status. The analysts have never altered the diagnosis of samples and no individuals were excluded because of diagnosis.

DNA methylation and hydroxymethylation analyses

Quality control

450K data (485,577 probes) were quality assessed using MethylAid [54], which examines Red(R)/Green(G) signal intensity, bisulfite conversion, specificity, staining, extension, target removal and hybridization as well as overall performance of the assay. All samples passed quality control and were subsequently processed using the Chip Analysis Methylation Pipeline (ChAMP) version 1.8.0 [55] and minfi version 1.16.0 [56] Bioconductor packages.

Probe filter

Upon loading raw IDAT files into ChAMP, probes were filtered by detection P value > 0.01, bead count < 3 in at least 5% of the samples, SNPs (minor allele frequency > 1% in the European population) and cross-reactivity as identified by Nordlund et al. [57] (pilot study) and Chen et al. [58] (cohorts 1 and 2). After filtering, the remaining probes for the pilot samples, cohort 1 and cohort 2 reached 374,756, 427,712 and 419,958, respectively. Notably, probes located on the X and Y chromosomes were also removed, as samples from both females and males were included in the study.

Between and within-array normalization

β-values of remaining probes were either between-sample quantile normalized followed by within-sample Beta-mixture quantile normalization (BMIQ)[59] as previously recommended [60] in ChAMP for the pilot study or within-sample normalized using the “Subset-quantile Within Array Normalization” (SWAN) method[61] as previously recommended for oxBS-450K data [22] in minfi for cohorts 1 and 2. Noticeably, within-sample normalization corrects for two different probe designs (type I and type II probes) included on the 450K BeadChip.

Hydroxymethylation (5hmC)

Filtering strategy and pipeline workflows are illustrated in Additional file 1: Figures S4 and S5. The champ. TrueMethyl function (ChAMP Bioconductor package version 1.8.0) [55] was applied to SWAN-normalized β-values to identify the “most variable positions” between oxBS and BS samples. The default Benjamini-Hochberg cut-off (B-H. adj. P value < 0.05) was used and filtered 110,666 sites, which predominantly encountered probes with mean hydroxymethylation levels around 0. Negative average hydroxymethylation sites (3383 probes, ~ 1%), were considered false positives and therefore also removed. P value distributions confirmed that mean 5hmC values > 0 had lower P values than mean 5hmC values < 0 (Additional file 1: Figure S4). Of the remaining 305,809 probes, 5hmC β-values were calculated by subtracting BS and oxBS β-values. Probes with > 1 negative value (28,421 probes) were filtered and since slide 4 only contained 2 neuronal nuclei DNA samples, no negative 5hmC values were allowed on slide 4 (causing 4505 probes to be removed) to allow for subsequent batch/slide correction with ComBat [62]. For comparison, we tested another method based on maximum likelihood estimation (MLE) available through the oxyBS version 1.0 Cran package [55] with default settings.

Correction for slide-effect

Slide effects were corrected using empirical Bayes methods [62] implemented in the ComBat function of the SVA Bioconductor package (version 3.18.0). Principal component analysis (PCA), combined with cofactor association testing before and after ComBat, confirmed that the slide-effect was successfully removed.

Neuronal subtypes deconvolution

DNA methylation sites previously identified to significantly differ between GABA and GLU (FDR < 0.05) were retrieved from [17] and filtered for |Δβ| > 0.7 which resulted in a total of 162 neuronal subtype-specific sites. In concordance with previous observations [17], the majority of the 162 examined neuronal subtype-specific CpG sites, were hypomethylated in GLU compared to GABA neurons (Additional file 1: Figure S2). Notably, due to probe filter (for details see above), the number of cell type specific probes was reduced to 144 for cohort 1 and 143 for cohort 2, respectively. GABA and GLU cell proportions were estimated from raw BS β-values using Houseman’s reference-based algorithm [63] as implemented in the projectCellType() function of the minfi Bioconductor package [56]. Estimated cell proportions were confirmed using the robust partial correlations (RPC) method from the EpiDISH (Epigenetic Dissection of Intra-Sample Heterogeneity) R package version 1.0.0 [64] (Additional file 1: Figure S2). While we found no significant differences in GABA and GLU neuronal proportions between MS and non-MS controls in neither cohorts (Additional file 1: Figure S2), Spearman’s rank correlation analysis revealed correlation with top PCs (Additional file 1: Figure S2), indicating that GLU proportions (the most prominent neuronal subtype in our samples) contribute to variation of DNA methylation in our samples. Sensitivity modestly increased after including GLU proportions as a co-variate in the linear model (Additional file 1: Figure S2). Heat maps and scatter plots of GABA and GLU specific probes were generated using the “ComplexHeatmap” Bioconductor package version 1.17.1 [65]. Putative differences in cell proportions between MS and non-MS controls were assessed with non-paired t tests and visualized with the ggboxplot() function of version 0.1.6 ggpubr R package. Venn diagrams were generated using the draw.pairwise.venn() function of version 1.6.20 VennDiagram R package.

Differentially methylated positions (DMPs) and regions (DMRs)

The Limma Bioconductor package version 3.26.3 [66] was used for detection of DMPs with M-values (Mi = log2(βi/(1 − βi))) as input as previously recommended [67]. The following covariates, as confirmed by PCA, were included in the model: Individual, Sex, Age, Lesion phenotype, Brain localization (according to antero-posterior axis) and GLU proportion for cohort 1 and Sex, Age, Lesion phenotype, Brain localization and GLU proportion for cohort 2. The influence of brain regionality has been investigated using association analysis, covariate regression and randomization and potential confounding effects of brain regionality have been excluded. DMRcate version 1.6.53 [68], which identifies DMRs based on kernel smoothing, was applied with default settings (λ = 1000, C = 2).

Meta-analysis

Meta-analysis of cohort 1 and 2 BS-derived 414,306 common probes was conducted using the inverse variance based method of the METAL tool, which weights the effect size for each study by their standard error [69]. Effect size estimates were retrieved as logFC (M-value based) and standard error estimates as sqrt (fit$s2.post) × fit$stdev.unscaled from Limma outputs, respectively. Simultaneous heterogeneity testing allowed for subsequent filtering of 149,177 heterogeneous probes based on an I2 threshold of 15% as previously suggested [70]. P values were adjusted for multiple comparisons using the Benjamini & Hochberg (B-H) method [71].

Gene annotation

Classical 450K annotations (TSS200, TSS1500, 1stExon, 5′UTR, Gene body, 3′UTR, CGI, Shelf, Shore and Open Sea) were derived from the “IlluminaHumanMethylation450kanno.ilmn12.hg19” version 0.2.1 and ChAMP version 1.8.0 Bioconductor packages. CpG islands (CGIs) were defined as GC content > 50%, observed/expected CpG ratio > 60%, and > 200 bp while CGI shore and shelves represent within and outside a 2-kb flanking region surrounding a CGI, respectively. Fisher’s exact test integrated in R (version 3.4.3) was used to estimate enrichment (alternative = “greater”) or depletion (alternative = “less”) of features of interest.

Locus-specific validation

“True” DNA methylation analysis of a CCGG motive located at chr1: 228503770 (included in OBSCN locus DMR chr1: 228503693–228503882) was carried out using methylation- and glucosylation-sensitive digestions of genomic DNA. Briefly, 100 ng of genomic DNA was mixed with UDP-glucose and 4 units of 5hmC-glucosytransferase (Quest 5-hmC Detection kit, ZymoResearch), allowing 5hmC to be glucosylated (glycosyl-5hmC). Mock glucosylation consists of all of the above with the exception of glucosyltransferase. Mock- and glucosylated-DNA were subsequently incubated with either HpaII, MspI (EpiJET, ThermoFisher Scientific) or in absence of enzyme (undigested). Unmethylated genomic DNA (EpiTect Control DNA, Qiagen) was used as a negative control. Methyl-sensitive (MSRE) HpaII enzyme cuts only unmodified CCGG motif whereas glucosyl-sensitive (GSRE) MspI enzyme cleaves all modified CCGG except glucosyl-5hmC. Detection was performed by qPCR on a BioRad CFX384 Real-Time Detection System with SYBR green fluorophore and the following primers: OBSCN_F: GCTGCTGCTCAAAAACTTGC and OBSCN_R: AATGCGGACGTCACCATATC. The percentage of total 5mC + 5hmC was quantified using the 2−ΔCt method comparing Mock- and HpaII-digested DNA with Mock- and undigested-DNA (with MspI as a positive control). Percentage of 5hmC was quantified by applying 2−ΔCt to Glucosylated- and MspI-digested DNA compared to Glucosylated- and undigested-DNA. True 5mC was expressed as % (5mC + 5hmC) − % 5hmC.

MiRNA target genes prediction

Differentially methylated loci mapping to miRNAs (UCSC Refseq annotations) from cohort 2 (including 5mC-, 5hmC- and BS-DMPs and − DMRs) were used as input for target prediction using mirDIP (version 4.1.8.2). Of note, corresponding mature miRNAs include: hsa-miR-339-5p, hsa-miR-1182, hsa-miR-1228-3p, hsa-miR-136-5p, hsa-miR-154-5p, hsa-miR-1909-3p, hsa-miR-202-3p, hsa-miR-377-3p, hsa-miR-431-5p, hsa-miR-432-5p, hsa-miR-496, hsa-miR-518, hsa-miR-548, hsa-miR-661, hsa-miR-1226-3p, hsa-miR-300, hsa-miR-802, hsa-miR-1251-5p, hsa-miR-19b-3p, hsa-miR-17-5p, hsa-miR-18a-5p, hsa-miR-19a-3p, hsa-miR-20a-5p, and hsa-miR-92a-3p. Two miRNA loci (MIR518A2 and MIR548) could not be found in the database. mirDIP integrates predictions from 30 independent resources and offers an integrative score which has been shown to provide more accuracy compared to confidence scores from the individual resources [72]. Among the predicted target genes, only the ones showing integrative scores > 0.7 were further examined using gene ontology analysis.

Gene ontology analyses

Genomic locations of 450K probes were obtained from the ChAMP Bioconductor package version 1.8.0 [55]. Gene ontology (GO) analysis was performed using ingenuity pathway analysis (IPA) (Qiagen), applying unbiased parameters for all criteria including tissues selection. The data were analysed focusing on canonical pathways. Right-tailed Fisher’s exact test was used to calculate P values. Enriched GO terms with adjusted P values < 0.05 (Benjamini-Hochberg, B-H) were considered statistically significant. Of note, analyses of BS-generated DMPs from cohorts 1 and 2 did not result in significant pathways after B-H adjustment, and, in this case, enriched pathways with P values < 0.05 were shown in Additional file 4 and mentioned in the text as “modest enrichment”. We further validated findings from IPA analyses from cohort 2 using overrepresentation analysis from the online software tool WebGestalt (www.webgestalt.org) [73] under default settings and summarized using REVIGO tool [74] based on multidimensional scaling of overrepresented GO terms with semantic similarities. STRING network was generated using STRING database version 10.5.

Brain sectioning and immunofluorescence

Brain blocks from 15 MS cases and controls (8 of them overlapping with DNA methylation analysis) were sectioned using a cryostat (Leica CM1850) at the Neurology clinic (Karolinska Hospital, Stockholm). The 14-μm-thick slices on SuperFrost slides were kept at − 20 degrees until further use. Immunofluorescence (IF) technique was utilized for examination of the transcription factor CREB in frozen MS and non-neurological control brain tissue. Rabbit monoclonal [E113] antibody specific for CREB phosphorylated on Serine 133 (1:200, Abcam #ab32096) was co-targeted with mouse monoclonal anti-NeuN antibody (Millipore #MAB377) in the samples with comparable (posterior) brain localization containing adjacent white and grey matter portions. Targets were visualized using fluorescently-labelled secondary antibodies Alexa Fluor 488 (Abcam, #ab150106) and 555 (Jackson Immunoresearch #111-545-003), respectively. NeuN+ white and grey matter neurons were examined for nuclear expression of phosphorylated CREB in high magnification using Zeiss LSM 700 confocal laser microscope. Captured group-representative IF- images are shown in the Fig. 5. Amount of NeuN+/CREB (phospho S133)+ WM neurons is presented as a percentage of the total amount of all detectable NeuN+ cells in the WM.

Statistical analysis

Details of genome-wide analysis of methylation data are provided in the sections above. All correlation analysis were performed using the Spearman’s rank test. Fisher’s and Chi-square test were used for enrichment and depletion analyses. Data with more than two groups were analysed using Kruskal-Wallis test and Dunn’s test for multiple comparisons.

Abbreviations

450K:

Illumina Infinium Human Methylation 450K BeadChip

5hmC:

DNA hydroxymethylation

5mC:

DNA methylation

BS:

Bisulfite

CNS:

Central nervous system

CREB:

cAMP response element binding

DhMP:

Differentially hydroxymethylated position

DhMR:

Differentially hydroxymethylated region

DMP:

Differentially methylated position

DMR:

Differentially methylated region

DNMT:

DNA methyltransferase

GABA:

Gamma-aminobutyric acid

GLU:

Glutamate

GM:

Grey matter

IF:

Immunofluorescence

MS:

Multiple sclerosis

NAWM:

Normal-appearing white matter

oxBS:

Oxidative bisulfite

TET:

Ten-eleven translocation

TSS:

Transcription starting site

UTR:

Untranslated region

WM:

White matter

References

  1. Compston A, Coles A. Multiple sclerosis. Lancet. 2008;372(9648):1502–17.

    Article  CAS  PubMed  Google Scholar 

  2. International Multiple Sclerosis Genetics C, Wellcome Trust Case Control C, Sawcer S, Hellenthal G, Pirinen M, Spencer CC, et al. Genetic risk and a primary role for cell-mediated immune mechanisms in multiple sclerosis. Nature. 2011;476(7359):214–9.

    Article  Google Scholar 

  3. Olsson T, Barcellos LF, Alfredsson L. Interactions between genetic, lifestyle and environmental risk factors for multiple sclerosis. Nat Rev Neurol. 2017;13(1):25–36.

    Article  CAS  PubMed  Google Scholar 

  4. de Groot V, Beckerman H, Uitdehaag BM, Hintzen RQ, Minneboo A, Heymans MW, et al. Physical and cognitive functioning after 3 years can be predicted using information from the diagnostic process in recently diagnosed multiple sclerosis. Arch Phys Med Rehabil. 2009;90(9):1478–88.

    Article  PubMed  Google Scholar 

  5. Barkhof F, Calabresi PA, Miller DH, Reingold SC. Imaging outcomes for neuroprotection and repair in multiple sclerosis trials. Nat Rev Neurol. 2009;5(5):256–66.

    Article  PubMed  Google Scholar 

  6. Dineen RA, Vilisaar J, Hlinka J, Bradshaw CM, Morgan PS, Constantinescu CS, et al. Disconnection as a mechanism for cognitive dysfunction in multiple sclerosis. Brain. 2009;132(Pt 1):239–49.

    Article  CAS  PubMed  Google Scholar 

  7. Francis PL, Chia TL, Jakubovic R, O’Connor P, Lee L, Feinstein A, et al. Extensive white matter dysfunction in cognitively impaired patients with secondary-progressive multiple sclerosis. AJNR Am J Neuroradiol. 2014;35(10):1910–5.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  8. Meijer KA, Muhlert N, Cercignani M, Sethi V, Ron MA, Thompson AJ, et al. White matter tract abnormalities are associated with cognitive dysfunction in secondary progressive multiple sclerosis. Mult Scler. 2016;22(11):1429–37.

    Article  PubMed  Google Scholar 

  9. Tallantyre EC, Bo L, Al-Rawashdeh O, Owens T, Polman CH, Lowe JS, et al. Clinico-pathological evidence that axonal loss underlies disability in progressive multiple sclerosis. Mult Scler. 2010;16(4):406–11.

    Article  PubMed  Google Scholar 

  10. Kalkers NF, Bergers E, Castelijns JA, van Walderveen MA, Bot JC, Ader HJ, et al. Optimizing the association between disability and biological markers in MS. Neurology. 2001;57(7):1253–8.

    Article  CAS  PubMed  Google Scholar 

  11. Popescu V, Agosta F, Hulst HE, Sluimer IC, Knol DL, Sormani MP, et al. Brain atrophy and lesion load predict long term disability in multiple sclerosis. J Neurol Neurosurg Psychiatry. 2013;84(10):1082–91.

    Article  PubMed  Google Scholar 

  12. Mahad DH, Trapp BD, Lassmann H. Pathological mechanisms in progressive multiple sclerosis. Lancet Neurol. 2015;14(2):183–93.

    Article  CAS  PubMed  Google Scholar 

  13. Sanchez-Mut JV, Heyn H, Vidal E, Delgado-Morales R, Moran S, Sayols S, et al. Whole genome grey and white matter DNA methylation profiles in dorsolateral prefrontal cortex. Synapse. 2017;71(6).

    Article  Google Scholar 

  14. Davies MN, Volta M, Pidsley R, Lunnon K, Dixit A, Lovestone S, et al. Functional annotation of the human brain methylome identifies tissue-specific epigenetic variation across brain and blood. Genome Biol. 2012;13(6):R43.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  15. Lister R, Mukamel EA, Nery JR, Urich M, Puddifoot CA, Johnson ND, et al. Global epigenomic reconfiguration during mammalian brain development. Science. 2013;341(6146):1237905.

    Article  PubMed  PubMed Central  Google Scholar 

  16. Zheleznyakova G, Piket E, Marabita F, Pahlevan Kakhki M, Ewing E, Ruhrmann S, et al. Epigenetic research in multiple sclerosis: progress, challenges and opportunities. Physiol Genomics. 2017;49(9):447–61.

    Article  CAS  PubMed  Google Scholar 

  17. Kozlenkov A, Wang M, Roussos P, Rudchenko S, Barbu M, Bibikova M, et al. Substantial DNA methylation differences between two major neuronal subtypes in human brain. Nucleic Acids Res. 2016;44(6):2593–612.

    Article  PubMed  Google Scholar 

  18. Tasic B, Yao Z, Graybuck LT, Smith KA, Nguyen TN, Bertagnolli D, et al. Shared and distinct transcriptomic cell types across neocortical areas. Nature. 2018;563(7729):72–8.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  19. Huynh JL, Garg P, Thin TH, Yoo S, Dutta R, Trapp BD, et al. Epigenome-wide differences in pathology-free regions of multiple sclerosis-affected brains. Nat Neurosci. 2014;17(1):121–30.

    Article  CAS  PubMed  Google Scholar 

  20. Mellen M, Ayata P, Dewell S, Kriaucionis S, Heintz N. MeCP2 binds to 5hmC enriched within active genes and accessible chromatin in the nervous system. Cell. 2012;151(7):1417–30.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  21. Colquitt BM, Allen WE, Barnea G, Lomvardas S. Alteration of genic 5-hydroxymethylcytosine patterning in olfactory neurons correlates with changes in gene expression and cell identity. Proc Natl Acad Sci U S A. 2013;110(36):14682–7.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  22. Stewart SK, Morris TJ, Guilhamon P, Bulstrode H, Bachman M, Balasubramanian S, et al. oxBS-450K: a method for analysing hydroxymethylation using 450K BeadChips. Methods0. 2015;72:9–15.

    CAS  Google Scholar 

  23. Houseman EA, Johnson KC, Christensen BC. OxyBS: estimation of 5-methylcytosine and 5-hydroxymethylcytosine from tandem-treated oxidative bisulfite and bisulfite DNA. Bioinformatics. 2016;32(16):2505–7.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  24. Neri F, Rapelli S, Krepelova A, Incarnato D, Parlato C, Basile G, et al. Intragenic DNA methylation prevents spurious transcription initiation. Nature. 2017;543(7643):72–7.

    Article  CAS  PubMed  Google Scholar 

  25. Mendizabal I, Zeng J, Keller TE, Yi SV. Body-hypomethylated human genes harbor extensive intragenic transcriptional activity and are prone to cancer-associated dysregulation. Nucleic Acids Res. 2017;45(8):4390–400.

    CAS  PubMed  PubMed Central  Google Scholar 

  26. Fan G, Beard C, Chen RZ, Csankovszki G, Sun Y, Siniaia M, et al. DNA hypomethylation perturbs the function and survival of CNS neurons in postnatal animals. J Neurosci. 2001;21(3):788–97.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  27. Rhee KD, Yu J, Zhao CY, Fan G, Yang XJ. Dnmt1-dependent DNA methylation is essential for photoreceptor terminal differentiation and retinal neuron survival. Cell Death Dis. 2012;3:e427.

    Article  PubMed  PubMed Central  Google Scholar 

  28. Hutnick LK, Golshani P, Namihira M, Xue Z, Matynia A, Yang XW, et al. DNA hypomethylation restricted to the murine forebrain induces cortical degeneration and impairs postnatal neuronal maturation. Hum Mol Genet. 2009;18(15):2875–88.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  29. Stampanoni Bassi M, Mori F, Buttari F, Marfia GA, Sancesario A, Centonze D, et al. Neurophysiology of synaptic functioning in multiple sclerosis. Clin Neurophysiol. 2017;128(7):1148–57.

    Article  PubMed  Google Scholar 

  30. Wei H, Feng Y, Liang F, Cheng W, Wu X, Zhou R, et al. Role of oxidative stress and DNA hydroxymethylation in the neurotoxicity of fine particulate matter. Toxicology. 2017;380:94–103.

    Article  CAS  PubMed  Google Scholar 

  31. Ramesh V, Bayam E, Cernilogar FM, Bonapace IM, Schulze M, Riemenschneider MJ, et al. Loss of Uhrf1 in neural stem cells leads to activation of retroviral elements and delayed neurodegeneration. Genes Dev. 2016;30(19):2199–212.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  32. Szulwach KE, Li X, Li Y, Song CX, Wu H, Dai Q, et al. 5-hmC-mediated epigenetic dynamics during postnatal neurodevelopment and aging. Nat Neurosci. 2011;14(12):1607–16.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  33. van Horssen J, van der Pol S, Nijland P, Amor S, Perron H. Human endogenous retrovirus W in brain lesions: rationale for targeted therapy in multiple sclerosis. Mult Scler Relat Disord. 2016;8:11–8.

    Article  PubMed  Google Scholar 

  34. Fagone P, Mangano K, Di Marco R, Touil-Boukoffa C, Chikovan T, Signorelli S, et al. Expression of DNA methylation genes in secondary progressive multiple sclerosis. J Neuroimmunol. 2016;290:66–9.

    Article  CAS  PubMed  Google Scholar 

  35. Mastronardi FG, Noor A, Wood DD, Paton T, Moscarello MA. Peptidyl argininedeiminase 2 CpG island in multiple sclerosis white matter is hypomethylated. J Neurosci Res. 2007;85(9):2006–16.

    Article  CAS  PubMed  Google Scholar 

  36. Singhal NK, Freeman E, Arning E, Wasek B, Clements R, Sheppard C, et al. Dysregulation of methionine metabolism in multiple sclerosis. Neurochem Int. 2018;112:1–4.

    Article  CAS  PubMed  Google Scholar 

  37. Beecham AH, Patsopoulos NA, Xifara DK, Davis MF, Kemppinen A, Cotsapas C, et al. Analysis of immune-related loci identifies 48 new susceptibility variants for multiple sclerosis. Nat Genet. 2013;45(11):1353–60.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  38. Werner P, Pitt D, Raine CS. Multiple sclerosis: altered glutamate homeostasis in lesions correlates with oligodendrocyte and axonal damage. Ann Neurol. 2001;50(2):169–80.

    Article  CAS  PubMed  Google Scholar 

  39. Clements RJ, McDonough J, Freeman EJ. Distribution of parvalbumin and calretinin immunoreactive interneurons in motor cortex from multiple sclerosis post-mortem tissue. Exp Brain Res. 2008;187(3):459–65.

    Article  CAS  PubMed  Google Scholar 

  40. Dutta R, Chang A, Doud MK, Kidd GJ, Ribaudo MV, Young EA, et al. Demyelination causes synaptic alterations in hippocampi from multiple sclerosis patients. Ann Neurol. 2011;69(3):445–54.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  41. Mandolesi G, Gentile A, Musella A, Fresegna D, De Vito F, Bullitta S, et al. Synaptopathy connects inflammation and neurodegeneration in multiple sclerosis. Nat Rev Neurol. 2015;11(12):711–24.

    Article  CAS  PubMed  Google Scholar 

  42. Azevedo CJ, Kornak J, Chu P, Sampat M, Okuda DT, Cree BA, et al. In vivo evidence of glutamate toxicity in multiple sclerosis. Ann Neurol. 2014;76(2):269–78.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  43. Baranzini SE, Srinivasan R, Khankhanian P, Okuda DT, Nelson SJ, Matthews PM, et al. Genetic variation influences glutamate concentrations in brains of patients with multiple sclerosis. Brain. 2010;133(9):2603–11.

    Article  PubMed  PubMed Central  Google Scholar 

  44. Gao Y, Deng K, Hou J, Bryson JB, Barco A, Nikulina E, et al. Activated CREB is sufficient to overcome inhibitors in myelin and promote spinal axon regeneration in vivo. Neuron. 2004;44(4):609–21.

    Article  CAS  PubMed  Google Scholar 

  45. Sakamoto K, Karelina K, Obrietan K. CREB: a multifaceted regulator of neuronal plasticity and protection. J Neurochem. 2011;116(1):1–9.

    Article  CAS  PubMed  Google Scholar 

  46. Lindberg RL, De Groot CJ, Certa U, Ravid R, Hoffmann F, Kappos L, et al. Multiple sclerosis as a generalized CNS disease--comparative microarray analysis of normal appearing white matter and lesions in secondary progressive MS. J Neuroimmunol. 2004;152(1-2):154–67.

    Article  CAS  PubMed  Google Scholar 

  47. Ren M, Han M, Wei X, Guo Y, Shi H, Zhang X, et al. FTY720 attenuates 6-OHDA-associated dopaminergic degeneration in cellular and mouse parkinsonian models. Neurochem Res. 2017;42(2):686–96.

    Article  CAS  PubMed  Google Scholar 

  48. Prins M, Schul E, Geurts J, van der Valk P, Drukarch B, van Dam AM. Pathological differences between white and grey matter multiple sclerosis lesions. Ann N Y Acad Sci. 2015;1351:99–113.

    Article  PubMed  Google Scholar 

  49. Lassmann H. Demyelination and neurodegeneration in multiple sclerosis: The role of hypoxia. Ann Neurol. 2016;79(4):520–1.

    Article  PubMed  Google Scholar 

  50. Loitfelder M, Fazekas F, Petrovic K, Fuchs S, Ropele S, Wallner-Blazek M, et al. Reorganization in cognitive networks with progression of multiple sclerosis: insights from fMRI. Neurology. 2011;76(6):526–33.

    Article  CAS  PubMed  Google Scholar 

  51. Chan D, Binks S, Nicholas JM, Frost C, Cardoso MJ, Ourselin S, et al. Effect of high-dose simvastatin on cognitive, neuropsychiatric, and health-related quality-of-life measures in secondary progressive multiple sclerosis: secondary analyses from the MS-STAT randomised, placebo-controlled trial. Lancet Neurol. 2017;16(8):591–600.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  52. Sedel F, Bernard D, Mock DM, Tourbah A. Targeting demyelination and virtual hypoxia with high-dose biotin as a treatment for progressive multiple sclerosis. Neuropharmacology. 2016;110(Pt B):644–53.

    Article  CAS  PubMed  Google Scholar 

  53. Matevossian A, Akbarian S. Neuronal nuclei isolation from human postmortem brain tissue. J Vis Exp. 2008;20.

  54. van Iterson M, Tobi EW, Slieker RC, den Hollander W, Luijk R, Slagboom PE, et al. MethylAid: visual and interactive quality control of large Illumina 450k datasets. Bioinformatics. 2014;30(23):3435–7.

    Article  PubMed  Google Scholar 

  55. Morris TJ, Butcher LM, Feber A, Teschendorff AE, Chakravarthy AR, Wojdacz TK, et al. ChAMP: 450k chip analysis methylation pipeline. Bioinformatics. 2014;30(3):428–30.

    Article  CAS  PubMed  Google Scholar 

  56. Aryee MJ, Jaffe AE, Corrada-Bravo H, Ladd-Acosta C, Feinberg AP, Hansen KD, et al. Minfi: a flexible and comprehensive Bioconductor package for the analysis of Infinium DNA methylation microarrays. Bioinformatics. 2014;30(10):1363–9.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  57. Nordlund J, Backlin CL, Wahlberg P, Busche S, Berglund EC, Eloranta ML, et al. Genome-wide signatures of differential DNA methylation in pediatric acute lymphoblastic leukemia. Genome Biol. 2013;14(9):r105.

    Article  PubMed  PubMed Central  Google Scholar 

  58. Chen YA, Lemire M, Choufani S, Butcher DT, Grafodatskaya D, Zanke BW, et al. Discovery of cross-reactive probes and polymorphic CpGs in the Illumina Infinium HumanMethylation450 microarray. Epigenetics. 2013;8(2):203–9.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  59. Teschendorff AE, Marabita F, Lechner M, Bartlett T, Tegner J, Gomez-Cabrero D, et al. A beta-mixture quantile normalization method for correcting probe design bias in Illumina Infinium 450k DNA methylation data. Bioinformatics. 2013;29(2):189–96.

    Article  CAS  PubMed  Google Scholar 

  60. Marabita F, Almgren M, Lindholm ME, Ruhrmann S, Fagerstrom-Billai F, Jagodic M, et al. An evaluation of analysis pipelines for DNA methylation profiling using the Illumina HumanMethylation450 BeadChip platform. Epigenetics. 2013;8(3):333–46.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  61. Maksimovic J, Gordon L, Oshlack A. SWAN: Subset-quantile within array normalization for illumina infinium HumanMethylation450 BeadChips. Genome Biol. 2012;13(6):R44.

    Article  PubMed  PubMed Central  Google Scholar 

  62. Johnson WE, Li C, Rabinovic A. Adjusting batch effects in microarray expression data using empirical Bayes methods. Biostatistics. 2007;8(1):118–27.

    Article  PubMed  Google Scholar 

  63. Houseman EA, Accomando WP, Koestler DC, Christensen BC, Marsit CJ, Nelson HH, et al. DNA methylation arrays as surrogate measures of cell mixture distribution. BMC bioinformatics. 2012;13:86.

    Article  PubMed  PubMed Central  Google Scholar 

  64. Teschendorff AE, Breeze CE, Zheng SC, Beck S. A comparison of reference-based algorithms for correcting cell-type heterogeneity in epigenome-wide association studies. BMC bioinformatics. 2017;18(1):105.

    Article  PubMed  PubMed Central  Google Scholar 

  65. Gu Z, Eils R, Schlesner M. Complex heatmaps reveal patterns and correlations in multidimensional genomic data. Bioinformatics. 2016;32(18):2847–9.

    Article  CAS  PubMed  Google Scholar 

  66. Ritchie ME, Phipson B, Wu D, Hu Y, Law CW, Shi W, et al. limma powers differential expression analyses for RNA-sequencing and microarray studies. Nucleic Acids Res. 2015;43(7):e47.

    Article  PubMed  PubMed Central  Google Scholar 

  67. Du P, Zhang X, Huang CC, Jafari N, Kibbe WA, Hou L, et al. Comparison of Beta-value and M-value methods for quantifying methylation levels by microarray analysis. BMC bioinformatics. 2010;11:587.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  68. Peters TJ, Buckley MJ, Statham AL, Pidsley R, Samaras K, R VL, et al. De novo identification of differentially methylated regions in the human genome. Epigenetics Chromatin. 2015;8:6.

    Article  PubMed  PubMed Central  Google Scholar 

  69. Willer CJ, Li Y, Abecasis GR. METAL: fast and efficient meta-analysis of genomewide association scans. Bioinformatics. 2010;26(17):2190–1.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  70. Gonnermann A, Framke T, Grosshennig A, Koch A. No solution yet for combining two independent studies in the presence of heterogeneity. Stat Med. 2015;34(16):2476–80.

    Article  PubMed  PubMed Central  Google Scholar 

  71. Benjamini Y, Hochberg Y. Controlling the false discovery rate: a practical and powerful approach to multiple testing. J Royal Stat Soc Series B (Methodological). 1995;57(1):289–300.

    Article  Google Scholar 

  72. Tokar T, Pastrello C, Rossos AEM, Abovsky M, Hauschild AC, Tsay M, et al. mirDIP 4.1-integrative database of human microRNA target predictions. Nucleic Acids Res. 2018;46(D1):D360–D70.

    Article  CAS  PubMed  Google Scholar 

  73. Zhang B, Kirov S, Snoddy J. WebGestalt: an integrated system for exploring gene sets in various biological contexts. Nucleic Acids Res. 2005;33(Web Server issue):W741–8.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  74. Supek F, Bosnjak M, Skunca N, Smuc T. REVIGO summarizes and visualizes long lists of gene ontology terms. PLoS One. 2011;6(7):e21800.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

Download references

Acknowledgements

We are grateful to A. van Vollenhoven for flow cytometry processing (Center for Molecular Medicine and Karolinska Institutet core facility), D. Fauvin for assistance with TrueMethyl and T. Morris for support in the BS/oxBS data pipeline. We thank R. Covacu for helpful feedback on data interpretation. We acknowledge GenomeScan/ServiceXS (Leiden, The Netherlands) for processing Illumina 450K array BS and oxBS data (cohort 2) and BEA core facility (Karolinska Institutet) for processing Illumina 450K array BS data (pilot samples and cohort 1). We thank the Multiple Sclerosis and Parkinson’s Tissue Bank (Imperial College London) for provision of brain tissue samples. Computations were performed on resources provided by SNIC through Uppsala Multidisciplinary Center for Advanced Computational Science (UPPMAX).

Funding

This work was supported by grants from the Swedish Research Council, the Swedish Association for Persons with Neurological Disabilities, the Swedish Brain Foundation, the Swedish MS Foundation, Åsa Vilhelmsson Foundation, Petrus and Augusta Hedlunds Foundation, the Stockholm County Council (ALF project) and Karolinska Institutet. L. Kular was supported by fellowship from the Margaretha af Ugglas Foundation. D. Gomez-Cabrero and J. Tegner were supported by EU FP7 306000 STATegra.

Availability of data and materials

The custom pipelines used in this study are detailed in the “Methods” section and specific custom codes are available from the corresponding author upon request. The 450K data that support the findings of this study are available in Gene Expression Omnibus (GEO) database under the accession number GSE119532.

Author information

Authors and Affiliations

  1. Department of Clinical Neuroscience, Center for Molecular Medicine, Karolinska Institutet, Stockholm, Sweden

    Lara Kular, Maria Needhamsen, Milena Z. Adzemovic, Tatiana Kramarova, Ewoud Ewing, Eliane Piket, Fredrik Piehl, Lou Brundin & Maja Jagodic

  2. Department of Medicine, Unit of Computational Medicine, Center for Molecular Medicine, Karolinska Institutet, Solna, Sweden

    David Gomez-Cabrero & Jesper Tegnér

  3. Mucosal and Salivary Biology Division, King’s College London Dental Institute, London, SE1 9RT, UK

    David Gomez-Cabrero

  4. Translational Bioinformatics Unit, Navarrabiomed, Complejo Hospitalario de Navarra (CHN), Universidad Pública de Navarra (UPNA), IdiSNA, Pamplona, Spain

    David Gomez-Cabrero

  5. Biological and Environmental Sciences and Engineering Division, Computer, Electrical and Mathematical Sciences and Engineering Division, King Abdullah University of Science and Technology, Thuwal, 23955, Saudi Arabia

    Jesper Tegnér

  6. Medical Genomics, UCL Cancer Institute, University College London, London, UK

    Stephan Beck

  7. Department of Neurology, Karolinska University Hospital, Stockholm, Sweden

    Fredrik Piehl & Lou Brundin

Authors
  1. Lara Kular
    View author publications

    You can also search for this author in PubMed Google Scholar

  2. Maria Needhamsen
    View author publications

    You can also search for this author in PubMed Google Scholar

  3. Milena Z. Adzemovic
    View author publications

    You can also search for this author in PubMed Google Scholar

  4. Tatiana Kramarova
    View author publications

    You can also search for this author in PubMed Google Scholar

  5. David Gomez-Cabrero
    View author publications

    You can also search for this author in PubMed Google Scholar

  6. Ewoud Ewing
    View author publications

    You can also search for this author in PubMed Google Scholar

  7. Eliane Piket
    View author publications

    You can also search for this author in PubMed Google Scholar

  8. Jesper Tegnér
    View author publications

    You can also search for this author in PubMed Google Scholar

  9. Stephan Beck
    View author publications

    You can also search for this author in PubMed Google Scholar

  10. Fredrik Piehl
    View author publications

    You can also search for this author in PubMed Google Scholar

  11. Lou Brundin
    View author publications

    You can also search for this author in PubMed Google Scholar

  12. Maja Jagodic
    View author publications

    You can also search for this author in PubMed Google Scholar

Contributions

MJ conceived, designed and supervised the study. LK, MZA and TK conducted the experiments. MN, LK, MJ, DG-C, EE and EP performed analyses. SB and JT assisted in data interpretation. LB aided in accessing the samples, techniques for nuclear isolation and scientific discussions. FP provided insight into clinical information. LK and MN wrote the manuscript with assistance from all authors. All authors read and approved the manuscript.

Corresponding author

Correspondence to Maja Jagodic.

Ethics declarations

Ethics approval and consent to participate

Brain tissues used for this study and obtained from the Multiple Sclerosis and Parkinson’s Tissue Bank (Imperial College London) were approved by local ethical guidelines. All research included in this manuscript conforms with the Declaration of Heksinki.

Consent for publication

Not applicable

Competing interests

The authors declare that they have no competing interests.

Publisher’s Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Additional files

Additional file 1:

Figure S1. Importance of separating neuronal fraction from brain samples. Figure S2. Neuronal subtype deconvolution. Figure S3. Signal distribution of BS and oxBS probes. Figure S4. Distribution of 5hmC β-values. Figure S5. Computational strategy of 5mC and 5hmC analyses. Figure S6. Comparison of 5hmC β-values between pipelines. Figure S7.Validation of OBSCN locus using methyl-sensitive (MSRE) and glucosyl-sensitive (GSRE) restriction enzymes. Figure S8. Correlation between BS and oxBS changes. Figure S9. Distribution of 5mC and 5hmC changes across features. Figure S10. Meta-analysis. Figure S11. Neuroanatomical localization of the samples. (PDF 2883 kb)

Additional file 2:

Table S1. Description of cohorts. (XLSX 102 kb)

Additional file 3:

Table S2. BS-DMPs (cohort 1, P-value <0.001). (XLSX 197 kb)

Additional file 4:

Table S3. Canonical pathways and upstream regulators from gene ontology analysis for cohort 1, cohort 2 and meta-analysis. (XLSX 244 kb)

Additional file 5:

Table S4. 5mC-DMPs, 5hmC-DhMPs and BS-DMPs (adjusted P-Value <0.05, cohort 2). (XLSX 1370 kb)

Additional file 6:

Table S5. 5mC-DMRs, 5hmC-DhMRs and BS-DMRs (cohort 2). (XLSX 508 kb)

Additional file 7:

Table S6. DMPs from meta-analysis of cohorts 1 and 2 (adjusted P-Value <0.05). (XLSX 1560 kb)

Additional file 8:

Table S7. Comparison of DNA methylation with gene expression (Huynh et al., 2014). (XLSX 24 kb)

Rights and permissions

Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated.

Reprints and permissions

About this article

Check for updates. Verify currency and authenticity via CrossMark

Cite this article

Kular, L., Needhamsen, M., Adzemovic, M.Z. et al. Neuronal methylome reveals CREB-associated neuro-axonal impairment in multiple sclerosis. Clin Epigenet 11, 86 (2019). https://doi.org/10.1186/s13148-019-0678-1

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: https://doi.org/10.1186/s13148-019-0678-1

Keywords

Clinical Epigenetics

ISSN: 1868-7083

Contact us