- Open Access
Neuronal methylome reveals CREB-associated neuro-axonal impairment in multiple sclerosis
Clinical Epigeneticsvolume 11, Article number: 86 (2019)
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.
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.
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.
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 . Although the cause of MS remains largely unknown, accumulating data support the notion of MS being a complex disease influenced by genetic  (primarily affecting immune genes) and environmental factors . 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 . Apart from focal lesions, neuroimaging  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 . 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 . 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 . Moreover, the brain and more specifically neurons are enriched in regulatory non-canonical DNA modifications such as CpG hydroxymethylation (5hmC) . A growing body of evidence suggests that in addition to its crucial role in brain development and function , DNA methylation could be implicated in MS brain pathology as well . 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.
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  and found that neuronal subtype proportion indeed had an impact on DNA methylation profiles (see the “Methods” section and Additional file 1: Figure S2).
We next aimed to investigate DNA methylation changes underpinning neuronal vulnerability in MS. Considering the extensive heterogeneity of GM methylome , likely due to pervasive phenotypic diversity of GM excitatory neurons , 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  (Additional file 3).
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 , 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) , as 5hmC β-values generated by the two methods strongly correlated (Additional file 1: Figure S6).
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).
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.
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 . 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).
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).
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 . 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  (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  (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.
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 , 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  and defects in neurite outgrowth and synaptic formation. 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 . Importantly, strong transcriptional activity of human endogenous retroviruses, including in CNS cells, has been detected in MS patients . 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  and brain  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 . 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 , 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 . Interestingly, CREB activation differed between WM and GM neurons, likely reflecting differences in neuronal subtypes and/or distinct susceptibility towards demyelination pathology , such as hypoxia-like metabolic injury . 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 . The aberrations in WM-neurons of MS patients identified here are interesting in light of neuronal circuit disconnection  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].
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.
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.
Fluorescence-activated cell sorting (FACS)-based neuronal nuclei isolation from dissected brain tissue was performed according to previously published protocol  (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
450K data (485,577 probes) were quality assessed using MethylAid , 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  and minfi version 1.16.0  Bioconductor packages.
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.  (pilot study) and Chen et al.  (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) as previously recommended  in ChAMP for the pilot study or within-sample normalized using the “Subset-quantile Within Array Normalization” (SWAN) method as previously recommended for oxBS-450K data  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.
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)  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 . For comparison, we tested another method based on maximum likelihood estimation (MLE) available through the oxyBS version 1.0 Cran package  with default settings.
Correction for slide-effect
Slide effects were corrected using empirical Bayes methods  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  and filtered for |Δβ| > 0.7 which resulted in a total of 162 neuronal subtype-specific sites. In concordance with previous observations , 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  as implemented in the projectCellType() function of the minfi Bioconductor package . 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  (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 . 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  was used for detection of DMPs with M-values (Mi = log2(βi/(1 − βi))) as input as previously recommended . 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 , which identifies DMRs based on kernel smoothing, was applied with default settings (λ = 1000, C = 2).
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 . 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 . P values were adjusted for multiple comparisons using the Benjamini & Hochberg (B-H) method .
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.
“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 22.214.171.124). 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 . 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 . 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)  under default settings and summarized using REVIGO tool  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.
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.
Illumina Infinium Human Methylation 450K BeadChip
Central nervous system
cAMP response element binding
Differentially hydroxymethylated position
Differentially hydroxymethylated region
Differentially methylated position
Differentially methylated region
Normal-appearing white matter
Transcription starting site
Compston A, Coles A. Multiple sclerosis. Lancet. 2008;372(9648):1502–17.
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.
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.
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.
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.
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.
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.
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.
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.
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.
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.
Mahad DH, Trapp BD, Lassmann H. Pathological mechanisms in progressive multiple sclerosis. Lancet Neurol. 2015;14(2):183–93.
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).
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.
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.
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.
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.
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.
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.
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.
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.
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.
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.
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.
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.
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.
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.
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.
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.
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.
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.
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.
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.
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.
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.
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.
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.
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.
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.
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.
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.
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.
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.
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.
Sakamoto K, Karelina K, Obrietan K. CREB: a multifaceted regulator of neuronal plasticity and protection. J Neurochem. 2011;116(1):1–9.
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.
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.
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.
Lassmann H. Demyelination and neurodegeneration in multiple sclerosis: The role of hypoxia. Ann Neurol. 2016;79(4):520–1.
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.
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.
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.
Matevossian A, Akbarian S. Neuronal nuclei isolation from human postmortem brain tissue. J Vis Exp. 2008;20.
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.
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.
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.
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.
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.
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.
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.
Maksimovic J, Gordon L, Oshlack A. SWAN: Subset-quantile within array normalization for illumina infinium HumanMethylation450 BeadChips. Genome Biol. 2012;13(6):R44.
Johnson WE, Li C, Rabinovic A. Adjusting batch effects in microarray expression data using empirical Bayes methods. Biostatistics. 2007;8(1):118–27.
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.
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.
Gu Z, Eils R, Schlesner M. Complex heatmaps reveal patterns and correlations in multidimensional genomic data. Bioinformatics. 2016;32(18):2847–9.
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.
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.
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.
Willer CJ, Li Y, Abecasis GR. METAL: fast and efficient meta-analysis of genomewide association scans. Bioinformatics. 2010;26(17):2190–1.
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.
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.
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.
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.
Supek F, Bosnjak M, Skunca N, Smuc T. REVIGO summarizes and visualizes long lists of gene ontology terms. PLoS One. 2011;6(7):e21800.
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).
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.
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
The authors declare that they have no competing interests.
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
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)
Table S1. Description of cohorts. (XLSX 102 kb)
Table S2. BS-DMPs (cohort 1, P-value <0.001). (XLSX 197 kb)
Table S3. Canonical pathways and upstream regulators from gene ontology analysis for cohort 1, cohort 2 and meta-analysis. (XLSX 244 kb)
Table S4. 5mC-DMPs, 5hmC-DhMPs and BS-DMPs (adjusted P-Value <0.05, cohort 2). (XLSX 1370 kb)
Table S5. 5mC-DMRs, 5hmC-DhMRs and BS-DMRs (cohort 2). (XLSX 508 kb)
Table S6. DMPs from meta-analysis of cohorts 1 and 2 (adjusted P-Value <0.05). (XLSX 1560 kb)
Table S7. Comparison of DNA methylation with gene expression (Huynh et al., 2014). (XLSX 24 kb)