Hypermethylation of MIR21 in CD4+ T cells from patients with relapsing-remitting multiple sclerosis associates with lower miRNA-21 levels and concomitant up-regulation of its target genes

Background: Multiple sclerosis (MS) is a chronic inflammatory disease of the central nervous system caused by genetic and environmental factors. DNA methylation, an epigenetic mechanism that controls genome activity, may provide a link between genetic and environmental risk factors. Objective: We sought to identify DNA methylation changes in CD4+ T cells in patients with relapsing-remitting (RR-MS) and secondary-progressive (SP-MS) disease and healthy controls (HC). Methods: We performed DNA methylation analysis in CD4+ T cells from RR-MS, SP-MS, and HC and associated identified changes with the nearby risk allele, smoking, age, and gene expression. Results: We observed significant methylation differences in the VMP1/MIR21 locus, with RR-MS displaying higher methylation compared to SP-MS and HC. VMP1/MIR21 methylation did not correlate with a known MS risk variant in VMP1 or smoking but displayed a significant negative correlation with age and the levels of mature miR-21 in CD4+ T cells. Accordingly, RR-MS displayed lower levels of miR-21 compared to SP-MS, which might reflect differences in age between the groups, and healthy individuals and a significant enrichment of up-regulated miR-21 target genes. Conclusion: Disease-related changes in epigenetic marking of MIR21 in RR-MS lead to differences in miR-21 expression with a consequence on miR-21 target genes.


Preparation of CD4+ T cells
For the discovery cohort peripheral blood mononuclear cells (PBMCs) were isolated and collected in sodium heparin tubes using a standard Ficoll (GE Healthcare) procedure. Cells were separated using density gradient centrifugation, collected from the interphase, washed twice in Dulbecco's phosphate buffered saline and prepared for cell sorting. Sorting of the CD4 + T cell population was performed by adding fluorochrome-conjugated antibodies against human CD4 (Becton Dickinson) and CD3 (BD Bioscience) using a MoFlo high-speed cell sorter (Beckman Coulter).
For the validation cohort PBMCs were isolated from peripheral blood using sodium citratecontaining cell preparation tubes (BD Vacutainer™ CPT™ Tube, Becton Dickinson), Sorting of the CD4+ T cell population was performed by adding microbeads against human CD4 using an autoMACS® cell separator (Miltenyi Biotec).
Directly after sorting, cell pellets were frozen and kept at -80C° until DNA/RNA extraction.

DNA extraction
Extraction of genomic DNA was carried out using a Gen Elute Mammalian Genomic DNA Miniprep kit (Sigma-Aldrich). The amount and quality of DNA was accessed by NanoDrop ND-1000 Spectrophotometer (NanoDrop Technologies).

DNA methylation analysis
DNA methylation was profiled using the Infinium HumanMethylation450 BeadChip (Illumina) arrays at the bioinformatics and expression analysis core facility (BEA) at Karolinska Institutet.
Methylation data was individually analyzed in R using the Minfi and ChAMP package, probes containing SNPs within 2 bp of the CpG were removed, and type 1 and type 2 probes were normalized using quantile normalization and BMIQ 5 . The sex of the patients was confirmed using the GetSex function from the Minfi package and the cell type was confirmed using the cell type deconvolution method from Minfi based on the Houseman algorithm 6 . Significant cofounders were identified using PCA. Batch effects were corrected using ComBat from the SVA package 7 . Differentially methylated positions (DMPs) were determined by using the limma package 8 applying a linear modeling that included MS status, age and sex as co-variates.
Differences were calculated between HC and RR-MS and RR-MS and SP-MS. Additionally any significant differences between all three groups were identified using the moderated F-statistic test included in eBayes 9 .

Meta-analysis
Meta-analysis on the 11 CpG probes was performed by combining the outcomes of the analyses between RR-MS and HC using -values from the cohorts from Sweden (n=24), Norway (n=30) and Australia (n=40). Comparison of RR-MS with HC was performed on -values and corrected for age in the Norwegian and Australian cohort (comprising females only), and age and sex in the Swedish cohort.
Due to the strong impact of age on the methylation in the MIR21 region only individuals with available age information and age-matched HC where retained for analysis. Similar data were obtained with all individuals (data not shown). We used two different meta-analysis methodologies: (1) the "summation of p value" method for combining p-values 10 , and (2) the effect-size based meta-analysis using a fixed effect model (estimated by restricted maximumlikelihood) since there was no evidence of significant heterogeneity 11 . Similar outcome to "summation of p value" was obtained using the Fisher's test for combining p-values (data not shown). To conduct the meta-analysis we used the R packages metaphor and metap.

RNA extraction
For subsequent qPCR analysis total RNA from the discovery and independent sample cohort was isolated using standard TRIzol protocol (Invitrogen) and Allprep Total RNA/DNA Kit (Qiagen), respectively, according to manufacturer´s recommendations. RNA concentration and purity were determined by measurement of A260/A280 ratios with a NanoDrop ND-1000 Spectrophotometer (NanoDrop Technologies). For subsequent transcriptome analysis RNA, was extracted using the miRNeasy Kit (Qiagen) according to manufacturer´s recommendations. The RNA integrity was determined using Bioanalyzer (Agilent Technologies). All samples had an RNA integrity number (RIN) above 8. All RNA samples obtained from the different preparations were immediately frozen and stored at −80°C until further use.

Quantitative real-time PCR analysis
Expression levels of miR-21 were determined by quantitative real-time PCR (qRT-PCR) using TaqMan  for every sample and quantified using the 2−ΔΔCT method. Statistical analysis between the groups was performed using Student's t-test (for two groups) and ANOVA with Bonferroni correction for selected groups (for more than two groups), and correlation between expression and methylation levels was assessed using Pearson's correlation test in GraphPad Prism 5 (GraphPad).

Genotyping
Allelic discrimination of rs8070345 (in the VMP1 locus) was performed using a predesigned

Association and correlation analysis
In the GOLDN cohort we used the lmekin function of the R kinship package to fit linear mixed effects models with the methylation β-value at each of the interrogated CpG sites as the outcome and the following predictors: SNP, age, sex, study site (Utah vs. Minnesota), current smoking (yes vs. no), body mass index, 4 principal components capturing T-cell purity 3 (as fixed effects), and family (as a random effect). We additionally fit linear regression models with the methylation β-value at each CpG site as the outcome and the non-genetic predictors, namely age, smoking, and BMI, and the T-cell purity principal components.
In the Multiple Sclerosis cohort we fit linear regression model with the methylation β-value at each of the interrogated CpG sites as the outcome and the following predictors: SNP, age and sex. We additionally fit linear regression model with the methylation β-value at each of the interrogated CpG sites as the outcome and the following predictors: age, sex, Multiple Sclerosis Severity Score (MSSS) and lymphocyte count. All analyses were done in Rcmd.

Transcriptome analysis
Total RNA (500ng) with RIN above 8. Sequencing was carried out on the Illumina HiSeq 2500 to generate 75bp Paired End Data with an average of 10M reads above Q30. The sequence reads were mapped to hg19 reference with Tophat2 12 and HTSeq 13 was used to quantify counts-per-gene. The expression data from CD4 + cells (34 samples) were normalized with Conditional Quantile Normalization (CQN) method 14 followed by batch reduction by ComBat 7 . Differential expression analysis was conducted using limma tool 8 in order to identify differences between disease types.

Ingenuity Pathways Analysis
Two list of genes, 1) differentially expressed genes between RR-MS and HC (p<0.05) and 2) TarBase7.0-predicted miR-21 target genes up-regulated in RR-MS compared to HC (p<0.05) from RNAseq analysis were uploaded to the Ingenuity Pathways Analysis platform (Qiagen). The first list was used to assess significant up-stream regulators that can explain observed gene expression changes estimated by the overlap p-value, which is calculated using Fisher's exact test and significance is generally attributed to p<0.01. The activation z-score is used to infer the activation states of predicted upstream regulators and z<-2 and z>2 indicate significantly inhibited and activated upstream regulators, respectively. The second list of genes was used to infer biological functions that might be affected by miR-21. Right-tailed Fisher's exact test was used to calculate a p-value determining the probability that each biological function assigned to that data set is due to chance alone. Benjamini-Hochberg correction for multiple testing was used to calculate significant p-values.