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.


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][4][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-brainbarrier 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][12][13][14][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][17][18][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][21][22][23][24][25][26][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 postmortem 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.

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 nonneurological 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.

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 � Ω À 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 ) 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 Pvalue 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.

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://mostafa vilab.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 singlecell 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].

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 freshfrozen 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).
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 −13 ) (Supplementary Table 4). The upper track represents the corresponding -log 10 (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.
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 −03 ) 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  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 −03 ) to 5.44 (P = 3 × 10 −06 ). 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.
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). 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).

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 singlecell RNA-sequencing data conducted in postmortem 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 (Chisquare 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).
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  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).
In addition to these shared cellular functions, distinct signatures could be observed for each cell typeassigned DMP-gene sets (Figure 5c). Astrocyteannotated DMP-genes were predominantly implicated in intracellular signaling such as TGF-β (TGFBR3, SMAD2, ITGA8, VASN) Table 3. Gene ontology analysis of differentially methylated genes (DMP with P adj < 0.05) associated with Multiple Sclerosis in glial nuclei.
Reactome clustering of biological interactions confirmed the GO findings and further delineated specific pathways within each biological process ( Supplementary Fig. 7, Supplementary Table 5).  [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. 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 typeannotated 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 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 DMPgenes were involved, among others, in proteinprotein 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 DMPgenes 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 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 −05 ). 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. 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][66][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, calciumactivated and inwardly rectifying K + channels in our cohort) and astrocyte-oligodendrocyte junctional coupling embedded in a panglial syncytium [71][72][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 telomeredriven 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, vesiclemediated 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][85][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.