Global DNA methylation changes in treated and untreated MS patients measured over time

Multiple sclerosis (MS) is an autoimmune, neurological disease. We investigated genome-wide DNA methylation profiles of CD4+ and CD8+ T cells from MS patients and healthy controls at baseline and a follow-up visit. Patients were all treatment-naïve at baseline, and either on treatment or remained untreated at the follow-up visit. MS patients show more changes in their T cell DNA methylation profiles as compared to healthy controls over time, with the most pronounced differences observed in the untreated MS patients. These findings underline the potential of DNA methylation as biomarkers in MS.


Introduction
Multiple sclerosis (MS) is a neurological condition characterised by varying degree of inflammation, demyelination and neurodegeneration of the central nervous system (CNS) (Dendrou et al., 2015). The underlying cause of MS is unknown; however, the onset of the disease is believed to result from genetic, epigenetic and environmental factors and interactions between these (Olsson et al., 2017). Genome-wide association studies (GWAS) during the last decades have identified 233 genetic risk variants associated with the disease (IMSGC, 2019a). The majority of these loci are mapped to genes of relevance for immune cell function, in both the adaptive and innate branch of the immune system (IMSGC, 2019b). T cells are considered important drivers for MS pathogenesis (Sawcer et al., 2011;Chitnis, 2007;Denic et al., 2013;Baecher-Allan et al., 2018), and many of the available disease modifying treatments (DMT) for MS target T cells (Linker et al., 2008). Both CD4 + and CD8 + T cells are found to be enriched in CNS lesions of MS patients (Hauser et al., 1986;Booss et al., 1983), and these cells are therefore considered to play a central role in the disease pathology.
Epigenetics is defined as heritable modifications of the DNA that do not affect the sequence of nucleotides. DNA methylation, a vastly studied epigenetic mark in somatic cells (Yin et al., 2017), is known to change with age, lifestyle, disease and medications (Csoka and Szyf, 2009;Martin and Fry, 2018;Bjornsson et al., 2008), suggesting that epigenetic factors could act as mechanisms through which environmental and genetic factors in MS interact. A major role of DNA methylation is the regulation of gene expression (Schübeler, 2015). Whilst a higher level of DNA methylation typically represses gene expression, there is also evidence that higher levels of DNA methylation can lead to increased gene expression (Neri et al., 2017). The majority of DNA methylation occurs at CpG sites of the DNA, and the functional effect of DNA methylation depends on the genomic location of the methylated residues and the cell type (Lorincz et al., 2004;Ball et al., 2009;Deaton et al., 2011;Aran et al., 2011).
Mapping the DNA methylation profiles at the genome-wide level of immune cells from MS patients longitudinally, could aid to increase the understanding of biological mechanisms involved in MS, and thereby support MS diagnosis, provide better guidance for therapeutic strategy and identify leads for the development of new therapeutic targets.
DNA methylation is extensively studied in MS (Huynh et al., 2014;Marabita et al., 2017;Maltby et al., 2017;Maltby et al., 2015;Kular et al., 2018;Rhead et al., 2018;Bos et al., 2015;Graves et al., 2014), and these studies support the hypothesis of the involvement of epigenetic factors in the disease pathology. These studies are cross-sectional (Huynh et al., 2014;Marabita et al., 2017;Maltby et al., 2017;Maltby et al., 2015;Kular et al., 2018;Rhead et al., 2018;Bos et al., 2015;Graves et al., 2014), there is currently little or no insight into how DNA methylation changes with disease duration and how MS treatment affects the DNA methylation profiles. We present a study of the genomewide DNA methylation changes in CD4 + and CD8 + T cells from relapsing-remitting MS (RRMS) patients over time, in both treated and untreated MS patients. To control for normal changes in DNA methylation, we included healthy controls sampled at two visits.

Study design
Thirty-five female RRMS patients were recruited from the Department of Neurology at Oslo University Hospital, and ten female healthy controls were recruited from research colleagues. Blood samples were collected at baseline (Rhead et al., 2018;Bos et al., 2015) and during the follow-up visit for each individual. MS patients were treatment-naïve at the baseline visit. The study is approved by the Norwegian Regional Ethical Committee for Medical and Health Research and participating MS patients and healthy controls have signed informed consent forms approved by the Norwegian Ethical Committee (REC).

Purification of CD4 + and CD8 + T cells
CD4 + and CD8 + T cells were purified using immunomagnetic cell separation kits (Human CD4 + negative selection kit and EasySep™ Human CD8 + positive selection kit and EasySep™, STEMCELL Technologies, Vancouver, Canada) as described in (Rhead et al., 2018;Bos et al., 2015). Purity of isolated T cells was confirmed to be at least 95% as measured by flow cytometry (Attune Acoustic Focusing Flow Cytometer, Life Technologies, Carlsbad, CA, USA). Antibodies used for purity assessment by flow cytometry were fluorescein isothiocyanateconjugated mouse anti-human CD4 (clone RTF-4 g, Southern Biotech, Birmingham, AL, USA), mouse anti-human CD8 (clone HIT8a, BD Biosciences, San Jose, CA, USA) and mouse IgG1 isotype control (15H6, Southern Biotech).

Assessment of DNA methylation, pre-processing and quality controls
The DNA methylation profiles for samples collected at baseline visits were obtained using the Infinium HumanMethylation450 (HM450) BeadChip (Illumina, San Diego, Ca, USA) and the Infinium Human-Methylation EPIC BeadChip (EPIC) as described in (Rhead et al., 2018;Bos et al., 2015). For the follow-up visit, DNA was extracted from the cells with QIAmp DNA mini kit (Qiagen, Hilden, Germany). Subsequently, 500 ng of genomic DNA was treated with sodium bisulphite prior to profiling the genome-wide DNA methylation on the EPIC BeadChip (Life & Brain GmbH, Bonn, Germany). The R package 'Minfi' (Aryee et al., 2014) was used for data processing. Raw methylation data from both visits was split by cell type for further processing. The function 'preprocessNoob' was applied for background subtraction and dye normalization, followed by 'preprocessQuantile' for quantile normalization. Probe positions with detection p-values >0.01 were set to missing, and probes with call rates of less than 95% were removed from analyses. Probes with common SNPs in the probe body or at the single nucleotide extension site, and probes predicted to cross-hybridize to other genomic regions (McCartney et al., 2016) were omitted from analysis.

Correction for confounders
The surrogate variable analysis (SVA) package in 'R' (Leek et al., 2012) was used to adjust for known and unknown sources of confounding, such as batch and chip effects. Following pre-processing and normalization of the data, SVA was applied on M-values (Du et al., 2010) of each sample group. In constructing the surrogate variables (SV), the visit (baseline or follow-up) was protected to avoid SVs capturing variation associated to this outcome variable. The constructed SVs were used as covariates in the subsequent downstream analyses and no other covariates were included in the model.

Statistical modelling and differential DNA methylation analysis
For each cell type, analysis of genome-wide DNA methylation changes between the two visits was performed for three separate groups: a) untreated patients, b) patients on treatment at follow-up and c) healthy controls. The statistical model considered the baseline and follow-up group as independent, and unpaired analyses were performed with the covariates estimated using the SVA method (Leek et al., 2012). Differentially methylated positions (DMPs) were identified using the empirical Bayes method in the R package 'limma' (Ritchie et al., 2015). The model used M-values for each CpG as the outcome with the visit (baseline or follow-up) as predictor and SVs as covariates. Bonferroni correction (Abdi, 2007) was applied to correct for multiple testing and adjusted p-values of <0.05 were considered significant.

Validation of differential methylation analysis by resampling
The healthy control sample group was smaller compared to the MS patient groups and therefore had less power to detect DMPs. In order to assess the extent of this difference in power on the number of the DMPs identified, a resampling method was used. For both cell types, each MS patient group was randomly resampled equal to the group size of the healthy controls (n = 10 for CD4 + T cells and n = 9 for CD8 + T cells). The process of calculating SVs and DMPs analyses for the two visits was repeated 100 times with the random matched selections of samples from MS patients collected at baseline and follow-up.

Genomic features
The 'UCSC_RefGene_Group' feature of the Illumina annotation file for HM450 (Hansen, 2016) was used for linking probes to genomic features. Probes missing a genomic feature annotation were assigned to the intergenic region (IGR).

Gene set testing on differentially methylated positions
Enrichment for KEGG (Kanehisa and Goto, 2000) pathways and Gene Ontotlogy (GO) categories (Ashburner et al., 2000;Nucleic Acids Res., 2021) were tested on the DMPs identified in our analyses. Gene set analysis was performed with the 'GOmeth' function (Maksimovic et al., 2021) available in the 'missMethyl' (Phipson et al., 2016) R package. The identifiers of significant DMPs were used as input in the pathway analysis whilst the complete list of all positions analysed was provided as background.

Results
We measured genome-wide DNA methylation of CD4 + and CD8 + T cells collected from RRMS patients (n = 35) and healthy controls (n = 10) at two visits. An overview of the study design and the analyses performed are illustrated in Fig. 1. Since DNA methylation changes over time, healthy controls were included as a normal reference (Bjornsson et al., 2008). The time differences in months between the baseline and follow-up visits for the patient groups (55.9 ± 16.2 and 32.5 ± 21.5 months for the untreated and treated patient groups, respectively) and healthy controls (65.7 ± 15.3 months) are shown in Fig. 2. An overall significant difference between the time interval between the three study groups (Kruskal-Wallis, p = 0.00073) was observed. Characteristics of the MS patients and healthy controls at the group level are provided in Table 1, and individual details are provided in table S1.
At baseline sampling, all MS patients were treatment-naïve, whereas at the follow-up visit, 20 patients (57%) had initiated DMT, and the remaining 15 patients (43%) were untreated. CD4 + T and CD8 + T cells were successfully isolated from respectively 17 and 19 treated patients. Two of the 15 patients in the untreated group had received DMT for brief periods between the two visits and discontinued treatment respectively CD4 + and CD8 + T cells were isolated from RRMS patients and healthy controls at two visits (baseline and follow-up). All collected DNA methylation data was split into either CD4 + T or CD8 + T cell data, prior to quantile normalization and differential methylation analysis in these analysis groups. I.S. Brorson et al. 5 years and 6 months prior to the follow-up visit. The remaining 13 patients stayed treatment-naïve between visits. Sample numbers for each group are listed in Table 2. No samples were omitted from analyses due to low call rate, and estimated gender based on methylation matched all participants recorded gender (female). Quality checks and normalization of DNA methylation data from CD4 + and CD8 + T cells were performed separately in the two data sets, prior to dividing into groups based on treatment status at the follow-up visit. Healthy controls' CD4 + and CD8 + T cells from baseline-and follow-up visit were treated as a separate group.

Differential methylated position analysis of baseline and follow-up samples
Analysis of differential DNA methylation changes were performed in CD4 + and CD8 + T cells from the three separate sample groups collected at two visits. DMPs with absolute mean Δbeta ranges between 0.045 and 0.087 were detected in all sample groups and both cell types between baseline and follow-up, listed in Table 3. The full lists of analysed CpGs for the three separate sample groups are listed in tables S2-S7. The lowest number of DMPs was found in CD4 + and CD8 + T cells from healthy controls, indicating that MS disease and/or DMT had a greater impact on DNA methylation than age. The highest number of DMPs was found in CD4 + and CD8 + T cells from the group of untreated MS patients at the follow-up visit, followed by the number of DMPs in CD4 + and CD8 + T cells from the treated MS patient group. The identified DMPs were further examined for overlapping sites across the MS patient groups, summarised in the Venn diagrams in Fig. 3. The highest proportion of unique DMPs was found in CD4 + T cells from untreated patients, which accounts for 57% of the detected DMPs in that group. The lowest proportion of unique DMPs was in CD4 + T cells from the treated patient group, representing 22% of the total number of DMPs in that group.

Validation of differential methylation analysis by resampling
To rule out that the differences in observed DMPs between the MS patients and healthy controls were attributable to lower sample size, and thus lower power to detect differences, we performed validation analyses. For both CD4 + and CD8 + T cells, 100 random selections of respectively ten and nine samples from MS patients were analysed. This analysis confirmed that a higher number of DMPs is detected in samples from untreated MS patients at the follow-up visit compared to treated MS patients. The number of significant DMPs detected in the full   HC: healthy controls, CD4 + T cells: CD4 + T cell DNA methylation data, CD8 + T cells: CD8 + T cell DNA methylation data. Numbers indicate samples collected at both visits. a Patients were treatment-naïve at baseline and untreated at follow-up. b Patients were treatment-naïve at baseline while on DMT at follow-up. c Healthy controls included at baseline and follow-up visit. analysis, and the mean number of significant DMPs obtained from the resampling analysis, are listed in Table 3. Of note, the number of DMPs in the healthy control group was higher than the mean number of DMPs detected in the treated MS group in these resampling analyses.

Proportion of hyper-and hypomethylated DMPs between baseline and follow-up
Since the level of DNA methylation affects gene expression, we examined the proportions of hyper-and hypomethylated DMPs (Fig. 4). We observed a higher proportion of hypermethylated DMPs in CD8 + T cells from MS patients at the follow-up visit (Fig. 4B), with respectively 36% and 32% hypermethylated DMPs in CD8 + T cells shared between the treated and untreated patient groups. The same was observed for CD4 + T cells from untreated MS patients, whereas DMPs in CD4 + T from treated patients showed no global differences in direction of effect (Fig. 4A). Next, we investigated the genomic features annotated to the identified DMPs, shown in Fig. 5 (Supplementary Figure 1 for the healthy controls). These genomic features were grouped based on their location relative to genes, i.e. within 1500-200 base pairs upstream of the transcriptional start site (TSS; TSS1500), 200 base pairs upstream from the TSS (TSS200), first exon of a gene (1stExon), 5 ′ and 3 ′ untranslated region (UTR), within gene body (body) and intergenic region (IGR). DMPs detected in CD4 + T cells from both the treated and untreated MS patient groups were more represented in the regulatory region of a gene, i.e. close to the TSS, in the 5'UTR as well as in the first exon of a gene. For both patient groups, these regions showed a higher proportion of hypermethylation. The CD8 + T cells from the treated MS patients (Fig. 5c), showed predominantly hypermethylation in the regulatory regions, the first exons as well as the IGR. The untreated patient group (Fig. 5D) showed the same trend, though not to the same extent as the treated patient group. For all subgroups, the lowest proportions of DMPs were observed within the gene body and the 3'UTR.

Gene set enrichment analysis of differentially methylated positions
The genes annotated to the identified DMPs were tested for enrichment of specific molecular pathways in each sample group. The 'GOmeth' (Maksimovic et al., 2021) function of the 'MissMethyl' (Phipson et al., 2016) package was used, and enrichment analysis was performed for KEGG pathways and GO categories. DMPs from the treated patient groups' CD8 + T cells were annotated to genes enriched in the T cell receptor signalling pathway (FDR = 0.03) and the Apoptosis pathway (FDR = 0.03), in addition to three GO categories, two representing cellular component processes and one biological process. No enrichment for KEGG pathways or GO categories in the treated groups' CD4 + T cell or healthy controls were found. DMPs identified in the  untreated patient groups' CD8 + T cells were annotated to genes enriched in the RNA transport pathway (FDR ≤ 0.01), as well as 31 GO categories, with the majority representing biological processes (n = 16), followed by cellular components (n = 12) and molecular functions. Untreated patients' CD4 + T cells DMPs' were annotated to genes enriched in the Cell cycle (FDR = 0.02) KEGG pathway, as well as four biological processes and four cellular component GO categories. The complete lists of pathways and GO categories analysed are provided in tables S8-S19.

Discussion
In this study, we performed analysis of DNA methylation data from MS patients and healthy controls sampled at two time points, with the covariates estimated using the SVA method (Leek et al., 2012). Unpaired analysis was performed as the surrogate variables effectively captured the pairwise dependency of the samples. We found that the genomewide DNA methylation profiles in T cells from MS patients exhibited significant changes over time. The treated MS patient group showed fewer significant DMPs between the visits compared to the untreated MS patients, further validated by resampling analysis. The MS patients were sampled at two visits, and all patients included were treatment-naïve at the baseline, which means that no immunomodulating treatment affected the epigenetic landscape at that time.
Whilst the patients were treatment-naïve at baseline, they have been diagnosed on average 8.3 years beforehand, which might have influenced the level of DNA methylation measurements at that time. Though studies show that chip-based DNA methylation measurements are biased in both sample placement on the chip as well as chip type (Nygaard et al., 2016), these confounding effects are corrected for in the present study by random placement of both patient and healthy control samples in each sampling batch, in addition to include surrogate variables as covariates in the regression models. Batch effects introduced by the separate chips, i.e. follow-up samples were run on a separate chip as compared to baseline samples, could affect measurements in the respective study groups. However, this would not affect directly the comparisons made, since both patients and healthy controls were present in every sampling batch. Future longitudinal studies investigating the DNA methylation effects of treatment in MS patients will benefit from samples before and after initiation of treatment, from study designs limiting batch effects with a uniform time interval between samplings. A strength of this study is the inclusion of healthy controls, which served as normal reference for age-related changes in DNA methylation profiles. Healthy controls displayed a lower number of significant DMPs than the MS patient groups at the follow-up visit compared to baseline. We investigated whether this observation was attributable to the larger sample size of MS patients, and thus more power to detect differences. In these resampling analyses, we confirmed a high number of significant DMPs in both cell types from the untreated MS patient group at the follow-up visit compared to baseline, indicating that this is of true biological relevance and not due to larger sample size in the patient groups. We observed a higher number of DMPs among the healthy controls as compared to the treated MS patient group in the resampling analysis. As DNA methylation levels have been shown to change over time (Bjornsson et al., 2008), the latter could be an outcome of the significantly longer time range between samplings in the healthy control group compared to the treated patient group. However, we cannot exclude that the observation is attributable to the presence of MS. The difference in the number of DMPs in the treated and untreated patient groups in both cell types at the follow-up visit cannot be explained by a higher power for the larger MS groups, as the random resampling of the MS patients matching the number of healthy controls showed the same effect.
These results indicates that the epigenetic configuration of T cells is associated with treatment status in MS patients. Albeit the differences in DNA methylation levels might indicate a biological effect due to the presence of MS, beyond the influence of time on DNA methylation, the possibility of confounding factors underlying the observed differences cannot be excluded. The difference in number of detected DMPs in the treated and untreated group can be attributable to the severe treatment heterogeneity and the time interval between the visits in the treated MS patient group compared to the untreated patient group.
Interestingly, there was an abundance of hypermethylation among significant DMPs at the follow-up visit compared to baseline in CD8 + T cells from both the treated and untreated patient groups. We have in previous studies identified genome-wide hypermethylation in CD8 + T cells from MS patients as compared to healthy controls (Rhead et al., 2018;Bos et al., 2015). In agreement with these studies, the current results indicate distinctive hypermethylation of MS patients' CD8 + T cells after time with the disease, regardless of treatment status. Previous studies have found hypermethylated DMPs in patients treated with dimethyl fumarate, with effects observed in CD4 + T cells (Ntranos et al., 2019;Maltby et al., 2018). In the present study, only one patient was on fumarate treatment, with the remaining patients receiving several types of DMT, and we did not observe more hypermethylated DMPs in CD4 + T cells in the treated group. Future studies on the methylome and treatment status in T cells from MS patients should be investigated in study groups that have received the same type of DMT, preferably for the same period of time.
Mapping of the identified DMPs to genomic regions in all patient groups and both cell types showed a higher proportion of DMPs in the gene regulatory regions as well as the first exon compared to the gene body and 3'UTR. In particular, the first exon had the most pronounced hypermethylation in all groups. The first exon is known for a distinct structure and is richer in CpG sites than other internal exons and introns (Kalari et al., 2006;Majewski and Ott, 2002). The least differential methylation across all sample groups was observed within the gene body and the 3'UTR. DNA methylation of promotor regions is associated with repression of gene expression, whereas intragenic DNA methylation is associated with induced gene expression (Neri et al., 2017). The features of differences in DNA methylation presented here suggests a functional role for the observed DNA methylation and should be further studied in larger longitudinal studies of MS. Furthermore, pathway analysis of the significant DMPs showed that treated MS patients' CD8 + T cells had an enrichment of DMPs in the T cell receptor signalling pathway, which is essential for an efficient immune response, and the apoptosis pathway, that ensures the elimination of damaged or redundant cells. Untreated MS patients' CD8 + T cells had an enrichment of DMPs in the RNA transport pathway, which is crucial for gene expression. The untreated patient groups' CD4 + T cells had enrichments of DMPs in the cell cycle pathway, in the classification of cellular processes of cell growth and death. DMPs from the treated patients' CD4 + T cells, as well as the healthy controls, had no enrichments in KEGG pathways or GO categories. There were no immune specific pathways in the GO categories identified by the enrichment of DMPs in treated and untreated MS patients' CD8 + T cells and the untreated MS patient CD4 + T cells. There is a potential for bias when performing gene set enrichment analysis on DNA methylation data due to the different numbers of CpGs per gene (Geeleher et al., 2013). 'GOmeth' takes into account that CpGs are not evenly distributed throughout the genome and corrects for multiple CpGs that occur within a single gene, as well as that single CpGs can be associated with multiple genes (Maksimovic et al., 2021). The test performed by 'GOmeth' weighs CpGs according to the number of genes they map to, which is meant to decrease the erroneous overrepresentation of associated pathways of significance.
We cannot exclude that the observed differences in DNA methylation are attributable to changes in the composition of subtypes for the CD4 + and CD8 + T cells. This can be addressed by longitudinal sampling of MS patients with a more detailed cellular phenotyping of the CD4 + and CD8 + cellular fractions. To better understand the dynamics in the observed changes in DNA methylation associated to the presence of MS, longitudinal studies with bigger study groups and multiple time intervals should be performed.

Conclusion
Our study suggests that the genome-wide DNA methylation landscape of MS patients is affected both by medications and the presence of MS itself. Furthermore, CD8 + T cells from MS patients are predominantly hypermethylated at the time of follow-up compared to baseline. These findings indicate an important role for DNA methylation of T cells in MS pathological processes and provides further evidence for a link between the environment and genetic makeup through DNA methylation.

Availability of data and materials
The raw datasets generated during the presented study cannot be deposited in publicly accesible databases due to restriction imposed by the Norwegian law on privacy. Requests for data access for the purpose of MS research can be addressed to the Oslo University Hospital grants@ ous-hf.no. A process of data transfer agreement signing will be initiated upon approval of the request by the study leaders and legal department.