DNA methylation signatures of Multiple Sclerosis occur independently of genetic risk and are primarily attributed to B cells and monocytes

doi:10.21203/rs.3.rs-2570500/v1

Download PDF

Research Article

DNA methylation signatures of Multiple Sclerosis occur independently of genetic risk and are primarily attributed to B cells and monocytes

https://doi.org/10.21203/rs.3.rs-2570500/v1

This work is licensed under a CC BY 4.0 License

Version 1

posted

You are reading this latest preprint version

Multiple sclerosis is a complex autoimmune disease that causes neuronal demyelination and debilitating physical and cognitive symptoms. Epigenetic factors can mediate genetic and environmental effects on disease risk. Here we profiled blood-based DNA methylation in 583 MS cases and 643 healthy controls representing 3 independent study groups. An epigenome-wide association study was performed, incorporating statistical deconvolution to identify immune cell-specific signals, which were validated in purified cell subsets. First, we elucidated the mediation effects of methylation at the HLA risk locus. Second, we identified methylation differences in MS that occur independently of known genetic risk loci and show that these differences more strongly differentiate disease than HLA genotype and the polygenic risk score. Finally, we show that methylation differences in MS occur predominantly in B cells and monocytes and indicate the involvement of cell-specific biological pathways. Overall, this study comprehensively characterizes the immune cell specific epigenetic architecture of MS.

Multiple Sclerosis

Epigenetics

Methylation

Epigenome-Wide Association Studies

Genetic Risk

Cell Deconvolution

We identified a distinct methylation pattern in whole blood of people with MS driven partly by genotype but also by environmental exposure, these changes are mainly present in B cells and monocytes.

Multiple Sclerosis (MS) is a chronic autoimmune and neurodegenerative disease of the central nervous system (CNS). It is caused by the combination of a complex genetic predisposition and environmental and lifestyle factors. Genome-wide association studies (GWAS) have identified over 200 single nucleotide polymorphisms (SNPs) associated with MS risk, the majority of which localize to genes relevant to immune function^1,2. The strongest risk locus is found within the major histocompatibility complex (MHC) region of chromosome 6 – the well characterized HLA-DRB1*1501 MS risk allele^3–5.

Despite this, the precise pathogenesis of MS remains unclear. There is good evidence for dysregulation of immune cell balance in the CNS, particularly with regards to the immune phenotypes^6,7. Studies suggest that altered ratios between the adaptive and innate immune system cells lead to disease activity^6,7. There are also a number of well-established environmental and lifestyle factors that contribute to disease risk such as, viral infections (specifically Epstein Barr Virus (EBV), infectious mononucleosis^8–10) and smoking¹¹, all of which interact with the genome. Other factors that have no direct genomic influence are considered to contribute to disease and include low ultraviolet radiation exposure, low vitamin D levels¹², and nulliparity¹³. The exact molecular mechanisms by which multiple sclerosis develops remain elusive, but epigenetic change is emerging as a contender that significantly influences disease risk.

One of the most well studied epigenetic modification is DNA methylation (DNAm), which refers to the addition of a methyl group onto the cytosine of a cytosine-guanine (CpG) dinucleotide. Most studies have focused on CpGs that are associated with gene expression, where increased levels of CpG methylation are generally associated with gene silencing and de-methylation with gene activation. Epigenetic modulation via CpG methylation is unique in reflecting environmental exposure that can be influenced by genetic, environmental and lifestyle factors, making it a plausible genomic modifier involved in MS onset and progression.

Epigenome-wide association studies (EWAS) of blood cell DNA have identified changes in DNA methylation profiles associated with MS^14–16. Although these studies have somewhat conflicting results, they all consistently identify differential methylation at HLA-DRB1*1501, which appears to partly mediate the relationship between the risk allele at this locus and MS¹⁶. Limited isolated cell type analyses suggests that each immune cell population has a unique epigenetic profile^17–20.

To date, the isolated cell type EWAS studies have involved modest sample sizes, whereas larger whole blood EWAS studies have likely missed cell type specific differential methylation due to the masking of signal from the mixed cell populations. Statistical methods now exist to estimate the relative proportions of major blood cell types using methylation data derived from mixed cell DNA samples (such as EpiDISH, reviewed in²¹). Methods to identify cell-specific differential methylation (csDNAm) was derived from these algorithms (CellDMC).

Herein we have conducted an EWAS of whole blood DNA from MS cases and controls and applied statistical deconvolution to reveal cell-specific methylation signatures associated with disease. Importantly, the discovery group (from the Ausimmune study) included cases with early MS (principally first clinical diagnosis of CNS demyelination) that have had a follow-up assessment of disease after 10 years²². This allowed us to identify causative epigenetic changes in specific cell subsets that occur prior to, or early in, disease onset rather than being a consequence of longstanding pathology or therapeutic intervention. Replication analysis confirmed our results in two independent studies.

Dna Methylation Analysis

To identify DNA methylation changes associated with MS, we performed EWAS on bisulfite treated genomic DNA isolated from whole blood (see Fig. 1 and methods section). A total of 583 people with MS (pwMS) and 643 controls from three independently collected groups were analyzed in this study. Further details of each group is in Supplementary table 1 and the methods section.

The discovery group was comprised of whole blood samples collected at baseline from The Australian Multi-Centre Study of Environment and Immune function (Ausimmune study)²². From this group, we used 208 MS cases and 402 matched non-MS controls. Most of this group entered the Ausimmune study after their first clinical diagnosis of CNS demyelination (FCD), were primarily treatment naïve at the time of sampling, and had mainly relapse onset phenotypes. We used whole blood samples from the baseline collection in cases who, at 10-years follow up, had converted to MS based on the 2017 McDonald criteria²³. Controls were matched for age, sex and region²².

The first replication group was from Sweden and consisted of 140 MS cases and 139 controls. This dataset consisted of previously published data from prevalent MS cases who were primarily on treatment at the time of blood sampling, mostly relapse onset phenotype, and had been profiled using the Illumina 450K arrays¹⁶.

The second replication group consisted of 235 MS cases sourced from four Australian MS specialist clinics: Melbourne, Newcastle, Hobart and Adelaide. The cases in this group were exclusively of relapse onset, but longer disease duration (median 16 years) and were on a variety of treatments at sample collection²⁴. This case group did not have controls so were matched with data from 102 non-MS controls from a publicly available dataset.

All DNA samples from the discovery groups and case samples from the second replication group were tested in-house for genome-wide methylation using the Illumina Infinium Human Methylation EPIC BeadChip (EPIC array) (Illumina Inc, SanDiego, CA, USA). The publicly available data set matched to the second replication group was also profiled using the EPIC array.

Whole blood EWAS identified major differential methylation at HLA genes

To identify DNAm changes in whole blood associated with MS, we first identified differentially methylated positions (DMPs) i.e., DNAm changes at single CpG sites. We used logistic regression with phenotype (case or control) as the outcome variable and incorporated major covariates in the model, such as age, sex, region, and cell-type proportion (see methods section). The mean difference in CpG methylation between the case and control group was used to assess the effect size (referred to as Δβ).

In the discovery group, we identified 3,218 DMPs, which met the genome-wide significance cut-off of P < 9.8x10^− 8. Of these, 159 DMPs also met a delta beta cut off of Δβ ≥ 0.02 (an effect size of 2% or greater). To confirm these results, we investigated if these DMPs were also present in our two replication groups using an empirical p value of < 0.05 and a Δβ > 0.02. We found 18 DMPs that overlapped in both replication groups with the same directionality as the discovery group. The majority of these DMPs (13/18) localised to HLA-DRB1, HLA-DRB5 or HLA-DRB6 genes (Table 1). Of the remaining 5 DMPs, 4 localized to genes related to interferon pathways and one DMP localized to an unannotated location (supplementary Table 2). Because there were 32 individuals in the discovery group who were treated with beta interferon prior to sample collection, we determined if the DMP signal at these four genes was originating from interferon beta treated individuals. To do this, we removed the 32 samples from the analysis and repeated it. After removal, their effect size reduced and two of the four DMPs no longer reached genome-wide statistical significance (located on CMPK2 and PARP9, see Supplementary Table 2).

Table 1

**The 13 DMPs at the HLA-DRB locus in all three groups.** All DMPs with Δβ > 2% and genome-wide significant in the discovery group as well as FDR ≤ 0.05 and Δβ > 2% in both replication cohorts.
CpG	Chr	Position	Gene	Feature	Discovery Cohort Δβ	Discovery Cohort p-value	Replication Cohort 1 Δβ	Replication Cohort 1 p-value	Replication Cohort 2 Δβ	Replication Cohort 2 p-value
cg01341801	6	32489203	HLA-DRB5	Body	0.192	3.43E-11	0.212	6.51E-07	0.276	1.94E-05
cg25140213	6	32522683	HLA-DRB6	Body	0.187	3.14E-10	0.189	1.63E-06	0.238	4.55E-05
cg17369694	6	32485396	HLA-DRB5	3'UTR	0.155	2.65E-11	0.174	1.90E-06	0.217	1.46E-04
cg09139047	6	32552042	HLA-DRB1	Body	-0.135	1.36E-08	-0.162	2.70E-06	-0.175	2.70E-03
cg14645244	6	32552205	HLA-DRB1	Body	-0.122	7.18E-10	-0.131	6.30E-07	-0.128	2.46E-04
cg10632894	6	32552453	HLA-DRB1	Body	-0.12	6.55E-08	-0.127	2.01E-06	-0.207	5.35E-05
cg08845336	6	32551891	HLA-DRB1	Body	-0.119	4.01E-09	-0.14	2.16E-06	-0.137	4.94E-05
cg10104420	6	32526414	HLA-DRB6	Body	-0.107	1.96E-09	-0.135	5.60E-07	-0.071	3.11E-02
cg19575208	6	32551888	HLA-DRB1	Body	-0.105	3.55E-09	-0.121	3.26E-06	-0.128	9.10E-05
cg17416722	6	32554385	HLA-DRB1	Body	0.096	4.67E-11	0.105	1.88E-06	0.14	2.72E-04
cg00440797	6	32493873	HLA-DRB5	Body	0.079	5.40E-10	0.07	5.95E-05	0.124	1.24E-04
cg00579921	6	32489991	HLA-DRB5	Body	0.064	7.06E-08	0.051	6.31E-04	0.101	6.14E-03
cg26036029	6	32552443	HLA-DRB1	Body	-0.061	4.72E-08	-0.042	5.30E-05	-0.076	8.69E-05
cg22930808	3	12228188	PARP9	5'UTR	-0.043	6.66E-08	-0.079	8.45E-08	-0.088	5.90E-04
cg01028142	2	7004578	CMPK2	Body	-0.036	6.28E-12	-0.039	4.77E-06	-0.062	1.57E-05
cg01533966	1	90363165	LRRC8D	5'UTR	0.03	1.41E-12	0.019	3.94E-02	0.03	2.47E-02
cg01311537	10	50386709	TMEM273	Body	-0.021	4.98E-08	-0.017	7.78E-03	-0.035	4.08E-03
cg10778971	14	94577101	IFI27	5'UTR	-0.02	4.43E-12	-0.014	1.15E-03	-0.014	9.12E-03
Δβ = delta beta

To increase statistical power, we combined all three datasets (see methods section). The combined dataset analysis identified 11,969 individual DMPs that met the genome-wide significance threshold of (P < 9.8x10^− 8). Of these, 190 also had a Δβ > 0.02 and 81 were located at the MHC Class II locus (Supplementary Fig. 1). Of the total 11,969 individual DMPs, 7803 (65%) were hypermethylated, 9234 (77%) mapped to the immediate vicinity of genes, and 2116 (18%) were located at the transcription start site of genes. These 11,969 DMPs represented a total of 5,560 genes.

From the DMP list of the pooled dataset, we identified differentially methylated regions (DMRs) using the DMRcate package by defining a DMR as three or more DMPs separated by less than 1000bp. To ensure maximum biological relevance we filtered the list to contain only DMRs where at least one DMP had a Δβ greater than 0.1 (10%). After this filtering, nine independent DMRs remained, all of which were located in HLA genes (Table 2). They contained between five and 40 DMPs with effect sizes ranging from 1–18% (Δβ 0.01–0.18). There were four DMRs in HLA-DRB1 (three hypermethylated and one hypomethylated), two in HLA-DRB6 (both hypermethylated), and one each in HLA-DRB5, HLA-DQA (both hypermethylated) and HLA-DQB (hypomethylated).

Table 2

DMRs in the combined cohort.
Chromosome	Start	End	Width	Number of CpGs	Maximum Difference	Mean Difference	Overlapping Genes	DMR Index
chr6	32548321	32550067	1747	10	0.161	0.051	HLA-DRB1	1
chr6	32551749	32552670	922	19	-0.167	-0.077	HLA-DRB1	2
chr6	32553920	32556093	2174	5	0.109	0.059	HLA-DRB1	3
chr6	32557422	32558459	1038	5	0.126	0.029	HLA-DRB1	4
chr6	32489203	32490444	1242	12	0.182	0.056	HLA-DRB5	5
chr6	32520615	32523136	2522	11	0.152	0.064	HLA-DRB6	6
chr6	32525805	32526702	898	12	0.117	0.034	HLA-DRB6	7
chr6	32604564	32610971	6408	24	0.113	0.043	HLA-DQA1	8
chr6	32632000	32636189	4190	40	-0.111	-0.012	HLA-DQB1	9

DNA methylation mediates genetic risk of HLA-DRB1

Given the strong genetic associations at HLA-DRB1 with MS we investigated if the four DMRs in this gene are influencing the relationship between genotype and phenotype (i.e., mediating the genetic effect through epigenetic alteration). Genotyping array data for cases and controls was available for both the discovery group and the first replication dataset to test this hypothesis. For each of the four DMRs we selected an index CpG (i.e., yielding the largest delta-beta) to serve as a proxy for the respective DMR. We first checked the methylation quantitative effect loci (mQTL) effect of HLA-DRB1 and found that regardless of status (case or control) the results reveal a striking relationship between methylation and risk haplotype at this locus (Supplementary Fig. 2).

We then performed mediation analysis with the causal inference test (CIT) method ²⁵ to test for causality between HLA-DRB1 risk haplotype (casual factor), DMR methylation (mediator) and CDMS as the phenotype. CIT analysis allows the establishment of a chain of causality between genotype, methylation and MS, ensuring that the genotype effect on MS is mediated through methylation. We used rs3135388 as the tagging SNP for the main HLA-DRB1 haplotype. Results of the CIT showed that the methylation differences identified between pwMS and controls are driven entirely by genotype in both the discovery and the first replication group (P < 0.05) (Fig. 2). We did not have the same genotype information from all cohorts for the HLA-DRB5 haplotype; however, in the discovery group, the CIT HLA-DRB5 reached statistical significance using the tagging SNP rs73405151 (data not shown).

Methylation differences act independently of genetic risk

Given the strong correlation between methylation and genotype at the MHC region, we wanted to assess epigenetic risk acting independently from known genetic risk of MS. To do this we computed EWAS models controlling for the effect of both HLA-DRB1 risk haplotype and non-MHC SNPs identified in the 2019 IMSGC GWAS², using a polygenic risk score (PRS). This was only performed in the discovery group given a) there was missing genetic data from other groups and b) the discovery group represents early-stage MS, which is an important attribute for risk assessment. From this analysis, we identified 1627 DMPs reaching statistical significance (P < 9.8x10^− 8), 123 of which had a Δβ above 2% (See supplementary File 1). Interestingly, none of the DMPs localized to the HLA-DRB1 gene, further indicating that the epigenetic effects we had identified at this locus are linked to genetic variation. In order to assess if identified DMPs were modulated by genotype or most likely the result of lifestyle and/or environmental exposure, we tested each CpG for cis-acting methQTLs. Only 104 out of 1627 (6.4%) produced statistic evidence of modulation by a nearby genetic locus (see Fig. 3d). Over representation analysis (ORA) on the list of 1092 unique genes identified after genotype correction revealed emerging pathways in MS such as NOTCH signaling pathways and axon guidance pathways (see Fig. 3c).

Using these DMPs we then constructed a genetic risk corrected methylation risk score (grcMRS) for each subject. This grcMRS represents epigenetic risk of MS independent of known genetic risk in this group. We then compared the grcMRS to HLA-DRB1 risk haplotype and the PRS to assess how well each index distinguishes MS from healthy controls. First, we show that PRS and HLA-DRB1 risk haplotype are not correlated, indicating independence of these genetic risk factors. We subsequently constructed a total genetic risk score, incorporating both the PRS and HLA-DRB1 risk haplotype. As expected, the grcMRS correlated very weakly with total genetic risk (see Fig. 3a).

To assess how those scores can act as classifiers for MS, we calculated AUC-ROC curves (see Fig. 3b). We show that the total genetic risk score (HLA-DRB1 + PRS) provides moderate discriminatory accuracy (AUC = 0.72, 95% CI: 0.66–0.76). By comparison the grcMRS performed better at distinguishing MS from healthy controls (AUC = 0.85, 95% CI: 0.82–0.89). Combining both genotype and methylation scores only marginally improved the accuracy (AUC = 0.87, 95% CI: 0.84–0.91) as too did the addition of other covariates such as cell proportions, sex and age (AUC = 0.89, 95% CI: 0.86–0.92).

Cell Specific Differential Methylation Is Mainly Attributed To B Cells And Monocytes

To identify if specific immune cell types are responsible for the differential methylation identified in the whole blood samples, we employed a two-step approach using the combined dataset. First, we calculated cell proportion estimates in the combined dataset using the R package EpiDISH (24). This revealed a statistically significant decrease in the proportion of Natural Killer (NK) cells (mean_case = 0.021, mean_control = 0.031, P = 1.34e^− 08) and CD8⁺ T cell proportions (mean_case = 0.095, mean_control = 0.106, P = 5.92e^− 05) and a statistically significant increase in CD4⁺ T cell proportions (mean_case= 0.081, mean_control = 0.071, P = 1.58e^− 03) in MS cases compared to controls. There were no differences in the proportion of monocytes, B cells and neutrophils between cases and controls (supplemental Fig. 3).

Second, we used these cell type proportions, in conjunction with CpG methylation (M) values, to calculate cell-type specific DMPs (csDMPs). Here we employed an adaptation of the functions in CellDMC ²⁶. Although it accounts for all cell types simultaneously, the base CellDMC regression model can become overburdened when incorporating many cell-types (i.e., many terms), thus reducing power and preventing the identification of cell specific effects. To overcome this, we constructed a "per-cell type" linear model (see methods section).

Genome-wide analysis of cell specific methylation patterns revealed an over-representation of csDMPs in both B cells and monocytes. From highest number of csDMPs to lowest were B cells > monocytes > CD4⁺ T cells > CD8⁺ T cells > NK cells then neutrophils (Fig. 4a). Interestingly, csDMPs were predominantly hypermethylated in B cells (over 60%) while monocytes had a more evenly distributed set of csDMPs (Fig. 4a). Consistent with whole blood analysis, there were marked differential methylation effects at HLA, but notably, this effect exhibited obvious variation across cell subtypes (Fig. 3b). More detailed analysis of the cell-specific modelling revealed that the strong differential signal observed at the HLA region in whole blood originates primarily from monocytes and B cells but was also present to a lesser extent in T lymphocytes (Fig. 4). It was notably absent or nearly absent in both NK cells and neutrophils (Supplementary file 2).

Comparison of the four cell types with highest number of csDMPs, that met both the significance and effect size cut offs, (monocytes, B cells, CD4⁺ T cells and CD8⁺ T cells) revealed that both B cells and monocytes had largely a different set of csDMPs. There were 1165 CpGs unique to B cells, 771 CpGs unique to monocytes and 347 common csDMPs between the two. When we investigated all the cell types we found that CD4 + T cells, CD8 + T cells, monocytes and B cells shared 28 common csDMPs, with the same directionality between them located on 7 different genes: HLA-DQB1, CMPK2, IGF2R, HLA-DRB1, HLA-DRB6, HLA-DRB5 and HLA-DQA1 (Supplementary Fig. 4 and supplementary table 3).

Our replication groups had an older mean age and a longer disease duration than the discovery group. To examine if the predominant methylation signal coming from B cells and monocytes occurs early in disease pathology or is a result of long-standing disease, we performed the same analysis described above in only the discovery group. We found the majority of DMPs appeared to be from B cells and monocytes, with absent or nearly absent signals from both NK cells and neutrophils (supplementary Table 3). Similarly, there was differential methylation predicted to a lesser extent in T lymphocytes.

DNA methylation at HLA-DRB1 correlates with gene expression in B cells and monocytes

We did not have accompanying RNA for the whole blood groups used in this study to examine gene expression. Therefore, to determine if methylation was related to gene expression, we examined the relationship between methylation and gene expression using RNA extracted from isolated monocytes, CD19⁺ B cells and CD4⁺ T cells from our previously published study¹⁶, one unpublished study, and unpublished dataset (Maltby et al, unpublished) (see methods for group details). All three were cross-sectional, unrelated MS case-control studies.

We used Pearson’s correlation to estimate the relationship between methylation and expression levels of HLA-DRB1 in all three cell types (Fig. 5) and found a strong negative correlation between gene expression and methylation at DMR-2 located in exon 2 of HLA-DRB1. This data indicates hypermethylation at this locus is associated with decreased gene expression in these immune cell subtypes (Fig. 5). As expected, based on the DNA methylation levels, there was a strong correlation between HLA-DRB1 expression and methylation levels, but this correlation was not as strong in the CD4⁺ T cells.

Cell types are listed in rows, CpGs are in columns, except for column 1, which represents gene-of-interest expression in each cell type. The correlation coefficient (R) is represented through the circle size and colour. Larger circle area indicates high absolute R value (maximum 1). Red indicates negative values; blue indicates positive values. The black rectangle represents the hypomethylated DMR 2 at HLA-DRB1 identified in whole blood and Cell specific analysis

Over-representation analysis shows enrichment of DMPs at immunological pathways

In the combined whole blood group, we identified 190 DMPs with Δβ > 2%. These DMPs mapped to 92 unique genes (Supplementary file 3). Over-representation analysis (ORA) using this gene list identified 28 pathways. The ORA showed an enrichment in immunological pathways, including cytokine and interferon-related pathways as well as antigen presentation through class I and II MHC (supplementary Fig. 6).

Our cell specific analysis identified 1740 csDMPs in B cells and 1212 csDMPs in monocytes (Fig. 3a). We created a gene list for B cells and monocytes which contained genes associated with at least one csDMP in either cell type. We then used these resulting gene lists in combination with the whole blood gene list to perform over-representation analysis. All the pathways identified in B cells were also present in both whole blood and monocytes. However, the gene list from monocytes generated unique pathways. For example, cytokine signaling and IFN alpha/beta pathways identified in whole blood were only found in monocytes (Fig. 6).

Despite the significant progress in developing treatments for MS, the pathogenic molecular events that occur early in the disease process are not understood. The combined effects of genetics and environmental risk factors suggest that epigenetic mechanisms are likely to be involved. Numerous studies have now been published which investigate the genome-wide DNA methylation changes in whole blood and isolated cell types; however, results are inconsistent and limited by heterogeneous cell populations, disease heterogeneity, use of a wide variety of disease modifying therapies (DMTs) or small sample size ^{14,16,18,19,27}. Using a large dataset of MS cases early in their diagnosis, with age/sex/location matched controls, we have confirmed previous findings that DNA methylation acts as mediator of the HLA-DRB1*1501 allele ¹⁶ and have provided evidence that this occurs early in disease pathology rather than as a result of treatment use or long-standing disease. Furthermore, we have revealed that differential methylation is primarily associated with B cells and monocytes, with a smaller contribution from T lymphocytes (CD4⁺ and CD8⁺). This is consistent with our previous observations in these four isolated cell types ¹⁸. We have also correlated methylation change with gene expression at the HLA-DRB1 locus thereby substantiating its key role in MS pathogenesis, which we replicated in two independent patient cohorts. However, outside of the MHC region, significant methylation changes appear to be unrelated to meaningful genotype changes thought to modulate MS risk (identified through the 2019 IMSGC GWAS).

The lack of discriminatory (risk) accuracy emerging from both PRS and HLA-DRB1 haplotype underlines the limited role played by genetics in MS risk². Correcting for known strong MS risk loci identified robust CpGs that are largely not modulated by genotype. We provide evidence that methylation alone is a stronger discriminator of MS compared to previously described genetic risk factors. A chain of causation still needs to be established, identifying if these changes in methylation are causative or a consequence of MS. Since these findings occur in an early stage MS patient group (most at FDE) we can rule out that these effects are a consequence of chronic disease. As DNA methylation is substantially influenced by environmental factors this result highlights the importance of environmental exposures in MS risk. Further studies need to systematically evaluate the effect of environmental exposure toward methylation change and MS risk such as seen in vitamin D²⁸, smoking²⁹ or EBV. Interestingly, the controlling for the influence of genotype also unveiled non-immune-related pathways disrupted in pwMS such as NOTCH signaling or axon guidance pathways, both emerging as key pathways in MS^30,31.

A major strength of our study is the discovery group, collected early in diagnosis and followed prospectively over a ten-year period. This allowed us to examine the methylation profiles of patients at baseline collection, who, at their 10-year follow up visit had fulfilled the criteria for MS. We have demonstrated that B cells and monocytes account for the largest contribution of whole blood DNA methylation signal in MS and that this is not simply a result of long-standing disease. This is consistent with previous work on isolated cell subtypes from people with MS ^17,18,32 and with current knowledge of B cells playing an important role in the pathogenesis of disease, highlighted by the effectiveness of treatments targeting the CD20 surface antigen predominantly expressed on B cells ³³. There has been the suggestion that these therapies specifically target the memory B cell compartment ³⁴ so future studies investigating the methylation of memory B cells will be of interest. In contrast to previous work, the csDMPs identified in B cells were predominantly hypermethylated.

Intriguingly, monocytes, which also have antigen presenting properties, likely play an important role in early stage disease by mediating both pro- and anti-inflammatory responses ^35–37. Furthermore, our ORA and GAT analysis revealed that monocyte differential methylation signals included unique pathways. This suggests monocytes may play an early role in MS pathogenesis.

Increased proportions of CD4⁺ T cells ^38,39, and deficit of CD8⁺ T cells ⁴⁰ in MS cases compared to controls, as seen in our estimations, has been previously documented. The increased proportions of CD4⁺ T cells is likely driven by Th17 helper cells, whereas the deficit of CD8⁺ T cells is probably due to reduced effector cell numbers ^38–40. Our cell proportion estimates are consistent with these past findings. The proportions of total NK cells has not previously been evaluated in MS; however, we have previously demonstrated subsets of NK cells are affected by disease modulating therapies and an increase of NK_bright cells is associated with stable disease ⁴¹.

To overcome issues of heterogeneous cell populations and small sample size, we used an adapted CellDMC method to identify which cell types contribute the most to the differential methylation signal. We were able to demonstrate that the main source of differential methylation in whole blood comes from B cells and monocytes, particularly the HLA-DR region. This is not surprising, given that HLA-DR is expressed primarily on antigen presenting cells such as B lymphocytes and monocytes ⁴². We also found a strong differential methylation signal from both CD4⁺ and CD8⁺ T cells. HLA-DR is normally absent on resting T lymphocytes; however, it has been proposed as a marker of immune system activation ⁴³. In systemic lupus erythematosus, CD8⁺ T and CD4⁺ T cells have increased HLA-DR expression during active disease leading to speculation that expression of HLA-DR represents an inflammatory response. It is plausible this is also happening in T lymphocytes of people with MS, but further studies are required to confirm this.

In addition to the 13 DMPs located at HLA genes, we also identified two robust signals that localized to interferon related genes. This suggests that these DMPs may be partially influenced by treatment. One small study found global methylation changes in monocytes after treatment with interferon beta ⁴⁴. The authors of this study, however, were unable to identify changes at specific CpGs ⁴⁴. A study of monozygotic twin pairs with MS found DMPs associated with INF-beta treatment however, we did not find any overlap with this study ⁴⁵. A larger sample size and prospective study design to confirm these findings is required.

The results of our csDMP analysis was compared to two previously reported studies, Ewing et al. ⁴⁶ and Ma et al. ³², which revealed 154 shared genes containing csDMPs in all three studies in B cells and two shared genes in monocytes (Supplementary Fig. 7), suggesting these may be very robust DMPs for these cell types. No shared DMPs were identified in either CD4⁺ or CD8⁺ T cell populations. Differences would be expected between studies due to differing assays (450K arrays (Ewing) and bisulfite sequencing (Ma)), different group statistical power and different group characteristics. The overlap between them likely represents replicated MS-associated DMPs.

One limitation of this study is that the combined dataset is heterogenous in terms of clinical characteristics. Although our discovery group was primarily treatment naïve and early in their disease course, the replication groups contained participants on mixed treatments and with long standing disease. Thus, we cannot draw conclusions regarding how treatment use or disease progression affects methylation change. Furthermore, despite the large size of the combined dataset, the statistical deconvolution methods are still only likely to be identifying major effect DMPs. Smaller effect CpGs probably exist but much larger studies are needed to identify them. Another limitation is that the expression data is from separate participant samples; therefore, we cannot assign causation between expression and methylation, we can only assign correlations. Future studies with larger groups of treatment naïve populations, followed prospectively over time, are necessary to confirm our DNA methylation and expression data to confirm the results presented herein.

This study provides compelling evidence that DNA methylation, independent of genotype, is strongly modified in pwMS, implicating the NOTCH and axon guidance pathways. The results of this study advance previous work implicating B-cells and monocytes in MS pathogenesis by revealing differential methylation is present at disease onset. The effect of differential methylation of B cells and monocytes is linked to gene expression differences in these cell types. While therapies effecting B cells have been shown to be extremely effective in MS, monocytes may also represent a future therapeutic target. Finally, we show that methylation, after correction for known genetic risk, outperforms genetics when discriminating between MS and non-MS individuals. This provides compelling evidence that environmental exposure and lifestyle play a larger role in MS onset than genetics.

See Fig. 1 for a study workflow diagram

Patients’ enrolment and blood collection

Sample Groups

Clinical characteristics have been listed in supplementary table 1.

Whole Blood methylation:

Discovery Group: details of the Ausimmune Study design have been previously published ²². Briefly, between 2003 and 2006, participants were recruited across four different centres – Brisbane city, Newcastle city and surrounds, Geelong city, the Western Districts of Victoria and the Island of Tasmania. Cases were aged 18–59 and entered the study around the time of their first demyelinating event (FDE). They have since been prospectively followed over 10 years. The Australian electoral roll was used to select controls from the general population at baseline. Non-MS controls were matched based on age (± 2 years), gender and region with a ratio of up to 4 controls for each case.

Epidemiological information was collected by face-to-face interview, self-completed questionnaires, and neurological examination. At the same time, biological samples were collected and stored. Detailed analysis of the information collected has been previously published ⁴⁷. For this study, we had access to whole blood from 208 MS cases and 402 matched controls. All 208 MS cases had converted to clinically definite MS (CDMS) by the 10 year follow up.

The first replication group was from a previously published case-control study ¹⁶ (group 2), and 450K array data was available from Gene Expression Omnibus (GEO) database under accession number GSE43976 and GSE106648. The dataset is comprised of whole blood DNA from 140 MS cases and 139 non-MS controls. Samples were collected in Sweden between 2005 and 2011, the cases were mostly on low-efficacy treatment (65%), and mostly relapsing-remitting phenotype.

The second replication group was comprised of 235 MS cases obtained from the specimen banks of 4 Australian specialist centres (Melbourne, Newcastle, Hobart and Adelaide) and 102 publicly available control idat files obtained through EWAS data hub (https://ngdc.cncb.ac.cn/ewas/datahub/repository) using the following criteria: Age between 23 and 73, sex = Female, platform = 850K, tissue = whole blood and sample type = control. The second replication group is publicly available on the GEO database (GSE106648) and was previously used in whole blood case-control EWAS (Campagna et al, submitted).

DNA extraction and Quality Control

For the discovery group and second replication dataset, DNA was extracted from whole blood using the QIAamp DNA Blood Mini Kit ™, Netherlands. Extracted DNA was quantified using the Qbit (Invitrogen™, USA) and integrity was assessed using the genomic DNA tapes on a TapeStation (Agilent™, USA) using the DNA Integrity Number (DIN) as a metric. All samples had a DIN ≥ 7, with minimal genomic DNA degradation.

Methylation and genotyping arrays

To determine SNP profiles, Global Screening Arrays 2.1 (Illumina™, USA) (hereafter referred to as GSA arrays) were used according to the manufacturers’ guidelines. For genome-wide methylation analysis, 500ng of genomic DNA was bisulfite-converted using the EZ-DNA Methylation™ Kit (Zymo) according to manufacturer’s converted DNA was hybridized to the Illumina Methylation EPIC BeadChip arrays (hereafter referred to as EPIC arrays). To avoid batch effects, samples were randomised on both GSA and EPIC arrays using the OSAT R package ⁴⁸. Both GSA and EPIC arrays were read using the iScan (Illumina™) to produce raw idat files.

Methylation analysis

For the discovery and both replication groups, analysis was performed using the ChAMP R package ^49,50. To summarise, idat files were loaded and filtered to remove low performing probes, probes mapping to multiples loci and low performing samples. Probes next to known polymorphisms and on XY chromosome were kept. Beta values were then normalised using the BMIQ method ⁵¹. Batch effects on both array and chip levels were corrected using the Combat algorithm ⁵².

The final discovery dataset contained 742,961 CpGs and 610 individuals (208 cases and 402 controls). The first replication dataset contained 478,687 CpGs across 279 samples. The second replication dataset contained 831,858 probes across 337 samples.

The final model used to identify differentially methylated positions (DMPs) was constructed using logistic regression, whereby phenotype (Case or Control) was considered the outcome variable and each CpG beta value and cell-type proportions (for natural killer cells, monocytes, B cells, CD8⁺ T cells, CD4⁺ T cells and neutrophils) were considered predictors. For CpG I:

$$glm({Outcome}_{\frac{Case}{Control}} ∼ {CpG}_{i}+Age+Sex+CellFraction+ Region )$$

The mean difference in methylation between case and control group was used to assess the effect size (ie. delta beta or Δβ).

When performing the genotype risk corrected EWAS (correcting for known genotype risk loci (HLA-DRB1 RISK HAPLOTYPE and aggregated polygenic risk score), we also incorporated genotype informations in the model to identify genotype risk independent DMPs:

$$glm({Outcome}_{\frac{Case}{Control}}∼ {CpG}_{i}+Age+Sex+CellFraction+ Region+HLA Haplotype+PRS )$$

DMRs were considered significant with Fisher multiple tests, Stouffer and the harmonic mean of the individual component FDRs p ≤ 0.05.

Sensitivity Analysis

Sensitivity analysis was performed on a range of covariates to assess their relative individual importance on CpG beta-value. Numeric continuous variables were assessed using linear models such as $lm({CpG}_{i}∼Covariate)$ and categorical variables using ANOVA models such as $AOV({CpG}_{i}∼Covariate)$.

Genotype Analysis

Raw idat files were loaded in GenomeStudio (Illumina) for filtering. Recommendations from Guo et al. 2014 ⁵⁴ were followed (except the ethnicity-specific steps) and the analysis ready dataset was exported to plink ⁵⁵ format. In short, QC was performed to filter haploid probes, then probes were filtered based on their GenTrain score (score ≥ 0.7), their cluster separation, SNP call rate (≥ 98%) and sample call frequencies (≥ 95%). After exporting to plink format, samples were tested for gender mismatch, relatedness, and Hardy-Weinberg equilibrium. Plink was also used to generate association statistics, using Fisher’s exact test. The HLA-DRB1*15:01 haplotype was detected using the rs3135388 tag SNP ⁵⁶.

For the CIT analysis, we used Haploview to identify the distribution of haplotypes. For each, gene and group, we identified the haplotype that is most significant in a Χ² test as a risk haplotype.

Scores and receiver operator characteristic curves

We constructed a genotype risk corrected Methylation Risk Score (grcMRS) only incorporating DMPs identified in the genotype risk corrected EWAS. The grcMRS was constructed with the following formula:

$${\text{M}\text{R}\text{S}}_{sample}= \sum {\beta value}_{ at CpGi} \times {\varDelta \beta }_{CpGi}$$

Where ${\beta value}_{ at CpGi}$ represents the beta-value of this particular sample at a cpg of interest and ${\varDelta \beta }_{CpGi}$ is a weight, using the value of delta-beta identified for this in the genotype risk independent EWAS (described in the methylation analysis section) .

For genotyping, we used two scores. First, the number of risk haplotype HLA-DRB1 RISK HAPLOTYPE (either 0, 1 or 2) and a polygenic risk score (PRS) generated using PRSice ⁵⁷ with all SNPs indetified throught the 2019 IMSGC GWAS².

ROC curves and area under the curve (AUC) were computed using the R package ROCR⁵⁸ using the scores generated with a logistic regression.

$$glm({Outcome}_{\frac{Case}{Control}} ∼ {Predictor}_{1}+\dots + {Predictor}_{n})$$

Where Outcome is the MS-Control outcome and Predictors are MRSs, PRS, cell proportion, sex, etc.

Immune cell deconvolution analysis

Cell proportion estimates were obtained using the R package EpiDISH ⁵⁹. Cell proportions, in conjunction with M-values, were used to calculate cell-type specific DMPs, using a method inspired by the CellDMC function ²⁶. CellDMC had two major limitations. First, results were biased based on initial cell proportion when using beta-values. Second, CellDMC modelling is overburdened when incorporating numerous cell-types, preventing the identification of effect happening in all or multiple cell types.

To compensate for this, we used a two-step approach using the combined dataset for maximum statistical power. First, we calculated cell proportion estimates using the R package EpiDISH. Second, we used the cell type proportions, in conjunction with CpG methylation (M) values, to calculate cell-type specific DMPs (csDMPs). Here we employed an adaptation of the functions in CellDMC ²⁶. Although it accounts for all cell types simultaneously the base CellDMC regression model can become overburdened when incorporating many cell-types (i.e. many terms), thus reducing power and preventing the identification of cell specific effects happening in multiple cell types.

Instead, for each cell-type, we constructed a linear model as follow:

$$lm({CpG}_{i} ∼ cellFrac+cellFrac:Phenotype+Covariates)$$

Where CpG_i represent the M-value at CpG_i, cellFrac is the cell-type fraction (ranging from 0 to 1) and cellFrac : Phenotype represents the interaction term between cell fraction and phenotype. CpG with an interaction $cellFrac:Phenotype$ estimate above 2 or under − 2 as well as a p-value under 9.8x10^− 8 were kept for further downstream analysis.

Expression and DNA methylation Analysis

Monocytes

Expression data was comprised of data from independent, cross sectional, case-control study. Expression data from monocytes from previously published work was acquired from the authors under a data transfer agreement and DNA methylation data 450K arrays was freely available from GEO ¹⁶. Monocytes were isolated from PBMCs by magnetic separation (Miltenyi) according to the manufactures’ instructions and expression of HLA-DRB1 was quantified by real-time PCR as previously described ¹⁶.

CD4⁺ T cells

Our group has previously published the results of the DNA methylation analysis using 450K arrays. From this study, we also had matched RNA available from 18 of the relapse-onset MS cases (unpublished results). The MS cases were all female, who were either treatment naïve or had been off MS-specific treatment for at least 6 months. They were matched with 33 age and sex matched non-MS controls.

The RNA from the CD4⁺ T cells was extracted using the RNAeasy mini kit (Qiagen, USA) following manufacturers’ protocols. The resulting RNA quality was assessed using the Agilent Bioanalyzer and met a quality cut off of a RIN number of 7 or higher. Samples meeting QC requirements were hybridized to the Illumina HT12 expression array cards (service provided by Diamantina Institute, Australia).

Expression analysis was carried using the R package beadarray ⁶⁰. In short, idat files were loaded using the readIdatFiles function, data was normalised using the normaliseIllumina with the quantile method. Finally, the limmaDE function was used to determine gene differential expression and we used the illuminaHumanv4.db package to annotate the results.

B cells

The B cell dataset was composed of a total of 73 samples, 23 controls (13 non-inflammatory neurologic disease controls and 10 healthy individuals) as well as 50 MS cases (39 RRMS and 11 SPMS).Peripheral blood mononuclear cells (PBMCs) were isolated directly after collection using a standard Ficoll (GE Healthcare) and sodium citrate-containing preparation tubes (Becton Dickinson) procedures, respectively. CD19⁺ B cells were isolated using positive selection on MACS microbeads magnetic separation (Miltenyi), according to manufacturer’s instructions (> 95% purity). Cells where prepared within an hour after sampling using AutoMACS (Milteny Biotec, Bergisch, Germany) according to the manufactures standard protocol and stored at − 70°C. DNA and RNA was extracted simultaneously using Qiagen Allprep DNA/RNA kit (Qiagen, Venlo, Netherlands) according to the manufactures standard protocol and stored at − 70°C. The amount and quality of DNA were assessed with a NanoDrop ND-1000 Spectrophotometer (NanoDrop Technologies Inc). The samples with sufficient DNA amounts were used in further analysis. Processing of the samples for Infinium HumanMethylationEPIC arrays (Illumina), including bisulfite conversion, was done at the National Genomics Institutet (NGI) at Science for life lab (Uppsala, Sweden).

RNA-Seq libraries were generated from 500 ng of total RNA using the Illumina TruSeq mRNA Stranded Library Preparation Kit (cat. no. RS-122-2103) according to the manufacturer's protocol. Library quality was determined using the Agilent High Sensitivity DNA Kit (cat. no. 5067 − 4626) and NanoDrop ND-1000 Spectrophotometer (NanoDrop Technologies Inc). Libraries were sequenced on the Illumina HiSeq 2500 as per the manufacturer's instructions, and ∼20M 75 bp paired-end reads were generated per sample. Read quality was assessed before and after trimming using FastQC v0.11.4. Illumina adapters and low-quality nucleotides were trimmed using Cutadapt v1.9.1. Reads were aligned to the human genome (GRCh38), and read count per gene was determined using STAR.

All cell types:

Methylation profiles for every cell type were analyzed individually in R using the Minfi ⁶¹ and ChAMP package ⁶² following the pipeline according to Marabita el al. ⁶³. Briefly, type 1 and type 2 probes were normalized using quantile normalization and BMIQ. Sex of the samples was confirmed using the GetSex function from the Minfi package and the cell type identity was confirmed using the cell type deconvolution method from Minfi based on the Houseman algorithm ⁶⁴. Age was confirmed using the AgeP function from the WateRmelon package⁶⁵. Furthermore, sample identification was confirmed by running a 44 SNP panel at NGI, and genotypes were confirmed by overlapping SNP calls with data from our genetic studies². The following probes were filtered out: i) probes not passing the detection p-value cutoff of 0.01, ii) probes with known SNPs, and iii) X and Y chromosome probes, except for the data relating to pwMS using rituximab and DMF, because of the pairwise analysis. Batch effects were identified using principal component analysis (PCA) and corrected using ComBat from the SVA package ⁶⁶. The loading of methylation profiles was performed in this manner for each dataset used in this study. DMPs were determined with linear modelling using the limma package ⁶⁷ in a model that included age, sex and B cell purity from the Houseman algorithm as covariates.

Correlation analysis (DNA methylation and gene expression)

Correlation analysis was performed using the R package corrplot. All three cell types were analysed separately. In every case, a measure of expression (from RT-PCR, array or sequencing) was correlated to a beta-value for every sample. Pearson R coefficient is reported on the plots.

metQTL analysis

We tested for the influence of genotype on methylation in two instances. Each time, we performed an ANOVA test:

$$AOV(CpG ∼ {SNP}_{AA-AB-BB})$$

Where CpG represents the beta value at a specific CpG and ${SNP}_{AA-AB-BB}$ represents the genotype groups.

We tested both the influence if the HLA-DRB1*15:01 risk haplotype on genome wide methylation and tested all SNPs within ± 18kBp (36kBp window) for each DMPs identified in the genotype risk independent EWAS.

Over-Representation Analysis

ORA was performed using the online Reactome platform⁶⁸ to analyse various gene lists with the inclusion of IntAct⁶⁹ interactors into the analysis. Pathways were kept if the Entities FDR was under 0.05. Figures were generated using the R package ggplot2⁷⁰.

DMP	Differentially Methylated Position
DMR	Differentially Methylated Region
CIT	Causal Inference Test
MS	Multiple Sclerosis
CNS	Central Nervous System
MHC	Major Histocompatibility Complex
SNP	Single Nucleotide Polymorphism
EBV	Epstein-Barr Virus
DNAm	DNA Methylation
CpG	5'—C—phosphate—G—3'
EWAS	Epigenome-Wide Association Study
csDNAm	cell specific DNA methylation
pwMS	people with Multiple Sclerosis
FCD	First Clinical Diagnosis of CNS demyelination
mQTL	methylation Quantitative Trait Loci
PRS	Polygenic Risk Score
ORA	Over-Representation Analysis
MRS	Methylation Risk Score
grcMRS	genetic risk corrected Methylation Risk Score
AUC	Area Under the Curve
csDMP	cell specific Differentially Methylated Position
DMT	Disease Modifying Therapies
GEO	Gene Expression Omnibus

Acknowledgments: AusImmune Investigator team: Robyn Lucas, Ingrid Van der Mei, Patricia Valery, Keith Dear and Caron Chapman. The authors would also like to acknowledge Peri Noori and David Gomez-Cabrero for their contribution to the work on isolated B cells. Finally, the authors would like to thank all the pwMS and control subjects who participated in this research.

Funding: This study was supported by the following funding bodies: NHMRC (AusLong) grant APP1127819, NMSS grant (AusImmune) RG-1803-30499 and MSRA (severity study) 18-0424

Author contributions: RJS, RAL and JLS contributed to the conception of the work. RJS, RAL and JLS contributed to the design of the work. AX, JNT, MS, HB, IK, LK, VGJ, TK, MJ, LA, BT, RJS, ALP, RAL and JLS contributed to the acquisition of data. AX, VEM, EE, MPC, SMB and RAL contributed to the analysis of data. AX, VEM, RAL and JLS contributed to the interpretation of data. AX, VEM, RJS, RAL and JLS contributed to the manuscript drafting. VEM, HB, VGJ, TK, MJ, BT, RJS, ALP, RAL and JLS contributed to the manuscript substantial revision.

Competing interests: HB has received institutional (Monash University) funding from Biogen, F. Hoffmann-La Roche Ltd, Merck, Alexion, CSL, and Novartis; has carried out contracted research for Novartis, Merck, F. Hoffmann-La Roche Ltd and Biogen; has taken part in speakers’ bureaus for Biogen, Genzyme, UCB, Novartis, F. Hoffmann-La Roche Ltd and Merck; has received personal compensation from Oxford Health Policy Forum for the Brain. JLS received travel compensation from Biogen, Merck and Novartis; has been involved in clinical trials with Biogen, Merck, Novartis and Roche; her institution has received honoraria for talks and advisory board service from Biogen, Merck, Novartis and Roche. VEM has accepted honoraria for presentations and research funds from Biogen and Merck.

Data and materials availability: The datasets generated and analysed during the current study are not publicly available due to the General Data Protection Regulation but are available from the corresponding author on reasonable request.

Consent for publication: Not applicable

Ethics approval and consent to participate: All participants to this study gave their written informed consent. The Ausimmune Study was approved by nine regional Human Research Ethics Committees (study number 09/04/15/5.04 and HREC/09/HNE/139).

International Multiple Sclerosis Genetics, C. et al. Analysis of immune-related loci identifies 48 new susceptibility variants for multiple sclerosis. Nat Genet 45, 1353-1360, doi:10.1038/ng.2770 (2013).
International Multiple Sclerosis Genetics, C. Multiple sclerosis genomic map implicates peripheral immune cells and microglia in susceptibility. Science 365, doi:10.1126/science.aav7188 (2019).
Kellar-Wood, H. F. et al. Multiple sclerosis and the HLA-D region: linkage and association studies. Journal of neuroimmunology 58, 183-190 (1995).
Sawcer, S. et al. A high-density screen for linkage in multiple sclerosis. American journal of human genetics 77, 454-467, doi:10.1086/444547 (2005).
Patsopoulos, N. A. et al. Fine-mapping the genetic association of the major histocompatibility complex in multiple sclerosis: HLA and non-HLA effects. PLoS Genet 9, e1003926, doi:10.1371/journal.pgen.1003926 (2013).
Attfield, K. E., Jensen, L. T., Kaufmann, M., Friese, M. A. & Fugger, L. The immunology of multiple sclerosis. Nat Rev Immunol, doi:10.1038/s41577-022-00718-z (2022).
Bar-Or, A. & Li, R. Cellular immunology of relapsing multiple sclerosis: interactions, checks, and balances. Lancet Neurol 20, 470-483, doi:10.1016/S1474-4422(21)00063-6 (2021).
Ruprecht, K. The role of Epstein-Barr virus in the etiology of multiple sclerosis: a current review. Expert Rev Clin Immunol, doi:10.1080/1744666X.2021.1847642 (2020).
Dunn, N., Kharlamova, N. & Fogdell-Hahn, A. The role of herpesvirus 6A and 6B in multiple sclerosis and epilepsy. Scand J Immunol, e12984, doi:10.1111/sji.12984 (2020).
Bjornevik, K. et al. Longitudinal analysis reveals high prevalence of Epstein-Barr virus associated with multiple sclerosis. Science 375, 296-301, doi:10.1126/science.abj8222 (2022).
Hedstrom, A. K. Smoking and disability progression in multiple sclerosis. Expert Rev Neurother 20, 739-741, doi:10.1080/14737175.2020.1743176 (2020).
Lucas, R. M. et al. Sun exposure and vitamin D are independent risk factors for CNS demyelination. Neurology 76, 540-548, doi:10.1212/WNL.0b013e31820af93d (2011).
Ponsonby, A. L. et al. Offspring number, pregnancy, and risk of a first clinical demyelinating event: the AusImmune Study. Neurology 78, 867-874, doi:10.1212/WNL.0b013e31824c4648 (2012).
Maltby, V. E. et al. Genome-wide DNA methylation profiling of CD8+ T cells shows a distinct epigenetic signature to CD4+ T cells in multiple sclerosis patients. Clin Epigenetics 7, 118, doi:10.1186/s13148-015-0152-7 (2015).
Maltby, V. E. et al. Differential methylation at MHC in CD4(+) T cells is associated with multiple sclerosis independently of HLA-DRB1. Clinical Epigenetics 9, doi:ARTN 71 10.1186/s13148-017-0371-1 (2017).
Kular, L. et al. DNA methylation as a mediator of HLA-DRB1(star)15:01 and a protective variant in multiple sclerosis. Nat Commun 9, doi:ARTN 239710.1038/s41467-018-04732-5 (2018).
Maltby, V. E. et al. Genome-wide DNA methylation changes in CD19(+) B cells from relapsing-remitting multiple sclerosis patients. Sci Rep-Uk 8, doi:ARTN 17418 10.1038/s41598-018-35603-0 (2018).
Ewing, E. et al. Combining evidence from four immune cell types identifies DNA methylation patterns that implicate functionally distinct pathways during Multiple Sclerosis progression. EBioMedicine 43, 411-423, doi:10.1016/j.ebiom.2019.04.042 (2019).
Bos, S. D. et al. Genome-wide DNA methylation profiles indicate CD8+ T cell hypermethylation in multiple sclerosis. PLoS One 10, e0117403, doi:10.1371/journal.pone.0117403 (2015).
Rhead, B. et al. Increased DNA methylation of SLFN12 in CD4+ and CD8+ T cells from multiple sclerosis patients. PLoS One 13, e0206511, doi:10.1371/journal.pone.0206511 (2018).
Campagna, M. P. et al. Epigenome-wide association studies: current knowledge, strategies and recommendations. Clin Epigenetics 13, 214, doi:10.1186/s13148-021-01200-8 (2021).
Lucas, R. et al. Observational analytic studies in multiple sclerosis: controlling bias through study design and conduct. The Australian Multicentre Study of Environment and Immune Function. Mult Scler 13, 827-839, doi:10.1177/1352458507077174 (2007).
Thompson, A. J. et al. Diagnosis of multiple sclerosis: 2017 revisions of the McDonald criteria. Lancet Neurol, doi:10.1016/S1474-4422(17)30470-2 (2017).
Campagna, M. P. et al. Whole-blood methylation signatures are associated with and accurately classify multiple sclerosis disease severity. Clin Epigenetics 14, 194, doi:10.1186/s13148-022-01397-2 (2022).
Millstein, J., Zhang, B., Zhu, J. & Schadt, E. E. Disentangling molecular relationships with a causal inference test. BMC Genet 10, 23, doi:10.1186/1471-2156-10-23 (2009).
Zheng, S. C., Breeze, C. E., Beck, S. & Teschendorff, A. E. Identification of differentially methylated cell types in epigenome-wide association studies. Nat Methods 15, 1059-1066, doi:10.1038/s41592-018-0213-x (2018).
Graves, M. et al. Methylation differences at the HLA-DRB1 locus in CD4+ T-Cells are associated with multiple sclerosis. Mult Scler 20, 1033-1041, doi:10.1177/1352458513516529 (2013).
Ong, L. T. C., Schibeci, S. D., Fewings, N. L., Booth, D. R. & Parnell, G. P. Age-dependent VDR peak DNA methylation as a mechanism for latitude-dependent multiple sclerosis risk. Epigenetics Chromatin 14, 9, doi:10.1186/s13072-021-00383-x (2021).
Marabita, F. et al. Smoking induces DNA methylation changes in Multiple Sclerosis patients with exposure-response relationship. Sci Rep 7, 14589, doi:10.1038/s41598-017-14788-w (2017).
Christopoulos, P. F. et al. Targeting the Notch Signaling Pathway in Chronic Inflammatory Diseases. Front Immunol 12, 668207, doi:10.3389/fimmu.2021.668207 (2021).
Lee, W. S., Lee, W. H., Bae, Y. C. & Suk, K. Axon Guidance Molecules Guiding Neuroinflammation. Exp Neurobiol 28, 311-319, doi:10.5607/en.2019.28.3.311 (2019).
Ma, Q. et al. Specific hypomethylation programs underpin B cell activation in early multiple sclerosis. Proc Natl Acad Sci U S A 118, doi:10.1073/pnas.2111920118 (2021).
Sabatino, J. J., Jr., Zamvil, S. S. & Hauser, S. L. B-Cell Therapies in Multiple Sclerosis. Cold Spring Harb Perspect Med 9, doi:10.1101/cshperspect.a032037 (2019).
Bittner, S., Ruck, T., Wiendl, H., Grauer, O. M. & Meuth, S. G. Targeting B cells in relapsing-remitting multiple sclerosis: from pathophysiology to optimal clinical management. Ther Adv Neurol Disord 10, 51-66, doi:10.1177/1756285616666741 (2017).
Parisi, L. et al. Macrophage Polarization in Chronic Inflammatory Diseases: Killers or Builders? J Immunol Res 2018, 8917804, doi:10.1155/2018/8917804 (2018).
Bar-Or, A. et al. Analyses of all matrix metalloproteinase members in leukocytes emphasize monocytes as major inflammatory mediators in multiple sclerosis. Brain 126, 2738-2749, doi:10.1093/brain/awg285 (2003).
D'Amico, E. et al. Immunological Subsets Characterization in Newly Diagnosed Relapsing-Remitting Multiple Sclerosis. Front Immunol 13, 819136, doi:10.3389/fimmu.2022.819136 (2022).
Durelli, L. et al. T-Helper 17 Cells Expand in Multiple Sclerosis and Are Inhibited by Interferon-beta. Ann Neurol 65, 499-509, doi:10.1002/ana.21652 (2009).
Zhang, X. et al. IL-11 Induces Th17 Cell Responses in Patients with Early Relapsing-Remitting Multiple Sclerosis. J Immunol 194, 5139-5149, doi:10.4049/jimmunol.1401680 (2015).
Pender, M. P., Csurhes, P. A., Pfluger, C. M. M. & Burrows, S. R. Deficiency of CD8(+) effector memory T cells is an early and persistent feature of multiple sclerosis. Mult Scler J 20, 1825-1832, doi:10.1177/1352458514536252 (2014).
Caruana, P., Lemmert, K., Ribbons, K., Lea, R. & Lechner-Scott, J. Natural killer cell subpopulations are associated with MRI activity in a relapsing-remitting multiple sclerosis patient cohort from Australia. Mult Scler 23, 1479-1487, doi:10.1177/1352458516679267 (2017).
Rock, K. L., Reits, E. & Neefjes, J. Present Yourself! By MHC Class I and MHC Class II Molecules. Trends Immunol 37, 724-737, doi:10.1016/j.it.2016.08.010 (2016).
Salgado, F. J. et al. Interleukin-dependent modulation of HLA-DR expression on CD4and CD8 activated T cells. Immunol Cell Biol 80, 138-147, doi:10.1046/j.1440-1711.2002.01055.x (2002).
Diniz, S. N., da Silva, C. F., de Almeida, I. T., da Silva Costa, F. E. & de Oliveira, E. M. L. INFbeta treatment affects global DNA methylation in monocytes of patients with multiple sclerosis. J Neuroimmunol 355, 577563, doi:10.1016/j.jneuroim.2021.577563 (2021).
Souren, N. Y. et al. DNA methylation signatures of monozygotic twins clinically discordant for multiple sclerosis. Nat Commun 10, 2094, doi:10.1038/s41467-019-09984-3 (2019).
Ewing, E. et al. Combining evidence from four immune cell types identifies DNA methylation patterns that implicate functionally distinct pathways during multiple sclerosis progression. Mult Scler J 25, 439-439 (2019).
Taylor, B. V. et al. Latitudinal variation in incidence and type of first central nervous system demyelinating events. Mult Scler 16, 398-405, doi:10.1177/1352458509359724 (2010).
Yan, L. et al. OSAT: a tool for sample-to-batch allocations in genomics experiments. BMC Genomics 13, 689, doi:10.1186/1471-2164-13-689 (2012).
Morris, T. J. et al. ChAMP: 450k Chip Analysis Methylation Pipeline. Bioinformatics 30, 428-430, doi:10.1093/bioinformatics/btt684 (2014).
Tian, Y. et al. ChAMP: updated methylation analysis pipeline for Illumina BeadChips. Bioinformatics 33, 3982-3984, doi:10.1093/bioinformatics/btx513 (2017).
Teschendorff, A. E. et al. A beta-mixture quantile normalization method for correcting probe design bias in Illumina Infinium 450 k DNA methylation data. Bioinformatics 29, 189-196, doi:10.1093/bioinformatics/bts680 (2013).
Johnson, W. E., Li, C. & Rabinovic, A. Adjusting batch effects in microarray expression data using empirical Bayes methods. Biostatistics 8, 118-127, doi:10.1093/biostatistics/kxj037 (2007).
Peters, T. J. et al. De novo identification of differentially methylated regions in the human genome. Epigenet Chromatin 8, doi:Artn 6 10.1186/1756-8935-8-6 (2015).
Guo, Y. et al. Illumina human exome genotyping array clustering and quality control. Nat Protoc 9, 2643-2662, doi:10.1038/nprot.2014.174 (2014).
Purcell, S. et al. PLINK: a tool set for whole-genome association and population-based linkage analyses. Am J Hum Genet 81, 559-575, doi:10.1086/519795 (2007).
Zivkovic, M. et al. The tag SNP for HLA-DRB1*1501, rs3135388, is significantly associated with multiple sclerosis susceptibility: Cost-effective high-throughput detection by real-time PCR. Clin Chim Acta 406, 27-30, doi:10.1016/j.cca.2009.05.004 (2009).
Choi, S. W. & O'Reilly, P. F. PRSice-2: Polygenic Risk Score software for biobank-scale data. Gigascience 8, doi:10.1093/gigascience/giz082 (2019).
Sing, T., Sander, O., Beerenwinkel, N. & Lengauer, T. ROCR: visualizing classifier performance in R. Bioinformatics 21, 3940-3941, doi:10.1093/bioinformatics/bti623 (2005).
Teschendorff, A. E., Breeze, C. E., Zheng, S. J. C. & Beck, S. A comparison of reference-based algorithms for correcting cell-type heterogeneity in Epigenome-Wide Association Studies. Bmc Bioinformatics 18, doi:ARTN 105 10.1186/s12859-017-1511-5 (2017).
Dunning, M. J., Smith, M. L., Ritchie, M. E. & Tavare, S. beadarray: R classes and methods for Illumina bead-based data. Bioinformatics 23, 2183-2184, doi:10.1093/bioinformatics/btm311 (2007).
Aryee, M. J. et al. Minfi: a flexible and comprehensive Bioconductor package for the analysis of Infinium DNA methylation microarrays. Bioinformatics 30, 1363-1369, doi:10.1093/bioinformatics/btu049 (2014).
Morris, T. J. et al. ChAMP: 450k Chip Analysis Methylation Pipeline. Bioinformatics 30, 428-430, doi:10.1093/bioinformatics/btt684 (2014).
Marabita, F. et al. An evaluation of analysis pipelines for DNA methylation profiling using the Illumina HumanMethylation450 BeadChip platform. Epigenetics 8, 333-346, doi:10.4161/epi.24008 (2013).
Houseman, E. A., Molitor, J. & Marsit, C. J. Reference-free cell mixture adjustments in analysis of DNA methylation data. Bioinformatics 30, 1431-1439, doi:10.1093/bioinformatics/btu029 (2014).
Pidsley, R. et al. A data-driven approach to preprocessing Illumina 450K methylation array data. BMC Genomics 14, 293, doi:10.1186/1471-2164-14-293 (2013).
Leek, J. T., Johnson, W. E., Parker, H. S., Jaffe, A. E. & Storey, J. D. The sva package for removing batch effects and other unwanted variation in high-throughput experiments. Bioinformatics 28, 882-883, doi:10.1093/bioinformatics/bts034 (2012).
Ritchie, M. E. et al. limma powers differential expression analyses for RNA-sequencing and microarray studies. Nucleic Acids Res 43, e47, doi:10.1093/nar/gkv007 (2015).
Fabregat, A. et al. Reactome pathway analysis: a high-performance in-memory approach. Bmc Bioinformatics 18, 142, doi:10.1186/s12859-017-1559-2 (2017).
Orchard, S. et al. The MIntAct project—IntAct as a common curation platform for 11 molecular interaction databases. Nucleic Acids Research 42, D358-D363, doi:10.1093/nar/gkt1115 (2013).
Villanueva, R. A. M. & Chen, Z. J. ggplot2: Elegant Graphics for Data Analysis, 2nd edition. Meas-Interdiscip Res 17, 160-167, doi:10.1080/15366367.2019.1565254 (2019).

Competing interest reported. HB has received institutional (Monash University) funding from Biogen, F. Hoffmann-La Roche Ltd, Merck, Alexion, CSL, and Novartis; has carried out contracted research for Novartis, Merck, F. Hoffmann-La Roche Ltd and Biogen; has taken part in speakers’ bureaus for Biogen, Genzyme, UCB, Novartis, F. Hoffmann-La Roche Ltd and Merck; has received personal compensation from Oxford Health Policy Forum for the Brain. JLS received travel compensation from Biogen, Merck and Novartis; has been involved in clinical trials with Biogen, Merck, Novartis and Roche; her institution has received honoraria for talks and advisory board service from Biogen, Merck, Novartis and Roche. VEM has accepted honoraria for presentations and research funds from Biogen and Merck.

Download PDF

Version 1

posted

You are reading this latest preprint version

DNA methylation signatures of Multiple Sclerosis occur independently of genetic risk and are primarily attributed to B cells and monocytes

Status:

Version 1

Abstract

Figures

One Sentence Summary

Introduction

Results

Dna Methylation Analysis

Methylation differences act independently of genetic risk

Cell Specific Differential Methylation Is Mainly Attributed To B Cells And Monocytes

Over-representation analysis shows enrichment of DMPs at immunological pathways

Discussion

Conclusion

Materials And Methods

Patients’ enrolment and blood collection

Sample Groups

Whole Blood methylation:

DNA extraction and Quality Control

Methylation and genotyping arrays

Methylation analysis

Sensitivity Analysis

Genotype Analysis

Scores and receiver operator characteristic curves

Immune cell deconvolution analysis

Expression and DNA methylation Analysis

All cell types:

Correlation analysis (DNA methylation and gene expression)

metQTL analysis

Over-Representation Analysis

Abbreviations

Declarations

References

Additional Declarations

Supplementary Files

Status:

Version 1