Abstract Background Myalgic Encephalomyelitis/Chronic Fatigue Syndrome (ME/CFS) is a complex condition involving multiple organ systems and characterized by persistent/relapsing debilitating fatigue, immune dysfunction, neurological problems, and other symptoms not curable for at least 6 months. Disruption of DNA methylation patterns has been tied to various immune and neurological diseases; however, its status in ME/CFS remains uncertain. Our study aimed at identifying changes in the DNA methylation patterns that associate with ME/CFS. Methods We extracted genomic DNA from peripheral blood mononuclear cells from 13 ME/CFS study subjects and 12 healthy controls and measured global DNA methylation by ELISA-like method and site-specific methylation status using Illumina MethylationEPIC microarrays. Pyrosequencing validation included 33 ME/CFS cases and 31 controls from two geographically distant cohorts. Results Global DNA methylation levels of ME/CFS cases were similar to those of controls. However, microarray-based approach allowed detection of 17,296 differentially methylated CpG sites in 6,368 genes across regulatory elements and within coding regions of genes. Analysis of DNA methylation in promoter regions revealed 307 differentially methylated promoters. Ingenuity pathway analysis indicated that genes associated with differentially methylated promoters participated in at least 15 different pathways mostly related to cell signaling with a strong immune component. Conclusions This is the first study that has explored genome-wide epigenetic changes associated with ME/CFS using the advanced Illumina MethylationEPIC microarrays covering about 850,000 CpG sites in two geographically distant cohorts of ME/CFS cases and matched controls. Our results are aligned with previous studies that indicate a dysregulation of the immune system in ME/CFS. They also suggest a potential role of epigenetic de-regulation in the pathobiology of ME/CFS. We propose screening of larger cohorts of ME/CFS cases to determine the external validity of these epigenetic changes in order to implement them as possible diagnostic markers in clinical setting.

Citation: Trivedi MS, Oltra E, Sarria L, Rose N, Beljanski V, Fletcher MA, et al. (2018) Identification of Myalgic Encephalomyelitis/Chronic Fatigue Syndrome-associated DNA methylation patterns. PLoS ONE 13(7): e0201066. https://doi.org/10.1371/journal.pone.0201066 Editor: Dajun Deng, Beijing Cancer Hospital, CHINA Received: March 23, 2018; Accepted: July 7, 2018; Published: July 23, 2018 Copyright: © 2018 Trivedi et al. This is an open access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited. Data Availability: All relevant data are within the paper and its Supporting Information files. All raw data are available from the Gene Expression Omnibus (GEO) database under the accession number GSE111183. Funding: This study was supported by the grant R15NS087604-01A1 from the National Institute of Neurological Disorders and Stroke of the National Institutes of Health and by Presidential Faculty Research and Development Grant from Nova Southeastern University (both awarded to LN). The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript. Competing interests: The authors have declared that no competing interests exist.

1. Introduction Myalgic Encephalomyelitis/Chronic Fatigue Syndrome (ME/CFS) is a condition that is characterized by an abrupt or delayed onset of persistent/relapsing symptomatology including memory and other neurological problems, muscle and joint pain, gastrointestinal issues, hormonal imbalance, immune dysfunction and debilitating fatigue. Moreover, such symptoms are usually unresolved with bed rest and are severe enough to impair average daily activity below 50 percent of usual activity level, lasting for a period of at least six months [1]. While the mechanism of ME/CFS remains unclear and diagnostic methods exclusively rely on symptomatology presentation and exclusion of laboratory findings, research efforts have demonstrated that ME/CFS impacts the endocrine, neurological, immune and metabolic processes resulting in impaired physiological homeostasis [2–4]. Statistical studies estimate the prevalence of ME/CFS at 0.23 to 0.41 percent [5, 6] of the general population, with a female to male ratio of 6:1 [7]. With this prevalence, annual costs to the United States economy have been estimated at $9 billion in lost productivity and up to $24 billion in health care expenditures [8–10]. Therefore, it seems that ME/CFS not only impacts an individual’s overall well-being and quality of life, but it also has far reaching effects on the society and economy and constitutes a significant public health concern. Currently, treatment of ME/CFS relies only on the management of symptomatology [11] and improvement in quality of life due to a lack of understanding of the mechanisms underpinning disease onset and progression, limiting treatment options to partial and/or temporary relief of symptoms [11]. While some advances have been made in identifying molecular changes associated with ME/CFS, its complexity and the involvement of multiple organ systems have hindered the exact causes of the disease [12]. An improved understanding of the key molecular mechanisms of ME/CFS and dysfunction within regulatory systems will translate into appropriate diagnostic methods and management of cases, providing more targeted approaches to treatment. Disruption of epigenetic mechanisms is linked to various immune, neurological and endocrine diseases [13–15]. Furthermore, DNA methylation patterns were found to be altered in several diseases often reported as comorbid to ME/CFS such as fibromyalgia (FM) and irritable bowel syndrome (IBS) [16, 17]. With respect to ME/CFS, we are aware of only a few studies, which examined differences in DNA methylation patterns between ME/CFS cases and controls [18–20]. These studies used Illumina Human Methylation450 BeadChip microarrays, which allow to analyze over 450,000 methylation sites per sample at single-nucleotide resolution. Other two additional studies limited the analysis to specific gene promoter regions using a site-specific approach for measuring DNA methylation in ME/CFS cases [21, 22]. The recently released Illumina MethylationEPIC microarrays allow analysis of DNA methylation changes on over 850,000 CpG sites [23]. This additional coverage should facilitate uncovering additional changes in transcription regulation present in ME/CFS subjects. In addition, by validating DNA methylation patterns associating with ME/CFS in 64 participants from two geographically distant independent cohorts using pyrosequencing, we identified consensus CpG hypomethylated sites which could be used in future screenings of associations of these epigenetic changes in ME/CFS. Therefore, this study intended to validate and build upon the previous analysis of genome-wide DNA methylation changes in ME/CFS and extend such findings by utilizing larger coverage of CpG sites. Future validation studies in larger cohorts of ME/CFS cases are warranted to provide reliable epigenetic markers in ME/CFS.

2. Materials and methods 2.1. Sample collection and processing 2.1.1. Cohort recruitment. ME/CFS cases and controls were recruited from two locations: the Miami/Fort Lauderdale (Florida) area and Valencia, Spain as part of larger biomarker-oriented studies. All subjects had comparable age and body mass index (BMI) (Table 1). Because females are affected by ME/CFS more than males (6:1 ratio) [7], we recruited only female participants to reduce potential confounding effects. For inclusion we utilized the Fukuda [24] and Canadian [25] case definitions. Fukuda requiring fatigue of greater than 6 months duration and 4 of 8 symptom criteria including exercise induced relapse, myalgia, arthralgia, headache of a new and different type, non-restorative sleep, cognitive complaints, sore throat and tender lymph nodes [24], and Canadian requiring exercise induced relapse and symptoms form at least of 3 categories (immune, autonomic, endocrine) [25]. PPT PowerPoint slide

PowerPoint slide PNG larger image

larger image TIFF original image Download: Table 1. Demographic information and SF-36 results for all ME/CFS cases (n = 33) and controls (n = 31) that participated in DNA methylation analysis and validation. *—p<0.05, Student’s t-test, ME/CFS cases versus controls subjects. Data are shown as mean ± standard error of mean. https://doi.org/10.1371/journal.pone.0201066.t001 All ME/CFS study subjects had the Medical Outcomes Study 36-item short-form survey (SF-36) summary physical score below the 50th percentile, based on population norms. All subjects were between 30 and 60 years old when samples were collected. Exclusion criteria: all ME/CFS subjects, selected by their medical history and physical examination, had no history of malignancy or other systemic disorders including the listed psychiatric exclusions required for a diagnosis of ME/CFS, as clarified by Reeves [26], using the Composite International Diagnostic Instrument (a computerized interview format assessment of comorbid or exclusionary psychiatric diagnoses) [27]. Additional exclusion criteria included: active smoking or alcohol history, taking medications that could potentially impact immune function (e.g. steroids, immune-suppressors, etc.), taking diuretics, beta blockers, calcium channel blockers, ACE inhibitors, as well as subjects with a BMI > 35. Controls were self-defined as healthy, sedentary (no regular exercise program, sedentary employment), and matched to ME/CFS cases by age (+/- 5 years), race/ethnicity and BMI (+/- 5). All subjects were asked to complete the Gynecologic Questionnaire to assess routine gynecologic parameters and were asked to come for the assessment and collection of blood during the first two weeks of their menstrual cycle if premenopausal. All subjects signed an informed consent approved by the Institutional Review Board of the Nova Southeastern University. Ethics review and approval for the subjects from Valencia was obtained from Hospital de La Plana (Villarreal, Spain) Ethics Committee. All subjects from Valencia (Spain) signed an informed consent before they could be included in the study. The following questionnaires were completed for all enrolled subjects (cases and controls): the Multidimensional Fatigue Inventory (MFI) [28] to measure fatigue, the SF-36 [29] to assess health-related quality of life, the Krupp Fatigue Severity Inventory (Krupp FSI) [30] to measure perception of fatigue severity, the Sickness Impact Profile [31] to measure the impact of symptoms on the activities of daily life. Student t-test was used to assess statistical difference between ME/CFS cases and controls in SF-36 scores for physical and mental health. 2.1.3. Isolation of PBMCs and DNA extractions. Blood samples (8 ml) collected and kept at room temperature in K 2 EDTA (Becton Dickinson) were processed within 2 h by dilution at 1:1 (v/v) ratio in phosphate-buffered saline solution (PBS), layering on top of 1 volume of Ficoll-Paque Premium (GE Healthcare) and separation by density centrifugation at 20°C at 500 g for 30 min (brakes off). The PBMC layer was isolated and washed with PBS. The isolated PBMC pellets were resuspended in 1 volume of red blood cell lysis buffer (155 mM NH 4 Cl, 10 mM NaHCO 3 , pH 7.4, and 0.1 mM EDTA), kept on ice for 5 min, and centrifuged (20°C, 500 g, 10 min). The washed pellets were resuspended in freezing medium (90% FBS, 10% DMSO) adjusting their concentration to 107 cells/ml, aliquoted and frozen at -150°C or liquid nitrogen until use. Genomic DNA was extracted from rapidly thawed PBMCs using DNeasy Blood & Tissue Kit (Qiagen), according to manufacturer´s instructions. In brief, pelleted PBMCs were resuspended in 200 μl of PBS supplemented with proteinase K and RNase A and genomic DNA was obtained from treated lysates using the kit binding columns and binding/wash/elution standard conditions. DNA quality and concentration were assessed by Agilent TapeStation 4200 (Agilent Technologies). All DNA samples had DNA Integrity Number (DIN) above 8 (data not shown). 2.2. Genomic DNA methylation profiling For global epigenetic measurements MethylFlash Methylated DNA 5-mC Quantification Kit (EpiGentek) and MethylFlash Methylated DNA 5-hmC Quantification Kit (EpiGentek) have been used with 100 ng of genomic DNA from each sample. DNA methylation and DNA hydroxy-methylation were quantified using 5-methylcytosine and 5-hydroxymetylcytosine respectively by an enzyme-linked immunosorbent assay-like reaction as previously described [32, 33]. The levels of methylated DNA were calculated using the absorbance on a microplate reader at 450 nm. Results were normalized against a standard curve prepared according the manufacturer’s protocol using reference 0–100% methylated standards. For genome-wide DNA methylation assessment 500 ng of genomic DNA from each sample was submitted to the Center of Genome Technology of the John P. Hussman Institute for Human Genomics in the Miller School of Medicine University of Miami. Genomic DNA was bisulfite converted using the EZ DNA Methylation Kit (Zymo Research), and after processing, according to Illumina’s specifications, DNA was hybridized to Illumina MethylationEPIC microarrays [34]. 2.3. Analysis of genomic DNA methylation data Signal intensity data (IDAT files) were imported into R/Bioconductor (v.3.3.1) package RnBeads (v.1.2.0) [35]. Methylation values for each of the probes on the Illumina MethylationEPIC microarray were produced as beta-values, calculated as the ratio of methylated probe intensity over total intensity (methylated and unmethylated) for each probe, which range from 0 to 1, roughly corresponding to the methylation percentage of the probe [35]. Low-quality data (probes with detection p-value > 0.01) were removed from the analysis. Probes outside of the context (non-CpG) and invariable through all sample probes were filtered out. Probes overlapping SNPs were removed using RnBeads function “rnb.execute.snp.removal” with the option “snp =“any””. The resulting data set (649,836 probes) was normalized using BMIQ normalization of RnBeads package [36] and used for exploratory and differential methylation analysis. Differential methylation analysis was performed by using two strategies: genomic region-based and site-based. While the genomic region-based approach looks for the average methylation level of the genomic region (CpG islands, promoters, genes, and genome-wide tiling regions) and then compares the methylation status of the regions across samples, site-based approach analyses the CpG loci are analyzed one at a time. To detect regions or loci that are differentially methylated in ME/CFS cases compared to controls we used R/Bioconductor limma package [37] incorporated into RnBeads analysis. CpG sites were considered to be differentially methylated if they met the following selection criteria: the absolute beta-value difference between the mean beta-values of cases and controls was greater than 0.05, and false discovery rate (FDR) ≤ 0.05 using the Benjamini-Hochberg procedure [38] for microarray-based analysis and p≤ 0.05 according to the Student t-test for pyrosequencing-based analysis. Promoter regions (1500 bp upstream of the transcription start site (TSS) and 500 bp downstream of TSS) were considered to be differentially methylated if they met the following criteria: FDR ≤ 0.1 according to Benjamini-Hochberg procedure [38] and the absolute beta-value difference between the mean beta-values of cases and controls was greater than 0.05. RnBeads calculates the combined rank score for differential DNA methylation to each genomic region. This combined rank is defined as the maximum (worst) of three individual rankings of 1) by absolute difference in mean DNA methylation levels in the genomic region, 2) by the relative difference in mean DNA methylation levels, which is calculated as the absolute value of the logarithm of the quotient of mean DNA methylation levels, and 3) by the CpG-based or genomic region-based p-value that is combined from CpG-specific p-values of the region using an extension of Fisher’s method. We present our data in the Supplementary Tables ordered by the combined rank. The raw microarray data have been deposited in Gene Expression Omnibus (GEO) at NCBI (https://www.ncbi.nlm.nih.gov/geo/) under the accession number GSE111183. 2.4. Validation of differentially methylated promoters by pyrosequencing For validation of the data obtained with the Illumina MethylationEPIC microarrays by pyrosequencing, genomic DNA was treated with EZ DNA Methylation-Gold Kit (Zymo Research) according to the manufacturer’s protocol. Methylation of DNA was quantified with bisulfite treatment of DNA and simultaneous PCR. A 50-μl PCR was performed in 25 μl of GoTaq Green Master mix (Promega), with 1 pmol of biotinylated forward primer, 1 pmol of reverse primer and 50 ng of bisulfite treated genomic DNA. The primer sequences for TABPB were: forward—5’- GTGGATTTTGGGTGGGGATTA-3’ and reverse—5’- AACTCCCACCAAACCATCCTTACC-3’; for ZBTB18 were: forward—5’- GAGTTTGAGGAGATGTATTTGATATT-3’ and reverse—5’- AACTTTTCAACCAATTTATAAATCTTTTCT-3’. The PCR was followed by pyrosequencing, using conditions previously described [33]. Results were reported as the percentage of the 5-methylated cytosines with respect to total, methylated and unmethylated, cytosines. Additionally, non-CpG cytosine residues were used as internal controls to validate bisulfite conversion. 2.5. Ingenuity pathway analysis Ingenuity Pathway Analysis software (IPA, Qiagen) was used to annotate genes with differentially methylated promoters (DMPs) and rank associated canonical pathways. Canonical pathways analysis identified the pathways from the IPA library of pathways that were most significant to the set of genes with DMPs. The significance of the association between a gene set and a canonical pathway was measured by determining the ratio between the number of genes in our list that mapped to the canonical pathway divided by the total number of all known genes that mapped to that pathway. Fisher’s exact test was used to calculate p‐values determining the probability that the association between the genes in the dataset and the canonical pathway was explained by chance alone.

4. Discussion DNA methylation plays an important role in the interplay between external (environment) and internal (gene expression) factors and thus may explain the late or stimulus-triggered on-set of multisystem complex diseases such as ME/CFS. However, only a few studies have reported changes in DNA methylation associated to ME/CFS [18–20]. Here, we used a newer and improved, technology developed by Illumina Inc. that was not available at the time the preceding studies were performed, with the intention of expanding previous findings. Illumina MethylationEPIC microarrays have almost twice the coverage of CpG sites compared to the older Illumina HumanMethylation450 microarrays, including 333,265 additional CpGs located in regulatory regions as identified by the ENCODE and FANTOM5 projects [36]. Here, we identified differentially methylated CpG sites and promoters in PBMCs of ME/CFS cases compared to controls in cohorts from two distant geographic locations. Association of hypomethylation of DMPs with immune regulatory pathway is in agreement with previous studies [18–20]. Those previous studies found immune cell regulation as the largest coordinated enrichment of differentially methylated pathways and reported genomic DNA hypomethylation of genes in immune pathways from PBMCs isolated from ME/CFS cases [18, 19] and also in CD4+ T cells from ME/CFS cases [20]. Specifically, we found significant differences in DNA methylation between ME/CFS cases and controls at 17,296 CpG sites in 6,368 genes overall (FDR ≤ 0.05), across promoters as well as other gene regulatory elements and within coding regions of genes. DNA hypomethylation was found in 98% of these sites in ME/CFS, whereas only 2% were hypermethylated compared to controls. Differential methylation of promoters can have an effect on the expression of the corresponding affected genes. Interestingly enough, less than half of the DMPs (47.6%) associated with protein-coding genes while over 50% affected non-protein coding DMPs corresponding to promoters of regulatory gene elements including short non-coding RNAs (including miRNAs, small nuclear RNAs, small nucleolar RNAs) (10.1%), long non-coding RNAs (including antisense RNAs and intergenic RNAs) (30.6%), and pseudogenes (11.4%) (Fig 5). Small nuclear RNAs and small nucleolar RNAs participate in RNA splicing [47] while long non-coding RNAs and pseudogenes play important roles in the regulation of transcription and mRNA stability by annealing to transcripts and sequestering miRNAs as molecular sponges [48], among other epigenetic mechanisms [49]. It is important to note that the highest beta difference (hypomethylation at 14%) in the DMP list was observed for the promoter of the gene that is related to natural killer (NK) cells function (S6 and S7 Tables). At the single gene level, our results showed hypomethylation of promoters of genes related to those of reference studies such as some solute carriers, RNA binding proteins, homeobox, and PHD finger proteins [20]. Among them we found the MED13L gene which is associated with muscular hypotonia and neurocognitive impairment [50], symptoms that associate with ME/CFS. We also found DMPs linked to the expression of ribosomal proteins (RPL23A and RPL7) as well as to the 5S ribosomal RNAs pseudogenes (RNA5SP245, RNA5SP77 and RNA5SP97) (S5 Table) suggesting possible disturbance of protein synthesis. In regards to the methylation status of glucocorticoid response-associated genes, we found significant hypomethylation of the SGK3 (serum glucocorticoid regulated kinase) gene promoter [18, 19]. Transcriptional profiling studies in ME/CFS cases have pointed out perturbations in T-cell and B-cell activation and dysregulation in gene expression broadly related to immune responses [51–56], and such changes fit with other studies showing altered production of interleukin and interferon cytokines as features of ME/CFS immune dysfunction [57]. In the current study, we report hypomethylation of DMPs of genes known to regulate the adaptive immune response such as immunoglobulins and pseudogenes (for example, IGDCC3, IGHD7-27, IGHJ1, IGHJ3 and other) (S6 and S7 Tables) which is in agreement with the reported improvement of ME/CFS cases following B cell depletion therapy [58, 59]. We also observed hypomethylation of promoters of MMP14, MAP4K4, MAPK12 and CREB5 (S6 and S7 Tables) possibly activating the TNF signaling pathway, that fits with the reported over-expression of pro-inflammatory cytokines [60] in ME/CFS. In addition, we also show hypomethylation of miRNA-148a promoter that could potentially contribute to its overexpression. It is known that the miRNA-148-152 family can impair the innate immune response and antigen presentation of TLR-triggered dendritic cells [61]. Furthermore, miRNA-148a can contribute to DNA hypomethylation in lupus CD4+ T cells by repressing DNA methyltransferase 1 expression [62]. It is noteworthy to mention that amongst the DMPs associated to miRNA expression in this study (S5 Table), miRNA-10a, miRNA-181 and miRNA-4710 were also reported to be affected at the DNA methylation level in at least one previous ME/CFS study [20]. Specifically, miRNA-10a has been associated with ME/CFS in NK cells [63] and miR-181 is related to TLR-mediated inflammatory response [64], which is associated with ME/CFS [65]. Hypomethylation of the IL21R gene promoter (S6 and S7 Tables) may indicate increased IL21 receptor expression. The ligand binding of this receptor leads to the activation of multiple downstream signaling molecules, including JAK1, JAK3, STAT1, and STAT3. Knockout studies of the ortholog mouse counterpart suggest a role for this gene in regulating immunoglobulin production [66]. IL21 regulates both innate and adaptive immune responses, and also exerts major effects on inflammatory responses that promote the development of autoimmune diseases and inflammatory disorders [67]. Importantly, IL21 signaling is critical for induction of spontaneous experimental autoimmune encephalomyelitis [68]. The DMPs reported in our studies are not only consistent and validate the previous studies of gene regulatory elements in genes within the immune cell regulation cluster [14–20], but also provide an improved understanding using advanced technology as well as validation using cohorts from distant geographic locations. It is noteworthy to mention that it is still inconclusive from previous studies or our current study, whether the significant epigenetic modifications found to associate with ME/CFS indicate a compensatory homeostatic mechanism or result from an adaptive immune response towards environmental inducers. However, these results indicate that DNA methylation constitutes a potential gene regulatory mechanism capable of mediating long-term changes in ME/CFS cases, as previously noted by other authors [18–20]. In summary, the association of differential DNA methylation with ME/CFS definitely suggests a potential role for epigenetic alterations in the pathophysiology of ME/CFS.

5. Conclusions To our knowledge, this is the first study that has explored genome-wide epigenetic changes associated with ME/CFS using the new Illumina MethylationEPIC microarrays covering over 850,000 CpG sites and that has validated the results in two geographically distant cohorts of ME/CFS cases. Our findings build on previous preliminary reports showing association of altered methylation profiles of genes with immune functions and confirm DNA methylation as an epigenetic regulatory mechanism associated with ME/CFS. Future validation studies in larger cohorts of ME/CFS cases are needed to confirm these findings and evaluate the effects of such methylation patterns on gene expression.

Acknowledgments We acknowledge employees of the Center for Genome Technology at the John P. Hussman Institute for Human Genomics, University of Miami Miller School of Medicine, Ioanna Konidari, Daniela Martinez and Martin Magurno for running Illumina MethylationEPIC microarrays. We would like to thank Beth Gilbert, M.S., for the help in preparation of manuscript.