Identification of Myalgic Encephalomyelitis/Chronic Fatigue Syndrome-associated DNA methylation patterns

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.


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][3][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][9][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][14][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][19][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.

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].
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 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. , 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.

ME/CFS Cases
The Student t-test was used to assess statistical difference between ME/CFS cases and controls in SF-36 scores for physical and mental health.

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 10 7 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).

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

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 Methylatio-nEPIC 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, sitebased 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 microarraybased 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.

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'-AACTCCCACCAAAC CATCCTTACC-3'; for ZBTB18 were: forward-5'-GAGTTTGAGGAGATGTATTTGATA TT-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.

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.

Study design and participant demographics
Flow chart of the study design is depicted on Fig 1. In brief, genomic DNA from 13 ME/CFS cases and 12 controls ("experimental cohort") was hybridized to the Illumina MethylationEPIC Genomic DNA was isolated from ME/CFS cases (N = 13) and controls (N = 12) and analyzed using Illumina Infinium Human MethylationEPIC BeadChip microarrays and this was defined as "experimental cohort". Following data analysis results were validated on all the samples from the "experimental", as well as from a secondary cohort that included samples from ME/CFS cases (N = 8) and controls (N = 8) from Miami/Fort Lauderdale area, as well as ME/CFS cases (N = 12) and controls (N = 11) from Valencia (Spain), using pyrosequencing analysis.
https://doi.org/10.1371/journal.pone.0201066.g001 microarrays, as detailed in Methods. These arrays can assay over 850,000 CpG sites and cover different regions across genes, including promoter, 5' UTR, first exon, gene body, and 3' UTR, as well as and distant regions. Coverage includes 99% of RefSeq genes with multiple probes per gene, 96% of CpG islands from the UCSC (University of California, Santa Cruz) database and additional content imputed from whole-genome bisulfite sequencing data [23]. Following data analysis, the obtained results were validated by pyrosequencing on the samples from the experimental cohort, as well as on additional 8 ME/CFS cases and 8 controls from Miami/Fort Lauderdale area and 12 ME/CFS cases and 11 controls from Valencia (Spain), as described in Methods.
It is worth noting that previous studies have had similar number of cases and controls using different Illumina Human Methylation450 microarrays. For example, Brenu et al [20] included 25 ME/CFS cases with 18 controls and measured methylation in T-cells; whereas De Vega et al [18] included 12 ME/CFS cases and 12 controls and measured methylation status in PBMCs. Hence, our sample size is similar to these previous publications for statistical comparison. Moreover, the major strength of the current findings is depicted in the validation of specific CpG sites in an expanded cohort of ME/CFS cases and controls from a different geographical location i.e. Valencia, Spain. However, additional validation in a larger cohort is warranted. Table 1 shows age and BMI matching of all ME/CFS cases and controls in the experimental cohort and in the 39 additional samples from two distant geographic locations (used for pyrosequencing validation). Vitality, social functioning, role emotional and mental health were significantly lower in ME/CFS cases compared to controls. In addition, statistically significant differences in physical health were found when cases were compared to controls, according to SF-36 scores. S1 and S2 Tables show age, BMI and SF-36 scores for the two separate cohorts, one from Florida, USA, and the other from Valencia, Spain, correspondingly.

Global methylation in ME/CFS
We did not find statistically significant differences in global DNA methylation levels using ELISA-based assay methods to detect either 5-methylcytosine (5-mC) or 5-hydroxymethylcytosine (5-hmC) modifications on DNA isolated from PBMCs of ME/CFS cases and controls (S3 Table). This analysis, however, shows overall cytosine modifications but it is limited at detecting specific fluctuations at particular locations in the genome.
Since former studies had shown site specific methylation differences in ME/CFS [18][19][20][21], we proceeded to evaluate this possibility in our experimental cohort, however, we used an upgraded version of the microarrays for human genome-wide methylation studies (Illumina Infinium Human MethylationEPIC microarrays). This novel technology which has almost twice coverage compared to that of previous Illumina HumanMethylation450 microarrays [36] may allow for validation and additions to previous findings in independent, geographically distant, cohorts of ME/CFS cases.

Differentially methylated CpG sites in ME/CFS
Data from the microarray-based genome-wide methylation analysis of the experimental cohort showed hypomethylation of the differentially methylated sites (DMS) (98%) in PBMCs from ME/CFS cases compared to controls.
As many as 649,836 probes remained for normalization and further analysis upon filtering and preprocessing of the raw genomic DNA methylation data with the R/Bioconductor v.3.3.1 RnBeads package v.1.2.0. Analysis of these probes led to the finding of 17,296 DMS located in 6,368 annotated genes indicating epigenetic differences between the ME/CFS cases and controls in experimental cohort when we used FDR 0.05 criteria. S4 Table shows the list of DMS along with their respective gene references. Furthermore, Principal Component Analysis (PCA) of full methylomes clearly differentiated between ME/CFS cases and controls as shown in Fig 2. In addition, unsupervised hierarchical clustering of the 500 most variable CpG sites (according to RnBeads ranking, see Methods) shows clear segregation between ME/CFS cases and controls (Fig 3) confirming that methylation DNA patterns of PBMCs in ME/CFS individuals differ from those of controls.
Out of the 17,296 DMS, a total of 14,261 DMS (82%) were found within or proximal to genes (i.e., genic locations) indicating that they may have a role in regulation of gene expression. About 98% of the DMS within the genic regions were hypomethylated and 2% were hypermethylated overall (Fig 4A and S4 Table). These DMS were mostly equally distributed between gene bodies (43.68%) and gene regulatory elements (combined TSS1500, TSS200, 5' UTR and 3' UTR) (43.84%). Hypomethylation of ME/CFS DNA was observed on about 21.87% of DMS located in TSS regions; 17.21% in 5' UTRs; 3.08% in 1st exons and 4.25% in 3' UTRs. Hypermethylation, on the other hand, appeared on 0.16% of DMS located in TSS regions, 0.34% in 5 and 3' UTRs, and 0.79% in gene bodies (Fig 4A). Regarding to the location relation to CpG islands, out of 17,296 DMS, 46.8% were located in CpG islands, whereas 53.2% (8,089 DMS) belong to proximity of CpG islands and CpG islands. As shown in Fig 4B, the mapping of 8,089 DMS with respect to their location in CpG islands shows that 99.3% of these 8,089 DMS are hypomethylated and only 0.7% were hypermethylated. Fig 4B also shows that out of these 8,089 DMS, approximately 25.25% were located within CpG islands, 61.9% mapped 2 kb upstream and downstream of CpG islands (N, S Shores), and 12.8% of them localized to regions 2 kb upstream and downstream of CpG shores (N, S Shelves). Within these regions, 74.1% of shores and shelves were hypomethylated, while only 0.63% of CpG island shores and shelves were hypermethylated, consistent with the general finding of hypomethylation being associated to ME/CFS. More relaxed criteria (FDR 0.1 and mean beta differences > 0.05) allowed to identify 27,129 DMS as listed in S5 Table.

Differentially methylated promoters in ME/CFS
Since methylation of promoter regions directly affect levels of gene expression, we proceeded to analyze differentially methylated promoters (DMPs) in our experimental cohort. As described in Methods, promoters were defined in RnBeads package as regions located between 1,500 bp upstream of TSS and 500 bp downstream of TSS. Criteria for DMPs were FDR 0.1 with the mean beta differences > 0.05. Using these criteria, we found 307 DMPs in PBMCs of ME/CFS cases. S6 Table lists these DMPs ranked according to the RnBeads package [35]. 306 of these DMPs are hypomethylated in the ME/CFS group while only one appeared hypermethylated. Fig 5 shows the distribution of these promoters according to their gene biotype. Almost half of the identified DMPs belong to protein-coding genes (47.6%), and more than half (52.1%) belong to regulatory RNA genes. The last consisted of 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%). This distribution of DMPs supports our hypothesis that the differential methylation observed might lead to deregulation of gene expression in ME/CFS. We also identified 144 DMPs with more stringent criteria of FDR 0.05 with the mean beta differences > 0.05 as listed in S7 Table; and it represents similar gene-function patterns as DMPs identified with FDR 0.1.

Functional pathway analysis of DMPs
To further understand the potential functions and interactions between the 307 genes controlled by DMPs within the context of cell biological pathways, we performed Gene Ontology  N = 12). 500 DMS above and below mean level are represented by blue and yellow across all cases and were used to cluster cases. As shown above, 12 of 13 ME/CFS cases (blue on color-coded bar above the dendrogram) cluster together (left dendrogram branch) and 12 of the 12 controls (blue color-coded bar) cluster together (right dendrogram branch), resulting in a divergence of these sub-phenotypes.

Validation of Illumina MethylationEPIC microarray results by pyrosequencing
To validate the results of our analysis obtained with Illumina MethylationEPIC microarrays, we selected two DMPs that had high density of CpG sites and possible ME/CFS relevance by their role in the neuro-immune signaling axis. The first, corresponding to the Zinc Finger and BTB Domain Containing 18 (ZBTB18), is a zinc finger transcriptional repressor with a crucial role in brain development and neuronal differentiation. The second, TABPB, encodes a transmembrane glycoprotein which mediates interaction between newly assembled major histocompatibility complex (MHC) class I molecules and the transporter associated with antigen processing (TAP), which is required for the transport of antigenic peptides across the endoplasmic reticulum. ZBTB18 is in the upper top list of DMPs by order of RnBeads combined rank (position 8 of 307 in S5 Table); while TABPB is in the lower part (position 246 of 307 in S5 Table), thus we included DMPs with values at both ends of the range in the validation step. For each of the 2 selected genes, we assayed the level of DNA methylation of six CpG sites in the promoter area. For ZBTB18, three of the CpG sites corresponded to oligonucleotides printed on Illumina MethylationEPIC microarrays, and three additional CpG sites not included in microarray content but located in the gene promoter region between cg15896892 and cg19698993, were also tested. In the case of the TABPB gene, all the assayed CpG sites correlated with probes included in Illumina MethylationEPIC microarrays. As summarized in Fig 7, pyrosequencing analysis strongly validated the direction of methylation differences between ME/CFS and control subjects (p < 0.01) as identified by Illumina MethylationEPIC microarray (FDR 0.05) in all CpG sites assayed. Importantly, it was observed that the methylation levels detected by pyrosequencing in ME/CFS individuals were similar to those observed by MethylationEPIC microarrays in all DMS (Fig 7), and that the differences in methylation levels assayed by pyrosequencing were even more pronounced than the same differences assayed by microarrays (Fig 7).

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 Average methylation level of CpG sites within promoters of TABPB and ZBTB18 were assessed using pyrosequencing. Corresponding probes for TABPB: CpG1-cg04415168 (TSS1500), CpG2-cg10376053 (TSS1500), CpG3-cg14288848 (TSS1500), CpG4-cg17055704 (TSS1500), CpG5-cg14473643 (TSS1500), CpG6-cg14309283 (TSS1500). Corresponding probes for ZBTB18: CpG1-cg16399365 (TSS1500), CpG2 -cg15896892 (TSS1500), CpG6-cg19698993 (TSS1500). CpG3, CpG4 and CpG5 are not printed on the Illumina MethylationEPIC microarrays and are located in ZBTB18 promoter between CpG2 and CpG6. Ã = FDR 0.05 for EPIC; Ã = p 0.05 for PS. Error bars represent the standard error of the mean. 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][19][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][19][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 miR-NAs, 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][52][53][54][55][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][15][16][17][18][19][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][19][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.

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.  (12)