DNA methylation changes in glial cells of the normal-appearing white matter in Multiple Sclerosis patients

ABSTRACT

Multiple Sclerosis (MS), the leading cause of non-traumatic neurological disability in young adults, is a chronic inflammatory and neurodegenerative disease of the central nervous system (CNS). Due to the poor accessibility to the target organ, CNS-confined processes underpinning the later progressive form of MS remain elusive thereby limiting treatment options. We aimed to examine DNA methylation, a stable epigenetic mark of genome activity, in glial cells to capture relevant molecular changes underlying MS neuropathology. We profiled DNA methylation in nuclei of non-neuronal cells, isolated from 38 post-mortem normal-appearing white matter (NAWM) specimens of MS patients (n = 8) in comparison to white matter of control individuals (n = 14), using Infinium MethylationEPIC BeadChip. We identified 1,226 significant (genome-wide adjusted P-value < 0.05) differentially methylated positions (DMPs) between MS patients and controls. Functional annotation of the altered DMP-genes uncovered alterations of processes related to cellular motility, cytoskeleton dynamics, metabolic processes, synaptic support, neuroinflammation and signaling, such as Wnt and TGF-β pathways. A fraction of the affected genes displayed transcriptional differences in the brain of MS patients, as reported by publically available transcriptomic data. Cell type-restricted annotation of DMP-genes attributed alterations of cytoskeleton rearrangement and extracellular matrix remodelling to all glial cell types, while some processes, including ion transport, Wnt/TGF-β signaling and immune processes were more specifically linked to oligodendrocytes, astrocytes and microglial cells, respectively. Our findings strongly suggest that NAWM glial cells are highly altered, even in the absence of lesional insult, collectively exhibiting a multicellular reaction in response to diffuse inflammation.

KEYWORDS:

• Multiple Sclerosis
• glial cells
• DNA methylation
• motility
• Wnt
• TGF
• neuromodulation

Introduction

Multiple sclerosis (MS) is a chronic inflammatory demyelinating and neurodegenerative disease of the central nervous system (CNS). MS pathology presents with the occurrence of immune-induced demyelinating lesions arising throughout the CNS and translating into various neurological symptoms. The nature of injuries in MS is highly heterogeneous as lesions observe spatio-temporal diversity, i.e., varying across the CNS and disease stages [1]. Disability closely mirrors neuro-axonal deterioration, irrespective of the type of lesions and disease course [2]. While inflammation is unambiguously observed at all stages of the disease, compartmentalized chronic inflammation orchestrated by resident CNS cells dominate later stages of the disease, independently of immune infiltrates or pre-existing demyelination and accounts for the clinical trajectory [3–5]. Ultimately, persistent inflammation and recurrent damage to the myelin insulating axons will eventually exhaust the repair capacity of the CNS [6,7], a weakened functional recovery presaging transition to the progressive stage of the disease. Whereas major progress has been achieved in understanding and treating early phase of disease development, via targeting of peripheral immune cells, the CNS-confined mechanisms underlying the later stage of disease progression remain elusive. This is likely due to the difficulty to gather molecular evidence from the affected tissue itself, brain specimens being accessible post-mortem thereby restricting methodological applications in relatively small case–control cohorts of high-quality samples. This knowledge gap considerably cripples the care of progressive MS forms, leaving clinicians with a scarcity of therapeutic solutions and patients with relentless and untreatable disabilites.

The study of glial cell populations has gained particular interest due to the possibility to decipher the neurotoxic and pro-inflammatory processes precipitating neuronal vulnerability in progressive MS and thereby, to ultimately restore the brain repair capacities favouring remyelination and neuroprotection [8]. Due to their myelin-producing ability, oligodendrocytes provide a structural and trophic support to neurons by controlling the myelination of axons and are particularly abundant in the white matter (WM) compared to the grey matter (GM) [9]. Astrocytes are key homoeostatic regulators of the CNS and, as such, exert versatile functions in synapse refinement, neurotransmission, formation of the blood-brain-barrier and metabolic control of the microenvironment, among others, depending on their primary location. Microglia, populating less than 10% of glial cells, are highly responsive CNS-resident macrophages that assume various functions pertaining to their immunocompetent and phagocytic capacities, including but not limited to synaptic pruning, immunosurveillance and clearance. Due to their inherent migratory abilities, microglial cells are highly dynamic and vigilant cells at resting state, continuously surveying their microenvironment [10]. Overall, the CNS homoeostasis entails tightly controlled region-specific (e.g., WM vs. GM) regulations of cellular phenotypes through multidirectional glia-glia and neuro-glia signaling.

Investigation of the rodent brain in MS-like models and histopathological characterization of post-mortem brain sections of patients have been instrumental in unveiling the overlapping sequences of events occurring in the MS CNS, mostly under lesional immune insult. Emerging evidence from single-cell transcriptome studies further support CNS damage in MS to likely ensue from a complex interplay between various glial cell populations, with each cell type manifesting marked spatial and temporal phenotypic heterogeneity [11–15]. Importantly, diffuse abnormalities in myelination and neuroinflammation accumulate outside of the affected areas as well, namely in the Normal Appearing White Matter (NAWM) [16–19]. These changes reflect phenotypic dysfunction of glial cells in the absence of infiltrating leukocytes and macroscopic lesion and have been proposed as prominent early processes preceding newly forming lesions [20–27]. Importantly, alterations of the NAWM have been associated with cortical axonal loss, cognitive decline and clinical disability [28,29]. Yet, the molecular changes occurring in the NAWM remain elusive, most studies capturing highly dynamic transcriptional states of cells from MS lesions. In that regard, exploring the molecular layer of epigenetic changes, which orchestrate both durable and transitional states, could provide additional insight into the mechanisms underpinning brain pathology in progressive MS [30].

DNA methylation, the most studied epigenetic mark, relies on the stable deposition of a methyl group onto cytosine, primarily in a CpG context, and exerts regulatory action on gene expression that depends on its gene location insofar as methylation of gene promoter is associated with transcriptional repression while gene body methylation is likely connected to increased gene transcription [31,32]. Profiling DNA methylation genome-wide at single-base resolution enables probing the chromatin state and genome activity reliably in post-mortem tissue. Comparison of demyelinated and myelinated hippocampi of MS patients has identified aberrant epigenetic changes associated with deregulation of a small set of genes [33]. Previous case–control methylome analysis of bulk brain tissue has revealed numerous subtle changes in the bulk NAWM of MS patients compared to controls [34]. More recently, by conducting DNA methylation analysis of neuronal nuclei, we have identified functionally relevant changes, including reduced CREB transcription factors activity, associated with neuro-axonal impairment in MS patients compared to controls [35]. Here we aimed to elucidate the molecular alterations occurring in non-neuronal nuclei sorted from the NAWM of MS patients in comparison to the WM of non-neurological disease control individuals.

Materials and methods

Subjects, cohorts and ethics

Brain tissue used in this study was obtained from the Multiple Sclerosis and Parkinson’s Tissue Bank (Imperial College London), approved by local ethical guidelines. All research included in this manuscript conforms with the Declaration of Heksinki. The material comprises 38 snap-frozen brain tissue blocks collected within 33 h post-mortem from NAWM tissue of progressive MS patients (n = 8) and WM tissue of controls (n = 14) (Table 1). Further details are given in Supplementary Table 1. Control subjects were selected based on a non-neurological cause of death. The samples were further annotated according to brain location following antero-posterior axis and characteristics of the tissue using the human brain atlas sectional anatomy database (http://www.thehumanbrain.info) prior to dissection. Of note, different brain samples, coming from distinct regions, were used from the same individuals: two brain specimens per individual were used for 7 out of 8 MS and 5 out of 14 NNC individuals and three samples per individual were used for additional two NNC individuals. Samples reaching the following inclusion criteria have been included in DNA methylation analysis: (i) all available samples with sufficient DNA amount from WM non-neuronal nuclei, (ii) samples that passed DNA methylation quality control and (iii) cases with confirmed MS diagnosis and non-neurological controls without any signs of inflammation in the CNS.

Table 1. Description of the cohort.

Group	Cases (MS)	Controls (NNC)
N individuals	8	14
Sex ratio (F/M)	8:0^#	4:10
Age (mean ± SD)	63.1 ± 12.6*	76.4 ± 13.4
PMI (mean ± SD)^a	15.6 ± 8.6	21.8 ± 7.7
N brain samples Type lesion (N)	15 NAWM (15)	23 -
Fraction 1 (% ± SD) Fraction 2 (% ± SD) Fraction 3 (% ± SD)	36.6 ± 28.2* 39.2 ± 28.5* 23.8 ± 23.8	61.1 ± 24.3 11.0 ± 10.8 27.3 ± 23.6

MS, Multiple Sclerosis, NNC, non-neurological disease controls, N, number, F/M, female/male, PMI, post-mortem interval, NAWM, normal-appearing white matter, ^aPMI not available for 4 control individuals. Cell fractions were determined by RefFreeEwas cell type deconvolution. ^# P-value < 0.05 between cases and controls by chi-square (Χ²) test, * P-value < 0.05 between cases and controls generated by two-tailed independent samples T test.

Sample preparation

Fluorescence activated cell sorting (FACS)-based neuronal nuclei isolation from dissected brain tissue was performed as previously described [35]. Briefly, following resuspension of the homogenized brain tissue in hypotonic lysis buffer (0.32 M sucrose, 5 mM CaCl₂, 3 mM MgAc₂, 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 MgAc₂, 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 (MoFlo^TM high-speed cell sorter), a representative experiment is shown in Supplementary Figure 1. Non-neuronal nuclei (representing approximately 75%-95% of the sorted fractions) 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 until further use.

Illumina Infinium Human MethylationEPIC

We used Illumina Infinium Human MethylationEPIC BeadChip (Illumina, Inc., San Diego, CA, U.S.A, EPIC) for quantitative and genome-wide DNA methylation profiling. Genomic DNA samples 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 EPIC arrays were blinded to the MS disease status during the experiments.

DNA methylation analysis

The analytical pipeline is illustrated in Supplementary Figures 1 and 2.

Quality control. EPIC data was quality assessed using QC report from the minfi package. All samples passed quality control and were subsequently processed using the Chip Analysis Methylation Pipeline (ChAMP) version 2.9.10 [36] and minfi version 1.24.0 [37] R-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 European population) [38,39] and cross-reactivity as identified by Nordlund et al [40]. and Chen et al [41]. After filtering for probes located on X and Y chromosomes, 700,482 probes remained.

Between and within-array normalization. Probes were subjected to within-sample normalization (Noob), which corrects for two different probe designs (type I and type II probes) included on the EPIC BeadChip. Slide effects, as identified using Principal Component Analysis (PCA), were corrected using empirical Bayes methods [42] implemented in the ComBat function of the SVA Bioconductor package version 3.26.0.

Deconvolution. Reference-free cell-type deconvolution was performed using RefFreeEWAS [43] version 2.2, as no accurate cell-based reference models exist for deconvolution of DNA methylation in the glial fraction. The optimal number of fractions for deconvolution was determined by calculating the deviance-boots (epsilon value) over 100 iterations for the range of 1 to 6 fractions, with minimum deviance with 3 fractions. The fractions were obtained by solving the model $Y=M\times \mathrm{\Omega }-T$ (where Y = original beta methylation matrix, M = cell-type specific beta methylation matrix, Ω = cell proportion matrix and T is number of cell types to deconvolute) using the non-negative matrix factorization method where the fractions represent the estimate proportion of cells in the mixture belonging to this cell type, which was used as a covariate in the linear model.

Differentially methylated positions (DMPs) and regions (DMRs). The Limma Bioconductor package version 3.34.9 [44] was used for detection of DMPs with M-values ($M_i={log}_2\left(\frac{{\beta }_i}{1-{\beta }_i}\right)$) as input as previously recommended [45]. Since several samples with different brain location were used from the same donor, Limma was conducted using within-individual comparions in a linear mixed model (LMM) by estimating the average correlation within individual prior to using the function duplicateCorrelation which was then used in the lmFit step. The disease was tested as an exposure, empirical Bayes (eBayes) was used to moderate the standard errors towards a common value, P-values were corrected using a Benjamini-Hochberg P-value correction. The following covariates were included in the model: sex, age and cell type proportions. DMRcate version 1.6.53 [46], which identifies DMRs based on kernel smoothing, was applied with default settings (λ = 1000, C = 2). DMPs and DMRs are presented in Supplementary Tables 2 and 3, respectively.

Gene annotation. Classical EPIC annotations (TSS200, TSS1500, 1stExon, 5ʹUTR, Gene body, 3ʹUTR, CGI, Shelf, Shore and Open Sea) as well as Fantom-annotated enhancers were derived from the ”IlluminaHumanMethylationEPICanno.ilm10b2.hg19” package version 0.6.0. CpG Islands (CGIs) were defined as GC content > 50%, observed/expected CpG ratio > 60%, > 200bp, while CGI shores and shelves represent 2kb flanking regions within or outside CGIs, 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.

Transcription factor enrichment analysis

We addressed putative enrichment of transcription factor (TF) binding sites at the DMP-related sequences by mapping genomic location of recognition sequence to the CpGs targeted by the EPIC probes. We downloaded genomic coordinates (hg19) from JASPAR CORE database 2022 (http://expdata.cmmt.ubc.ca/JASPAR/downloads/UCSC_tracks/2022/) as a bed file, which we intersected, using BedIntersect function from bedtools, with a bed file of the EPIC design. We could assign transcription factor binding site to each CpG of the EPIC array (including our 700,482 probes), but only if the CpG base falls directly in the binding site. We next performed statistical enrichment testing using Fisher’s exact test to test for an enrichment of TFs binding in our DMP list (Supplementary Table 4).

Annotation of DMPs

Annotations of the DMPs are presented in Supplementary Table 2.

Age, sex, methylation Quantitative Trait Locus (meQTL). We assessed potential genetic influence on the identified changes by annotating the DMPs according to CNS-specific meQTL data from bulk prefrontal cortex samples (n = 468, http://mostafavilab.stat.ubc.ca/xQTLServe, version 2021) [47] and NeuN-negative glial fraction sorted from temporal cortex tissue (n = 22) [48]. A total of 280 and 9 DMPs were annotated as meQTL-CpGs from bulk and glial meQTL studies, respectively. Moreover, we marked DMPs that were previously reported as cross-tissue sex-associated probes [49] as well DMPs that display an absolute Spearman correlation value > 0.8 between methylation (fitted β-value) and the age of individuals in our cohort. Residual associations to sex and age could be seen for 51 and 14 probes, respectively, out of 1,226 DMPs. Of note, because findings from the aforementioned studies were generated by Illumina Infinium Human Methylation 450K BeadChip (450 K), a fraction of the identified DMPs could not be mapped. Therefore, since 685 of the 1,126 DMPs were only present in the 450K array, the bulk meQTL-, glia meQTL-, sex- and age-associated fractions represent 40%, 1.3%, 7% and 2% of the 450K-annotated DMPs, respectively. We performed statistical enrichment testing of meQTL-CpGs in our 450K-annotated DMP list using Fisher’s exact test.

Single-cell transcriptomics. We explored cell type-specific alterations at the 687 DMP-genes by annotating DMP-genes according to single-cell and cell type-immunopanned transcriptomic profiles of the healthy human brain [50,51]. We performed hierarchical clustering on each dataset separately (Supplementary Fig. 5, Supplementary Table 6). Given the discrepancies between findings of these two studies, likely due to the distinct methodologies employed, we selected consensus cell type-specific gene sets based on the union of these two independent analyses, i.e., group of DMP-containing genes with stronger basal expression in one of the three cell types compared to the others. Of note, 143 genes strongly expressed in more than one cell type were included in each respective cell type-specific gene sets and therefore overlap between cell type-specific groups and cell type-specific constitutive expression of 177 genes was unavailable. To further address cell type-specific differential expression in MS brain, we annotated our DMP-genes according to findings from three studies published to date reporting single-cell RNA-seq findings from MS lesions versus control brain tissue. The first study profiled single nuclei from cortical GM (with attached meninges) and subcortical WM of MS-lesions (n = 10) compared to control tissue (n = 9) [13]. The second study compared single-cell transcriptomes from WM samples (MS lesions and NAWM) of MS patients (n = 4) and WM of NNC individuals (n = 5) [14]. The third study profiled the edge of demyelinated WM lesions at various stages of inflammation of MS patients (n = 5) and WM samples of NNC (n = 3) [52].

Gene ontology analyses

Gene ontology (GO) analysis of the DMP-genes (P_adj < 0.05) was performed using overrepresentation analysis (ORA) from the online software tool WebGestalt (www.webgestalt.org)[53] under default settings. Findings from ORA were clustered using REVIGO (http://revigo.irb.hr/)[54] or GeneSetCluster package [55] version 1.2.1. STRING network analysis applied to DMPs (P_adj < 0.05) was generated using STRING database version 11.0.

Results

DNA methylation changes in non-neuronal cells from MS patients

We performed genome-wide DNA methylation profiling on bisulfite (BS)-treated genomic DNA isolated from NeuN-negative (hereafter referred to as glial) fraction sorted from the WM of fresh-frozen post-mortem tissue blocks. A total of 38 glial nuclei samples isolated from NAWM tissue of 8 MS patients and WM of 14 non-neurological disease controls (NNC) were profiled using Illumina HumanMethylationEPIC BeadChip (Table 1, Supplementary Table 1). We conducted reference-free deconvolution to account for varying proportions of cell types, which overcomes the current lack of existing reference methylome for distinct human glial cell types, using RefFreeEWAS tool [43]. This method identified optimal separation to be based on three fractions, of which two differ (P < 0.05) between MS and NNC samples (Table 1, Supplementary Fig. 2).

After correction for confounders, we identified 1,226 differentially methylation positions (DMPs) mapping to 687 annotated genes, between MS and NNC (adjusted P-value, P_adj < 0.05) (Figure 1a,b). Most DMPs (65%, 803/1226) exhibited differences exceeding |Δβ| > 0.05; the most significant ones are listed in Table 2 (complete data is presented in Supplementary Table 2). Slightly more than half of the identified changes (55%, 682/1226) exhibited hypermethylation in MS patients compared to controls. Markedly, a substantial fraction (44%) of the identified DMPs (541/1226) are unique to the Illumina EPIC array, hence not present in the previous Illumina 450K array. These DMPs represented ~45% (327/687) of the annotated genes, with minor overlap (10/327) with 450K-annotated DMP-genes. Of note, only three and seven DMPs were previously reported in MS vs. NNC bulk brain [34] and neurons [35], respectively (Supplementary Fig. 3). Annotation of the DMPs according to CNS-restricted meQTL data [47,48] indicated that methylation at nine DMPs, mapping to BEGAIN, TBKBP1, RPTOR, TRIO, SLC15A4, CCDC64B and C17orf57 genes, are controlled by genetic variation specifically in glial cells, with additional 271 DMPs exhibiting meQTL effect in bulk cortical tissue (Supplementary Table 2). Fisher’s exact test indicated significant enrichment of glia- and bulk brain-restricted meQTL CpGs within the DMPs (P = 1.06 x 10⁻⁰⁵ and P = 2.20 x 10⁻¹⁶, respectively). We examined whether changes clustered in differentially methylated regions (DMRs) and identified 276 DMRs mapping to 225 gene-annotated loci between MS patients and controls (mean Δβ_DMR ranging from −0.18 to 0.19). The majority of the DMRs (71%, 197/276) encompassed at least one DMP (Supplementary Table 3).

Table 2. Top significant differentially methylated positions associated with Multiple Sclerosis in glial nuclei.

Probe	Chr.	Position	Gene	Feature	cgi	MeanMS	MeanNNC	Δβ	P	P.adj
cg20810195	1	86,970,291		IGR	opensea	0.738	0.868	−0.114	1.7E-11	4.6E-07
cg09674340	1	202,509,286	PPP1R12B	1stExon	opensea	0.933	0.845	0.091	1.5E-15	6.4E-11
cg26919182	1	202,522,232	PPP1R12B	Body	opensea	0.935	0.598	0.337	3.7E-33	2.6E-27
cg15817705	1	209,406,063		IGR	shore	0.737	0.849	−0.104	3.1E-11	8.2E-07
cg12691488	1	243,053,673		IGR	Island	0.134	0.297	−0.152	2.0E-12	6.6E-08
cg15228509	1	243,074,717		IGR	opensea	0.840	0.555	0.292	2.4E-28	8.3E-23
cg19765154	2	191,524,409	NAB1	Body	opensea	0.715	0.885	−0.177	1.0E-25	2.4E-20
cg20262915	2	191,524,489	NAB1	Body	opensea	0.395	0.489	−0.104	4.7E-12	1.4E-07
cg06759085	2	191,555,076	NAB1	3ʹUTR	opensea	0.780	0.895	−0.109	6.3E-16	3.1E-11
cg11532947	2	191,555,480	NAB1	3ʹUTR	opensea	0.890	0.939	−0.047	1.4E-12	5.0E-08
cg00148935	3	16,398,839	RFTN1	Body	opensea	0.881	0.646	0.233	1.1E-20	9.2E-16
cg11643285	3	16,411,667	RFTN1	Body	opensea	0.954	0.837	0.117	1.7E-23	3.0E-18
cg17238319	3	16,428,391	RFTN1	Body	opensea	0.937	0.753	0.184	2.6E-22	3.1E-17
cg08034535	4	140,375,469	RAB33B	1stExon	Island	0.109	0.057	0.049	6.3E-12	1.8E-07
cg26516287	7	12,629,275	SCIN	Body	opensea	0.875	0.810	0.064	2.9E-11	7.9E-07
cg06513015	7	64,459,246	ERV3-1	5ʹUTR	opensea	0.864	0.681	0.190	8.2E-20	6.4E-15
cg07852945	9	84,303,915	TLE1	TSS1500	Island	0.230	0.080	0.144	8.4E-18	4.9E-13
cg02716779	9	131,016,143	DNM1	3ʹUTR	shelf	0.785	0.679	0.108	8.0E-16	3.7E-11
cg25294185	11	65,487,814	RNASEH2C	Body	Island	0.022	0.136	−0.116	1.6E-21	1.6E-16
cg09241427	11	70,860,363	SHANK2	5ʹUTR	opensea	0.906	0.817	0.095	1.4E-12	5.1E-08
cg09516963	12	68,042,445	DYRK2	TSS200	Island	0.676	0.317	0.377	1.2E-18	7.4E-14
cg25438338	14	24,701,525	GMPR2	TSS1500	shore	0.166	0.091	0.075	2.1E-12	6.6E-08
cg23719534	15	101,099,284		IGR	Island	0.969	0.880	0.093	6.6E-19	4.6E-14
cg04946709	16	59,789,030	LOC644649	Body	Island	0.659	0.814	−0.149	2.2E-15	9.0E-11
cg20299935	17	21,795,943		IGR	opensea	0.812	0.692	0.122	4.5E-13	1.7E-08
cg08906898	20	34,319,899	RBM39	Body	opensea	0.767	0.906	−0.130	3.9E-17	2.1E-12
cg17612569	21	27,107,221	GABPA	TSS200	Island	0.047	0.219	−0.170	7.3E-23	1.0E-17

Chr, chromosome, cgi, CpG Island, IGR, intergenic region, NNC, non-neurological control, MS, Multiple Sclerosis, Δβ, delta-beta values. Probes highlighted in bold are included in differentially methylated region.

Figure 1. DNA methylation changes in non-neuronal nuclei of Multiple Sclerosis patients. a. Circos plot illustrating differentially methylated positions (DMPs) in non-neuronal nuclei of the normal appearing white matter of Multiple Sclerosis (MS) compared to white matter of non-neurological controls (NNC). The outer track is a hg19 ideogram illustrating chromosome and cytoband information. The upper track represents the corresponding -log₁₀ (P-values), with dark red indicating sites passing significance threshold of Benjamini-Hochberg-adjusted P-value, P_adj < 0.05 and light red illustrating CpGs with P < 0.001. The inner track corresponds to the effect size of DNA methylation changes (Δβ-value) between MS and NNC, with red and blue depicting hyper- and hypomethylation, respectively. b Heatmap of the DMPs (P_adj < 0.05) with red and blue depicting z-score transformed hyper- and hypomethylation, respectively. c. Distribution of DMPs (P_adj < 0.05) across classical gene features including promoter-like features (TSS1500, TSS200, 1stExon, 5′UTR), gene body and 3′UTR (shades of green) and CpG Island (CGI)-related features including CpGs Islands, shores, shelves and open sea (shades of orange). Violin plots and barplot depict β-values (with the mean coloured in red) and DMPs distribution (in comparison to EPIC background), respectively. *P < 0.05, **P < 0.001 using Fisher’s exact test for enrichment and depletion analyses. All DMPs are listed in Supplementary Table 2.

We next explored the distribution of DNA methylation changes (DMPs, P_adj < 0.05) according to gene- and CpG Island (CGI)-related features and found overall higher β-values at DMPs in MS compared to controls across promoter-related features (TSS1500, TSS200, 5ʹUTR) while lower methylation was observed in the first exon and 3ʹUTR (Figure 1c). Fisher’s exact test further showed that DMPs were significantly enriched in distal promoter (TSS1500, P = 0.003), first exon (P = 0.008) and intergenic region (P = 0.004), while being depleted from gene bodies (P = 0.018), UTRs (5ʹUTR, P = 0.002; 3ʹUTR P = 0.031) and proximal promoter (TSS200, P = 0.024) (Figure 1c). DMPs were not significantly enriched in FANTOM 5-annotated enhancers. Enrichment analysis according to CGI-related features revealed an enrichment of DMPs in shores (P = 6.04 x 10⁻⁰⁷) and to a lesser extent in CGI (P = 0.001), along with depletion in open seas (P = 2.03 x 10⁻⁰⁵) (Figure 1c).

To further expand the informative nature of methylation changes in MS glia, we conducted transcription factors enrichment analysis by mapping JASPAR-annotated transcription factor recognition motifs to each genomic location targeted by the EPIC probes. Results from Fisher’s exact test revealed significant enrichment of two transcription factor binding sites (NR2C2 and SOX18) within the DMPs compared to the CpGs of the EPIC array after correction for multiple testing, the strongest evidence coming from NR2C2 nuclear receptor (P = 3.30 x 10⁻¹³) (Supplementary Table 4).

Altogether, these data imply noticeable DNA methylation changes in glial cells of MS patients compared to controls occurring in regulatory segments of genes.

DNA methylation changes affect genes involved in cytoskeleton, motility, signaling and metabolism

To gain insight into the biological relevance of the DNA methylation changes in glial cells of MS patients, we conducted a Gene Set Analysis of genes harbouring DMPs (P_adj < 0.05, 687 genes). Clustering of Gene Ontology (GO) terms indicated alteration of genes linked to six main functions: cytoskeleton organization, cell signaling, molecule transports, neurogenesis, cell motility and metabolic processes (Figure 2a, Supplementary Table 5). Accordingly, altered genes form a biologically interconnected gene network (interaction enrichment P-value = 1.9 x 10⁻⁰³) with the core DMP-genes involved in the GO categories (n = 254 genes) presented in Figure 2b. More specifically, as listed in Table 3, DMP-genes implicated in cellular motility encode cell-adhesion molecules, such as cadherin and integrin, together with players of chemotaxis and extracellular matrix (ECM) remodelling (e.g., collagen and hyaluronan dynamics). Cytoskeleton dynamics genes are further represented by cytoskeleton-associated proteins (e.g., RhoGTPases) as well as vesicle-mediated transport linked to endocytosis, exocytosis and intracellular vesicle trafficking. The most represented intracellular signaling pathways are connected to Wnt/β-catenin and TGF-β/SMAD signaling followed by inflammation-mediated signaling pathways. Neurogenic signaling processes are for example reflected by molecules involved in signaling through glutamate and GABA, along with players of axo-glial processes, including myelination. Metabolic processes encompass differential methylation at genes encoding proteins involved in mitochondria integrity and oxidative stress, nutrient homoeostasis, molecule degradation and ion transport (predominantly potassium channels). Transcriptional regulation comprises a plethora of DNA-binding transcription factors involved in nervous and immune cell fate, proliferation and cell arrest as well as DNA repair and chromatin regulation. Differential methylation also implicated genes involved in RNA synthesis including ribosomal and viral transcription. GO analysis of genes harbouring DMP conditioned to their genomic location showed that enrichment of these processes arise from DMPs located in most gene segments (Supplementary Fig. 4).

Table 3. Gene ontology analysis of differentially methylated genes (DMP with P_adj < 0.05) associated with Multiple Sclerosis in glial nuclei.

Category	GeneSet	Representative DMP-genes
Cellular motility and cytoskeleton dynamics	GO:0034446, GO:0006935, GO:0040011, GO:0034114, GO:0030036, GO:0097435, GO:0031346, GO:0030199, GO:1901028	Cell-adhesion and chemotaxis: CDH23, CDH4, ITGB2, ITGBL1, ITGA8, PF4, CSF3R, LTB4R2, PLGRKT, ANTXR1, MDGA1, PARVG, F11R, CGNL1, PALLD, PTPRJ
		ECM remodelling: HAS1, COL5A1, COLA3, COL11A2, LOXL3, LOXL4, PLGRKT, FAM20C, FBLN1, TMOD2
		Cytoskeleton rearrangement: ARHGEF7, PKN2, MCF2L, KANK1, TRIO, DST, SEPT9, SYNPO2, SWAP70, SORBS3, LPAR1, STMN1, EPS8
		Vesicle-mediated transport: TRAK1, CLN5, DNM1, LRP12, ATP8A2, GPRC5B, SCIN, PTPRN2, PACSIN2, RAB11B, CORO7, MYO18A
Signaling pathways and nervous processes	GO:0007267, GO:0046578, GO:0035023, GO:0071560, GO:1905114, GO:0022008, GO:0099536	Wnt signaling: WNT4, WNT6, WNT10A, WIF1, LRP5, TSPAN12, VAX2, TLE1, TCF4, CSNK1D
		TGF-β signaling: TGFBR3, SMAD2, ACVR2A, VASN
		Inflammation, apoptosis: STAT5A, MAPK8/JNK1, IL1RL2, MIF, IL17RA, AES/CD209, PTGER4, BNIP2, PPP1R13B, NTSR1
		Neurotransmission: GRIN2A, GRIN2D, SLC1A6, SLC17A7, CACNG7, GABRP, GABBR2, SLC6A13, SLC6A2, CHRNB2, HTR1A, NALCN, OPRK1, ADRA1A
		Axo-glial processes: MBP, PLLP, ARHGEF10, SHANK1, SHANK2, SHANK3, SYN2, SYN3, BEGAIN, OTOF, ICA1
		Other downstream signaling: PDE2A, GAB2, PLCG2, CIB1, MAPK7, DAB2IP
Metabolism and molecule transport	GO:0043269, GO:0030001, GO:2000191, GO:0051938, GO:0046391, GO:0045893, GO:0001678.	Oxidation and energy: NOX4, CYBA, GPX1, NAB1, GABPA, ACAA2, MUL1, IRS2, SLC5A6, KMO, FFAR2, RPTOR, DDX6, USP14
		Ion transport: KCNG2, ATP4A, KCNK9, KCNG1, KCNA1, KCNAB2, KCNN3, KCNMA1, KCNJ12, DPP6, ATP1A3, TRPM4, CACNA1E, ATP2B2, SLC9A3
		Transcriptional regulation: transcription (FOXA2, FOXJ2, GLI2, GLI3, EGR2, SOX7, BCL11B, ETS1, AES, HIVEP3), proliferation/cell arrest (CDK13, TFDP1, E2F4, MNAT1, ZNF16, PLAGL1, FBXO31, APBB1), chromatin regulation (KDM3A, PRMT5, PRDM16, CITED2, CREBBP, ZXDC)
		RNA synthesis: TAF2, PRPSAP1, PRPSAP2, HNRNPH1, WDR43, UBTF, NFIA, BRF1, ERV3-1

Figure 2. Functional annotation of DNA methylation changes associated with Multiple Sclerosis. a. Scattered plot showing clusters of gene ontology (GO) terms associated with differentially methylated genes (DMP with P_adj < 0.05) in glial nuclei of MS patients compared to NNC. Biological processes GO terms were obtained using over-representation analysis and clusters were visualized in two-dimensional space (assigning x and y coordinates to each term) by applying multidimensional scaling of the matrix of GO terms according to pairwise semantic similarity, using REVIGO [54]. The circle size represents −log10 (P-value), with small to big diameters ranging from 2.07 (P = 8 × 10⁻⁰³) to 5.44 (P = 3 × 10⁻⁰⁶). Top GO terms associated to different categories of Biological processes are presented in the right panel. b. Representation of the genes containing DMPs (P_adj < 0.05) involved in the biological processes using STRING network analysis. Colours represent different clusters (kmeans clustering set at 5 clusters). Grey gradient indicated the strength of data support (darker grey representing stronger evidence, dotted line showing lower level of evidence, i.e., combined interaction score 0.4–0.6). Full GO data are presented in Supplementary Table 5.

Functional annotation of DMRs confirmed GO analysis of DMPs with enrichment of ECM and migration (e.g., ITGB2, MANBA, NEU4, SCIN, STMN1), TGF-β and Wnt signaling pathways (GREM2, SMAD2, TLE1) and metabolism (NOX5, CRYZ, CYP1A1, LDHAL6A) among others (Supplementary Table 5). Examples of DMR genes implicated in these processes are illustrated in Figure 3. They encode cell surface integrin subunit (ITGB2) involved in adhesion, migration and phagocytosis, among others, mitochondrial cytochrome P450-mediated enzymes (GSTM5), developmental transcription factor antagonizing Wnt signaling (VAX2), key players in involved in neurogenesis and iron homoeostasis (SEZ6L2 and BOLA2, respectively), insertase to the ER membrane (WRB), mediator of necroptosis (MLKL) and IFN-induced dynamin-like GTPase (MX2).

Figure 3. Differentially methylated regions in non-neuronal nuclei of Multiple Sclerosis patients. Plots illustrating genes associated with differentially methylated regions (DMRs) identified in glial nuclei from Multiple Sclerosis (MS) cases compared to non-neurological controls (NNC). UCSC genome browser annotations are shown: the hg19 ideogram illustrating chromosome and cytoband information, the complete gene structure of the locus (blue track), CpG Island (CGI) when present (green track), DNAseI hypersensitivity and transcription factor (TF) binding clusters (grey gradient) as well as other regulatory properties exemplified in cell lines from chromatin state segmentation by hidden Markov model from ENCODE/Broad with the following colour coding for active promoter (dark red), weak promoter (light red), poised promoter (purple), strong/weak enhancer (orange/yellow), insulator (blue), strong/weak transcription (dark/light green) and strong/weak heterochromatin (dark/light grey). The DMR (black) location appears in the last track. Methylation (β-values) of single consecutive CpGs within the DMR for cases and controls is depicted in red and blue, respectively, with connecting lines indicating mean methylation ± SEM. All DMRs are listed in Supplementary Table 3.

DNA methylation changes are associated with gene expression differences in MS glial cells

We sought to further explore the putative functional impact of the identified methylation changes by examining expression of the corresponding genes from published bulk and single-cell RNA-sequencing data conducted in post-mortem brain MS and control samples. Comparison of differentially methylated genes with transcripts detected in RNA-seq data of bulk NAWM of MS versus WM of NNC individuals [34] showed that a minor fraction (~14%, 80/561) displayed significant transcriptional differences, as reported in the original analysis with P < 0.05 (Figure 4, Supplementary Table 2). The majority of them (51/80) were found downregulated in MS NAWM compared to control WM, this association being higher than expected for genes containing DMPs within gene body (Chi-square test P = 0.019) (Figure 4). Overall, the transcriptionally dysregulated DMP-genes encode proteins primarily involved in adhesion and migration of nervous cells (e.g., SLIT3, DAB2IP, GLI3, NRP1, CDH4 genes) as well as cytoskeleton remodelling and vesicle trafficking (e.g., ATP8A2, TMOD2, OBSCN, TRAK1, KIF5C, PARVG genes) and inflammatory response (e.g., DAB2IP, ITGB2, CYBA, LILRA4 genes).

Figure 4. Association of DNA methylation changes with gene expression in the normal appearing white matter of Multiple Sclerosis patients. Scatterplot illustrating association of DMPs (P_adj < 0.05) methylation values (Δβ) in glial nuclei samples with gene expression data (RNA-seq) reported in bulk NAWM of Multiple Sclerosis (MS) patients compared to WM of non-neurological controls (NNC) [34]. Red and blue colours indicate upregulation and downregulation, respectively, in MS versus NNC bulk brain tissue. The barplot represents the proportions (percentage) of upregulated (red) and downregulated (blue) genes for all gene-annotated DMPs and DMPs located in promoter (TSS1500, TSS200) or gene body. The dotted line indicating the expected proportion * P < 0.05 generated with the Chi-square test. logFC, log2-fold change, Δβ, difference in the beta-value.

To gain insight into cell type-specific dysregulation in MS, we annotated DMP-genes according to single-cell RNA-seq analyses of MS lesions compared to NNC brain tissue [13,14,52] (Supplementary Table 2). Results indicated that 16% (115/687) of DMP-genes were found differentially expressed in glial cell types of MS lesions compared to NNC brain tissue (Supplementary Table 2). Of note, some of the DMP-genes identified in MS lesions at single-cell resolution were also found differentially expressed in bulk NAWM tissue compared to WM of NNC individuals. These DMP-genes reported as differentially expressed both in the bulk NAWM and in glial cells from MS lesions reflect alterations in oligodendrocytes (DYSF, RBFOX1, PPP1R12B, KIF5C, PER3 genes), microglia (GPX1, HIVEP3, CDCP1, HERPUD1, HLA-DPA1, FYB, SCIN genes), astrocytes (PRDM16, GLI3, WIF1, SLIT3 genes) or in multiple glial cell types (KCNMA1, NTM, NRP1, CDH4, TRIO genes) in MS compared to NNC individuals (Supplementary Table 2).

These findings jointly suggest that methylation changes identified in glial nuclei of MS patients could, at least partly, be associated with transcriptional differences.

DNA methylation changes affect shared and distinct pathways among glial cells

DNA methylation changes detected in glial cells likely reflect molecular changes occurring in several glial cell types. To further disentangle their possible contribution in MS brain, we assigned these alterations to specific glial cell types based on their constitutive expression in the healthy human brain [50,51] (Supplementary Fig. 5, Supplementary Table 6). We explored the functional implication of epigenetic dysregulation at DMP-genes (P_adj < 0.05) assigned to astrocytes (n = 309 genes), microglia (n = 153 genes), and oligodendroyctes (n = 196 genes) using GO analysis. Overall, GO terms clustered into distinct groups showing both common and distinct functions between glial cell types (Figure 5a, Supplementary Table 5). MS-associated enrichment of biological functions related to cytoskeleton dynamics and ECM organization was shared by all glial cell types (Figure 5b). Both microglia and astrocytes manifested enrichment of terms linked to adhesion and differentiation while formation of cellular projection and neurogenesis could be associated to oligodendrocytes and astrocytes (Figure 5b).

Figure 5. Cell type-specific functional association of DNA methylation changes in glial cells of Multiple Sclerosis patients. a. Heatmap of significant GO terms related to Biological Processes (min 3 molecules) generated from genes that associate with cell type-annotated DMPs. GO terms were grouped into clusters using the relative risk (RR) between different pathways, which was calculated based on the number of overlapping genes per pathway [55]. The left panel reflects relative distance between GO terms and clusters. b. Significantly enriched biological processes shared by at least 2 out of 3 cell types. c. Scattered plot showing clusters of gene ontology (GO) terms associated with differentially methylated genes (DMP with P_adj < 0.05) in MS patients compared to NNC, annotated as expressed in astrocytes (left), microglial cells (middle) and oligodendrocyte (right). Biological Processes GO terms were obtained using over-representation analysis and clusters were visualized in two-dimensional space (assigning x and y coordinates to each term) by applying multidimensional scaling of the matrix of GO terms according to semantic similarity, using REVIGO [54]. The circle size represents −log10 (P-value), with small to big diameters ranging from 1.32 (P = 0.04) to 4.2 (P= 6.5 × 10⁻⁰⁵). The top biological processes are displayed in the lower panel. d. Representation of the cell type-restricted genes containing DMPs (P_adj < 0.05) involved in the top 50 biological processes from GO analysis in astrocytes (green), microglia (purple) and oligodendrocytes (orange), using STRING network analysis. Grey gradient indicated the strength of data support (darker grey representing stronger evidence, dotted line showing lower level of evidence, i.e., combined interaction score 0.4–0.6). Full GO data are presented in Supplementary Table 5. Astrocytes, microglial cells and oligodendrocyte are depicted in green, purple and orange, respectively.

In addition to these shared cellular functions, distinct signatures could be observed for each cell type-assigned DMP-gene sets (Figure 5c). Astrocyte-annotated DMP-genes were predominantly implicated in intracellular signaling such as TGF-β (TGFBR3, SMAD2, ITGA8, VASN) and Wnt (GLI3, WIF1, SOX, TSPAN12, AES, TLE1) signaling pathways, neurotransmission (reflected by several neurotransmitter receptors), along with mitochondria damage (ACAA2, MUL1, MAPK8/JNK1, PPP1R13B), DNA repair (FAM168A, BPHL, APBB1, MGMT, RTEL1) and regulation of transcription (represented by multiple transcriptional regulators). Microglial DMP-genes showed an overrepresentation of innate immune-related functions such as inflammatory response (IL1RL2, FFAR2, SLC15A4, MAPK8/JNK1, PTGER4, SWAP70, IL17RA, PLCG2), defence response (LILRA4, FCGRT, MX1/2, RNASE6, ALPK1, CDK13), oxidative stress (CYBA, GPX1), cell motility (CDH4, CSF3R, ITGB2, HAS1, COL5A1, USP14) and signaling (CITED2, CSNK1D). The analysis also underscored genes implicated in vesicle-mediated transport involved in endocytosis (CLTA, FCGRT, SCIN, PACSIN2, USP6NL, RIN2, RUFY1), ER-Golgi trafficking and autophagy (BET1, COPB1, RAB33B) as well as co-and cargo-transport (CNST, SLC15A4). Functions enriched in DMP-genes assigned to oligodendrocytes were connected to synaptic transmission (PXK, USP14, ADRA1A, SYN3, GRIN2A/2D, CACNA1E), potassium and sodium transport (DPP6, KCNG1, KCNA1, NALCN, FAM155A, PDE2A), cytoskeletal rearrangement (STMN1, LPAR1, RAP2A, ARHGEF7, KANK1, KIF5C), signaling adaptor protein (GAB2, APBB1) and myelination (ARHGEF10, MBP, NAB1, PLLP). Examples of gene networks associated to GO terms in each cell type are illustrated in Figure 5d. GO analysis of the cell type-assigned DMP-genes according to their gene location indicated that the majority of the enriched biological processes reflect methylation changes occurring in the gene body and further uncovered specific functions that are more enriched in promoter-occupying DMPs (Supplementary Fig. 6, Supplementary Table 5).

Reactome clustering of biological interactions confirmed the GO findings and further delineated specific pathways within each biological process (Supplementary Fig. 7, Supplementary Table 5). Accordingly, Wnt signaling found enriched in astrocyte-assigned genes arise primarily from the Repression of Wnt target genes set. Terms related to neurotransmission characterizing astrocyte- and oligodendrocyte-assigned DMP-genes were involved, among others, in protein-protein interaction at synapses pathways such as Neurexins and Neurogilins. Transcriptional processes found enriched in DMPs-gene assigned to all cell types referred to FOXO-mediated transcription, while astrocytic and microglial DMP-genes were additionally involved in Transcriptional regulation by AP-2 and TP53.

These findings collectively support alterations of shared processes prevailing in all cell types, such as cytoskeleton organization, as well as putative distinct cell type-specific functions converging to Wnt/TGF-β signaling in astrocytes, cellular motility and innate immunity in microglia and ion transport and neuromodulation in oligodendrocytes.

Discussion

We exploited the stable nature of DNA methylation, informing on the genome activity in post-mortem material, to investigate the still elusive molecular changes occurring in NAWM non-neuronal cells of MS patients in comparison to WM cells of NNC individuals. We found DNA methylation changes at genes involved in cellular motility, cytoskeleton rearrangement, cell-to-cell and intracellular signaling, such as Wnt and TGF-β signaling, neuromodulation and neuroinflammation, among others, a fraction of these genes were previously shown to be dysregulated in the brain of MS patients. Our findings strongly suggest that NAWM glial cells in MS are highly altered in the absence of lesional insult, collectively exhibiting a multicellular response to diffuse inflammation. Whether the observed changes reflect pre-lesional processes or tissue reaction to adjoining inflammatory damage remains to be elucidated.

Glial cells of MS patients exhibited epigenetic alterations affecting several genes implicated in cellular and intracellular motility, from substrate adhesion molecules and cell–cell junction to cytoskeleton dynamics and ECM remodelling. Differences in cellular migration and cytoskeleton dynamics have been described in developing or healthy CNS tissue and in the context of focal insult, such as MS demyelinating plaques [56,57]. The identification of such alterations in our cohort, composed of NAWM tissue, supports the occurrence of abnormalities in the unaffected areas as well, as previously suggested in NAWM bulk tissue and neurons [34,35] and might reflect tapering gradient of reactive gliosis or focal event prior to injury. In line with this, active microglia assume enhanced agility and rapidly converge to sites of damage, directional migration following a chemoattractant gradient of soluble molecules released by damaged cells and intertwined astrocytes [10,58]. Accordingly, dysfunctional astrocytes present with retracted processes and loss of cell–cell junctions, and ablation of reactive astrocytes hinders the recruitment of microglia to demyelinating lesions [24,59]. Such inhibition has been shown to impair debris clearance, oligodendrocyte maturation and remyelination [60]. The latter process also mobilizes intense cytoskeleton dynamics involved in the formation of growth cones and amyloid-like myelin-enriched protrusions ensheathing and enwrapping denuded axons [61].

Alteration of genes connected to cell-to-cell and intracellular signaling through Wnt and TGF-β families further underscore a coordinated response in the NAWM of MS patients due to to diffuse damage and confirmed previous observation of aberrant Wnt and TGF-β activity in the MS brain [27,62,63]. A pivotal role of TGF signaling in experimental autoimmune, neuroinflammatory or demyelinating diseases is already highly recognized [64]. Notably, the differentially methylated genes of Wnt/β-catenin pathways operate at multiple levels of the signaling cascade, from complex formation between extracellular Wnt ligand (WNT4, WNT6, WNT10A), Wnt ligand antagonist (WIF1) and receptors (LRP5), to subsequent β-catenin stabilization (TSPAN12) and activation of TCF/LEF-mediated transcriptional programme (TLE1, TCF4, VAX2). Dysfunctional Wnt pathways in MS-like animal models has further shown to enhance astrogliosis, hinder remyelination, and impair cell migration to the demyelinated lesions [62,65–67].

Interestingly, glial cells of the MS-NAWM displayed epigenetic alterations of molecules involved in neuronal integrity and plasticity, compared to WM glial cells of controls. Methylation changes affected genes encoding synapse scaffolding and clustering proteins, such as postsynaptic density proteins (PSD), synapsins and neurexin, and can be exemplified by the robust hypermethylation (DMP with Δβ > 0.15), spanning over a DMR at exon 2 of the PSD95-associated BEGAIN gene. Hypermethylation of the very same region has been identified in bulk brain tissue and glia of patients affected by neuropsychiatric disorders and inversely correlated with expression changes [68]. More generally, genetic aberrations of many of the differentially methylated genes in our study, i.e., TRIO, NRXN2, SYN2, CACNA1E, DNM1 and RAB11B, have been associated with movement disorders, cognitive and visual impairments and brain abnormalities (e.g., atrophy, white matter apoplasia) [69]. Methylation changes at genes encoding modulators of synaptic transmission, in particular potassium/calcium transport and glutamate transmission further reinforces the role of glial cells in supporting neuronal functions in general. Accordingly, the glial control of neuronal activity strongly relies on extracellular potassium control, i.e., both astrocytes and oligodendrocytes respond to neuronal discharge and reciprocally coordinate neuronal excitability by buffering/dispersing extracellular potassium [70]. They operate via potassium channels (illustrated by differentially methylated genes encoding voltage-gated, calcium-activated and inwardly rectifying K⁺ channels in our cohort) and astrocyte-oligodendrocyte junctional coupling embedded in a panglial syncytium [71–73]. Additionally, astrocytes govern neuromodulatory networks and promote excitatory signaling via sheathing of a vast number of tripartite synapses, expression of neuromodulatory receptors, glutamate uptake and secretion of gliotransmitters [74,75]. While this finding reflects the dynamic features of glial cells of the NAWM, the circumstances underlying such alterations in the seemingly unaffected tissue, i.e., as a cause or consequence of brain disconnectivity, warrant further investigation.

While this study is the first to examine methylation changes in non-neuronal NAWM cells, the small sample size of the cohort and the difference in potential confounders, such as age, sex and genetic background, between MS and controls represent limitations. Moreover, because epigenetic marks are highly tissue- and cell type-specific, especially in the CNS [76,77], the analysis of the NeuN-negative fraction sorted from brain tissue could undeniably be biased by cellular heterogeneity arising predominantly from mixed glial cell types, with an additional minor contribution of non-glial (e.g., endothelial cells, peripheral immune) cell types. To mitigate the impact of such potential confounders, we analysed sorted nuclei from the WM exclusively and further applied reference-free deconvolution accounting for varying proportions of cell types. However, one cannot exclude the possibility of missed signals and, whether the identified changes reflect large differences in one cell type and/or shared variations in several cell types remain to be elucidated. The poor overlap between the DMPs identified in the present study and the DMR-CpGs reported in the bulk NAWM [34] is likely due to the noticeable disparities (e.g., normalization, batch effect correction, statistical testing) in the analytical strategies. Overall, our approach likely captures the molecular signature of interwoven processes, as reported at the transcriptomic level [27] and as suggested in our study by the strong enrichment of DMPs with sequence binding to the ubiquitous and pleiotropic NR2C2 transcription factor. The nuclear hormone receptor NR2C2 regulates pivotal processes as various as neuronal and astrocytic development [78], myelination [79], lipid metabolism [80,81] or DNA damage response involved in telomeres maintenance and telomere-driven genome instability [82,83]. We attempted to delineate common and cell type-specific process by integrating single-cell transcriptomic data and subsequently derive biological interpretations, which certainly warrants validation. Cell type-specific functional annotation unambiguously attributed abnormalities of cytoskeleton dynamics and ECM organization to all glial cell types. GO analyses implied a pervasive role of oligodendrocytes and astrocytes in supporting neuronal integrity via cell-cell signaling and neuromodulation, among others and unveiled more specific functions in microglia, converging to cell migration, vesicle-mediated transport, RNA metabolism and neuroinflammation such as antiviral processes. This can be exemplified by large (Δβ> 0.15) and wide (DMR) hypermethylation at the ERV3-1 gene, encoding human endogenous retrovirus (HERV)-R element. HERVs have gained large interest in neurodegenerative diseases in general and in MS in particular, notably as disease markers and putative targets for MS therapy [84–86], with a still unclear role of HERV-R [87,88]. Whether changes at the ERV1-3 locus in glial cells of MS patients reflect the physiological role of HERV-R (per se or via regulation of neighbouring genes) as reported [89] or rather a pathogenic event, as described for other HERVs remain to be ascertained.

In conclusion, DNA methylation analysis of non-neuronal glial nuclei sorted from post-mortem WM of MS patients and controls unravelled molecular changes affecting motile, signaling and neuromodulatory properties of NAWM glial cells. These alterations likely compile multiple discrete focal insults impinging the CNS in absence/prior to lesion or alternatively reflect compensatory mechanisms circumventing the gradual and global deterioration of the neural circuitry. Overall, our findings portray the NAWM as a highly dynamic tissue in response to the persistent inflammation endured by the CNS in MS patients. The identification of epigenetic abnormalities, which are modifiable by nature, affecting theses processes provides additional evidence in favour of therapeutic strategies aiming at rehabilitating the functional capacities of the CNS in progressive MS patients by targeting CNS resident cells.

Contributors

MJ, LK and MN conceived and designed the study. MJ supervised the study. LK and RC performed sample isolation and FACS sorting. MPK conducted DNA extraction. EE carried out DNA methylation analysis and MN and DGC contributed to the analytical pipeline with input. LK conducted additional analyses. LB aided in accessing the samples. LK wrote the manuscript with assistance from all authors. All authors read and approved the manuscript.

Data availability

The EPIC data that support the findings of this study are available in Gene Expression Omnibus (GEO) database under the accession number GSE166207.

Acknowledgments

We are grateful to A. van Vollenhoven for flow cytometry processing (Center for Molecular Medicine core facility). We acknowledge GenomeScan/ServiceXS (Leiden, The Netherlands) for processing Illumina EPIC arrays. 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).

Disclosure statement

No potential conflict of interest was reported by the author(s).

Supplementary material

Supplemental data for this article can be accessed here

Additional information

Funding

This study 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, the Stockholm County Council - ALF project, the European Union’s Horizon 2020 research and innovation programme (grant agreement No 733161) and the European Research Council (ERC) (grant agreement No 818170), the Knut and Alice Wallenberg Foundation grant, Åke Wilberg Foundation, Hedlund Foundation, Norlins Foundation, Bergvalls Foundation and Karolinska Institute’s funds. LK is supported by a fellowship from the Margaretha af Ugglas Foundation. MPK is supported by McDonald Fellowship from Multiple Sclerosis International Federation (MSIF). The funders of the study had no role in study design, sample acquisition, data collection, data analysis, data interpretation or writing of the manuscript.