DNA Methylation Changes in Fibromyalgia Suggest the Role of the Immune-Inflammatory Response and Central Sensitization

Fibromyalgia (FM) has been explained as a result of gene-environment interactions. The present study aims to verify DNA methylation differences in eleven candidate genome regions previously associated to FM, evaluating DNA methylation patterns as potential disease biomarkers. DNA methylation was analyzed through bisulfite sequencing, comparing 42 FM women and their 42 healthy sisters. The associations between the level of methylation in these regions were further explored through a network analysis. Lastly, a logistic regression model investigated the regions potentially associated with FM, when controlling for sociodemographic variables and depressive symptoms. The analysis highlighted significant differences in the GCSAML region methylation between patients and controls. Moreover, seventeen single CpGs, belonging to other genes, were significantly different, however, only one cytosine related to GCSAML survived the correction for multiple comparisons. The network structure of methylation sites was different for each group; GRM2 methylation represented a central node only for FM patients. Logistic regression revealed that depressive symptoms and DNA methylation in the GRM2 region were significantly associated with FM risk. Our study encourages better exploration of GCSAML and GRM2 functions and their possible role in FM affecting immune, inflammatory response, and central sensitization of pain.


Introduction
Fibromyalgia (FM) is a long-term pain syndrome characterized by chronic widespread pain (CWP) and a constellation of additional comorbidities, mainly including fatigue, sleep impairment, depression, and cognitive dysfunctions. Its prevalence is 2-4% in the general population, with higher frequencies in women than in men [1]. Thanks to the development of diagnostic criteria [2], and the recently approved ICD-11 (International Classification of Diseases 11th Revision) coding system [3], the diagnosis of FM has improved in the last years. However, one of the main problems remains the lack of objective markers. 2 of 15 The pathogenesis of this disease is unclear. Central sensitization, with alterations in nociceptors, neurons, and glia processing pain signals, was proposed to explain CWP [4]; consistently, FM has recently been characterized as nociplastic pain according to a mechanistic description [5]. Other hypotheses for FM pathogenesis include altered inflammatory mechanisms and altered immune system response [6,7], abnormalities of hypothalamicpituitary-adrenal (HPA) axis [8], or altered dopamine response to pain [9].
A gene-environment interaction, predisposing or protecting against FM risk, model seems to be the best explanation for these pathological phenotypes. The combinations of polymorphisms in the serotoninergic and catecholaminergic pathways [10,11] and environmental factors [12] might be mediated or affected by epigenetic changes. Epigenetics includes heritable changes in the gene function that cannot be explained by changes in the DNA sequence, influencing both gene expression and phenotype [13,14]. In particular, DNA methylation is an epigenetic mark involved in gene expression regulation, catalyzed by a family of DNA methyltransferases that transfer a methyl group from S-adenyl methionine onto the DNA cytosine to form 5-methylcytosine [15]. Studies proposed that DNA methylation might reflect or contribute to the complex gene-environment interplay inducing FM pathogenesis [11]. However, the processes by which epigenetics might fine-tune the relationship between the genetic background, experienced stress and the development of FM remain a major challenge, particularly in humans. The few studies evaluating DNA methylation as a potential biomarker of FM [11] often considered a heterogeneous population with incomplete phenotypic descriptions, and no attention for comorbidities that might highly impact the epigenetic signatures.
Therefore, the present study aims to explore DNA methylation in a group of FM women and their healthy sisters, subjected to extensive clinical evaluation. We used a twofold approach. First, we compared the two groups on the methylation sites independently. Second, to better understand the complex disease state, we compared siblings' methylation using a network-based modeling, an undirected graph holding vertices as methylation sites and edges as the association between them. Because sociodemographic and clinical variables might affect DNA methylation, we also tested the simultaneous influence of DNA methylation levels and these variables on FM development. Exploring epigenetic variations could reveal new insights in FM pathogenesis or biomarkers of the disease.

Subjects
Eighty-four Caucasian participants (42 female patients with FM and 42 related healthy sisters) were selected for the present study from a cohort of 543 families in which at least one member was diagnosed with FM. The more severely FM affected patients (group 5, previously described in a pilot study [16]) were recruited and diagnosed by a primary care physician or by a professional specialist in rheumatology or neurology, according to the ACR (American College of Rheumatology) 2010 criteria. For each FM patient, a related healthy sister was selected as control. The study design is shown in Figure S1.

Demographic and Clinical Assessment
All the participants, subjects and healthy controls (HCs), were submitted to a clinical interview about demographic data and to the following scales and questionnaires (Table S1): Fibromyalgia Impact Questionnaire (FIQ) [17,18] and Visual Analog Scales (VAS) to assess the core symptoms of FM [19]; Beck Depression Inventory (BDI) [20,21] and Pittsburgh Sleep Quality Inventory (PSQI) [22,23] to assess, respectively, depressive symptoms and sleep disturbance.

Samples Collection
As previously reported [16], peripheral whole blood collection, two tubes of 10 mL per subject, was performed via venipuncture and leukocytes were separated through a washing protocol. DNA purification protocol from the isolated leukocytes was performed using QIAamp DNA Blood Midi/Maxi Kit (Spin Protocol, QIAGEN) at the Galician Public Foundation of Genomic Medicine of the University of Santiago de Compostela (Spain). Aliquots of the genomic DNA extracted were sent to Aalborg University for the present epigenetic study.

DNA Methylation Analysis
We analyzed DNA methylation level in eleven genomic regions (Table 1), including CpGs islands, promoters, and transcription start sites. These regions of interest were representative of the main symptoms and hypotheses associated to FM pathogenesis (pain perception, inflammatory response, stress and immune system, dopaminergic pathway), as identified in a previous pilot study [16]. 27325757 [37] Targeted NextGen Bisulfite Sequencing was conducted by EpigenDx, Inc. (Hopkinton, MA, USA) in four main steps. (i) Bisulfite conversion: extracted DNA samples (500 ng) were bisulfite modified using Zymo EZ-96 DNA Methylation™ Kit (Zymoresearch, CA, USA); (ii) PCR amplification: the bisulfite-treated DNA were subsequently amplified with separate multiplex or simplex PCRs using Qiagen HotStarTaq (0.5 units), 0.2 µM primers, and 3 µL of in a final volume of 20 µL. Quality and quantity of the PCR products were checked using the QIAxcel Advanced System, and then they were pooled and purified using QIAquick PCR Purification Kit columns (Qiagen); (iii) libraries preparation: a custom Library Preparation method created by EpigenDx was used and then library molecules were purified using Agencourt AMPure XP beads (Beckman Coulter, Brea, CA, USA) and quantified using the Qiagen QIAxcel Advanced System. Barcoded samples were then pooled in an equimolar fashion before template preparation and enrichment were performed on the Ion Chef™ system (Thermo Fisher, Waltham, MA, USA) using Ion 520™ and Ion 530™ ExT Chef reagents; (iv) sequencing: following this, enriched, templatepositive library molecules were then sequenced on the Ion S5™ sequencer using an Ion 530™ sequencing chip (Thermo Fisher, Waltham, MA, USA).

Bioinformatic Data Processing
Bisulfite treated single-end sequencing reads were mapped to the human reference genome (GRCh38/hg38), using BS-Seeker2 [38] with default parameters and bowtie2 aligner [39]. Conversion of the alignment files to the CGmap format has been performed using CGmapTools [40]; finally, Metilene [41] was used for the identification of differentially methylated regions and cytosines.

Statistical Analyses
Metilene software, with a binary segmentation algorithm combined with a twodimensional statistical test [41], was used for the detection of differentially methylated regions (DMRs, test's parameters: -f 2 -m 1 -d 0.01) using the whole-genome methylation matrix and supplying the genomic coordinates of our regions of interest. To call differentially methylated cytosines (DMCs), Metilene was supplied with a subset of the whole-genome methylation matrix containing only our regions of interest (test parameters: -f 3 -m 1 -d 0.01). In particular, the software assesses the statistical significance of potential DMRs by a two-dimensional version of the Kolmogorov-Smirnov test (KS-test) [42] and an independent Mann-Whitney U test (MWU-test). The software Metilene assesses each cytosine for differential methylation (DMCs test) by using the Mann-Whitney-U test. The corresponding p-values are reported in the output.
The associations between the different methylated genome regions have been explored through a network analysis. This approach aims to clarify the structure of a large number of connections that otherwise would seem irrelevant or too complicated. Specifically, being interested in the differences between the two groups (i.e., FM patients and HCs), two networks were built using partial correlations. A network is a graphical representation of the relationships (named edges) between variables (named nodes). Each network was first compared by the centrality measures and then the two networks' structures were compared using the edge correlations. This approach defines in each network the magnitude of association between variables (edge weights). A significant correlation between these two sets of variables suggests a similarity between the two networks. The networks have been plotted with JASP (v 0.11.1.0), whereas the comparison analyses were performed using R (v 3.6.2).
Concerning the socio-demographic data and other information such as medication, the two samples were analyzed using T-tests or chi-squared tests, depending on whether they were continuous or categorical variables. The descriptive statistics for all the variables included mean, standard deviation, frequencies, and percentages.
Subsequently, a logistic regression model was estimated to test the concurrent effect of mean methylation levels in the candidate genome regions on the risk to develop FM. As the presence of depression in FM patients is a potential confounding factor, depression scores measured through BDI were also included in the model. In addition, because socio demographic variables might affect the methylation levels, the model also controlled for the following socio-demographic variables: age in years, BMI in Kg/m 2 ; job status (coded as 0 = unemployed, 1 = employed); education (coded as 0 = none, 1 = primary, 2 = secondary, 3 = professional and 4 = university); living status (coded as 0 = alone and 1 = with someone). Additionally, current use of the following medications was coded as 0 = none and 1 = at least one: GABAergic medications (such as benzodiazepines), gastrointestinal medications (such as proton pump inhibitors), cardio-vascular medications (such as statins), and antihistaminic medications. Antidepressant medications and painkillers were not included because they represent the gold-standard treatment for FM and hence would have perfectly split the two samples. Since the sample size did not allow the inclusion of all the co-variates, a two-step Cluster Analysis (including both continuous and categorical variables) was applied to reduce the dimensionality of the covariates (IBM SPSS 26.0). The distance between the variables was the log-likelihood. The number of clusters extracted was automatic and based on the Bayesian Information Criterion (BIC).
The clusters have been thus added in the regression model with the candidate genome regions and depression.
Robust standard errors were applied to the regression models in order to reduce the possible bias introduced in the estimations by heteroscedasticity. The regression analyses were conducted in Stata/IC 15.1 (StataCorp, College Station, TX USA).
For all the statistical analyses, results were considered statistically significant for p ≤ 0.05.

Characteristics of the Study Participants
The forty-two FM women were aged 22-75 years (mean age 50 ± 10 years) and the 42 related healthy sisters were aged 28-72 years (mean age 47 ± 10 years). Table 2 shows the patients' characteristics. Questionnaires' and scales' scores related to depression (BDI), sleep impairment (PSQI), and the main FM symptoms of pain (FIQ, WPI, SSS, VAS) were all significantly higher in patients than controls (p < 0.000). Among the other variables, the number of participants consuming anti-anxiety, antidepressant, anti-inflammatory medications, and opioids painkillers were significantly higher in FM women compared with controls (p < 0.000). The average time since diagnosis was about 9 years, reflecting the long-term course of FM management.

DNA Methylation Analysis Comparing Cases and Controls
Two types of analyses were applied: the differentially methylated regions (DMRs) test to determine methylation by grouping neighboring cytosines, and the differentially methylated cytosines (DMCs) test to reveal methylation at single cytosine level (Table S2).
Concerning DMRs (Table 3), only one differentially methylated region (methylation difference ≥ 1%), evidenced by Metilene software, showed significant differences by using Mann-Whitney U-test (MWU-test) and the Kolmogorov-Smirnov test (KS-test): the region ( Figure S2) was related to GCSAML gene, the Germinal center associated signaling and motility like gene, a putative signaling protein with a potential function in the immune response. The GCSAML region (chr1: 247518380-247518621) showed a level of methylation significantly higher in FM patients (mean methylation 0.168) than controls (mean methylation 0.158) with both tests (p (MWU) = 0; p (KS) = 0.007). Concerning DMCs, differences in the level of methylation in the single cytosines of the included regions comparing FM women with their healthy sisters were verified by using Mann-Whitney U-tests. The corresponding p-values reported in the output (Table 4) showed seventeen differentially methylated cytosines, belonging to six genes, GCSAML, DRD3, TRPA1, IL25, OXT, and MCF2, that reached statistical significance (p < 0.05). In particular, five cytosines were evidenced in the GCSAML gene (methylation difference ≥ 5.06%), confirming once again the possible correlation of this gene with FM. Three DMCs were related to the region of the TRPA1 gene, encoding a receptor involved in pain detection (methylation difference ≥ 1.91%). Another five cytosines were evidenced in the OXT gene (methylation difference ≥ 1.8%), the oxytocin hormone, involved in stress, cognition, and complex behavior. Two DMCs resulted in the MCF2 gene region (methylation difference ≥ 4.8%); this gene encodes an oncogenic protein, a member of the DBL family of Rho GEFs (Rho GDP-GTP exchange factors) that exerts control over some members of the Rho family small GTPases. Finally, one differentially methylated cytosine was identified in the DRD3 gene (chr3: 114178637-114178638; methylation difference = 3.6%), and one in the IL25 region (chr14: 23372248-23372249; methylation difference = 3.8%). Thirteen out of these seventeen identified cytosines were higher methylated in FM women compared with their healthy sisters. Four cytosines (GCSAML, chr1: 247518586-247518587; TRPA1, chr8: 72076406-72076407; OXT, chr20: 3071336-3071337/chr20: 3071468-3071469) resulted in lower methylated in FM women. However, applying the correction for multiple comparisons, the significance of the DMRs and DMCs disappeared, except for one cytosine (chr1: 247518586-247518587) related to the GCSAML gene.

Network Analysis
A visual inspection evidenced that the two extracted networks (Figure 1a) showed different centrality measures (Figure 1b): for the patients, DNA methylation level in the GRM2 region represented the most central and connected node, whereas for the HCs group, the most central node is the MCF2 methylation level. The analysis of correlation between edges showed a non-significant correlation between the networks (r = 0.1169553; p = 0.3940608), confirming a different structure for each sample. Comparison of the centrality plots of the two networks in FM patients (coded as 1, light blue) and their healthy sisters (coded as 0, orange). The betweenness quantifies the number of times a node acts as a bridge along the shortest path between two other nodes. The closeness of a node is the average length of the shortest path between the node and all other nodes in the graph. The degree is defined as the number of links incident upon a node (i.e., the number of ties that a node has). On the y axis the methylation sites are reported, and on the x axis the strength of these three dimensions (betweeness, closeness, and degree) is evidenced. Higher scores suggest a greater importance of the methylation site in the network structure. As confirmation of the inspection level, GRM2 had greater levels of betweeness, closeness, and degree, having the highest value among all the methylation sites in the FM group only. This was not true for the healthy sisters, in which SYT2 and MCF2 were more central to the network.

Cluster Analysis of the Socio-Demographic Variables
The cluster analysis extracted two clusters with a sufficient quality and balance (ratio between the two clusters: 2.23). The main three predictors of the clusters were employment, education, and GABAergic medications. The first cluster (n = 22, 31%) was constituted mainly by subjects that were employed (95.5%), had a university education (40.9%), were not taking benzodiazepines (95.5%) nor other medications, and were slightly younger (mean age 43.64 years). On the other hand, the subjects in the second cluster (n = 49; 69%) were older (mean age 52.12), mostly had only primary education (71.4%), were housewives or unemployed (67.3%), and nearly half were taking benzodiazepines (53.1%). The two clusters were almost overlapping regarding living status and BMI. When comparing the two clusters on the diagnosis of FM, these were not significantly different (chi-squared = 3.896; p = 0.72). Comparison of the centrality plots of the two networks in FM patients (coded as 1, light blue) and their healthy sisters (coded as 0, orange). The betweenness quantifies the number of times a node acts as a bridge along the shortest path between two other nodes. The closeness of a node is the average length of the shortest path between the node and all other nodes in the graph. The degree is defined as the number of links incident upon a node (i.e., the number of ties that a node has). On the y axis the methylation sites are reported, and on the x axis the strength of these three dimensions (betweeness, closeness, and degree) is evidenced. Higher scores suggest a greater importance of the methylation site in the network structure. As confirmation of the inspection level, GRM2 had greater levels of betweeness, closeness, and degree, having the highest value among all the methylation sites in the FM group only. This was not true for the healthy sisters, in which SYT2 and MCF2 were more central to the network.

Cluster Analysis of the Socio-Demographic Variables
The cluster analysis extracted two clusters with a sufficient quality and balance (ratio between the two clusters: 2.23). The main three predictors of the clusters were employment, education, and GABAergic medications. The first cluster (n = 22, 31%) was constituted mainly by subjects that were employed (95.5%), had a university education (40.9%), were not taking benzodiazepines (95.5%) nor other medications, and were slightly younger (mean age 43.64 years). On the other hand, the subjects in the second cluster (n = 49; 69%) were older (mean age 52.12), mostly had only primary education (71.4%), were housewives or unemployed (67.3%), and nearly half were taking benzodiazepines (53.1%). The two clusters were almost overlapping regarding living status and BMI. When comparing the two clusters on the diagnosis of FM, these were not significantly different (chi-squared = 3.896; p = 0.72).

The Concurrent Effect of DNA Methylation, Depression, and the Clustered Sociodemographic Data on FM Risk
We performed a binary logistic regression to predict the risk to develop FM, including in the model as independent variables the mean methylation levels of the candidate genome regions, depression, and the sociodemographic clusters. The diagnosis of FM was included as a dependent variable (coded as 0 = no and 1 = yes) ( Table 5). Among the methylated regions inserted in the model, only the level of methylation in GRM2 gene was significantly and negatively associated with FM diagnosis; in particular, the DNA methylation increase in the GRM2 region (chr3: 51706813-51707270) resulted to confer 39% lower risk to develop FM (OR = 0.614; 95% CI = 0.388-0.971; p = 0.037). Among the other variables entered in the model, only depressive symptoms remained associated with FM: a unit increase in BDI scale (used to measure depression) corresponded to 1.3 times higher risk of suffering from FM (OR = 1366; 95% CI = 1.170-1.594; p < 0.000). It is important to note that the variables inserted in the present model resulted to account for about 56% of the total variability of the dependent variable FM.

Discussion
The present study analyzed the DNA methylation levels in eleven genome regions of 42 women with fibromyalgia compared with their 42 healthy sisters, using a siblings approach that reduces the genetic heterogeneity and differential prenatal or early-life exposures: this represents a powerful design to investigate the association of DNA methylation with FM. The comparison revealed a slight but significant different methylation level in the GCSAML gene region and identified differentially methylated cytosines in GCSAML, DRD3, TRPA1, IL25, OXT, and MCF2 genes regions. Nevertheless, except for one cytosine in the GCSAML gene, all the differences disappeared, applying the correction for multiple comparison. The networks analysis revealed a significantly different structure of methylation sites comparing the two groups, with GRM2 methylation representing a central node only in the FM group. When testing the simultaneous effects of the mean methylation levels in the candidate genome regions together with depression and the clustered sociodemographic clinical data, GRM2 methylation was significantly and negatively associated with FM risk, while depression was positively associated with it.

DNA Methylation Analysis
Our findings in the differentially methylated regions (DMRs) test brings the focus on the GCSAML gene, with an increased methylation level observed in FM women compared with their healthy sisters. The GCSAML gene, which was observed to be subjected to epigenetic regulation, encodes a signaling molecule thought to be associated to the sites of proliferation and differentiation of mature B lymphocytes [28]. This result points out the connection with potential immune system dysfunction in FM patients. Previous investigations revealed altered expression of immune pathways and markers of tissue destruction in FM women [30,43,44]. In chronic fatigue syndrome, differentially methylated [45] and differentially expressed genes related to the immune response [46,47] were identified. A recent study also showed specific transposable elements overexpressed in the immune cells of FM patients [48]. Taken as a whole, these studies, together with our finding, appear to support a possible immune system dysregulation in FM. However, most of the epigenetic studies on myalgic encephalomyelitis, chronic fatigue syndrome, or chronic pain reported a hypomethylation trend [45,49], in contrast with our results, indicating possible different epigenetic regulation in FM compared with other chronic pain states.
Additionally, in the DMCs test, of the seventeen differentially methylated cytosines identified, only four were found less methylated in FM patients, while all the others were shown to be higher methylated in FM patients compared with controls, indicating the importance of the single base resolution technology as a valid approach to detect different trends in methylation within the same region. The genes evidenced in the DMCs test are related to neuronal development, dopaminergic and pain pathways, inflammation, and sociability. In particular, five differentially methylated cytosines confirmed GCSAML as possibly involved in FM pathogenesis. The DRD3 and IL25 genes, in which only one DMC was found, recall the involvement of the dopaminergic pathway and inflammation, respectively, in FM patients. Dopamine was shown to significantly influence pain perception, with striatal dopamine release associated with pain inhibition and DRD3 Ser9Gly polymorphism related to thermal pain perception in CWP patients [34,50]. The inflammatory cytokines IL25 was previously found to be up-regulated in FM [30]. The three cytosines identified in TRPA1 bring the focus on the results of Bell and coworkers, in which the promoter region methylation was inversely associated with both heat pain and pressure pain thresholds [26]. Two DMCs were evidenced in the proto-oncogene MCF2, which modulates the activity of small GTPases, and it is involved in dendrite elongation and neurite outgrowth [51]. Interestingly, MCF2 genetic and epigenetic variants were associated with FM and many other pathological and psychiatric diseases [52]. Five significantly differentially methylated (>1%) cytosines were instead detected in the oxytocin gene. Oxytocin is relevant for the perturbations in the HPA axis observed in FM patients, because it was shown to induce adreno-corticotropin-hormone release at the anterior lobe of the pituitary [8]. In healthy subjects, oxytocin decreases cortisol release and anxiety in response to social stress [53]; its anti-nociceptive, analgesic, anxiolytic, and sedative effects are well known [54,55].

Concomitant Risk Factors on FM Risk
As shown by the logistic regression model, testing the simultaneous incidence of DNA methylation changes, depression, and the clustered sociodemographic data on the risk to develop FM, GRM2 DNA methylation and depression were confirmed to increase FM risk. Interestingly, our results put the attention on the role of GRM2 gene methylation, which was also evidenced in our previous pilot study [16], although with not consistent results. GRM2, encoding the Glutamate Metabotropic Receptor 2, affects glutamate release, the major excitatory neurotransmitter in the central nervous system CNS, and thus might be involved in both central sensitization and immune/inflammatory pathological mechanisms [56,57]. In addition, even though the groups comparison did not yield any significant result after controlling for multiple testing, the network analysis allowed a more fine-grained interpretation of the results. In fact, we showed how GRM2 methylation represented a central node only in the FM sample, suggesting its relevance in the pathogenesis of the disorder. Moreover, this methylation site was not equally important in their healthy sisters that had a completely different structure. Nevertheless, the small sample size in our study and the correlational nature of the network analysis limits the generalizability of the results and does not allow inferring causality. Further modelling and investigations of networks' structural difference between siblings are needed.
Our results might support the hypothesis of an altered immune system response in FM. We propose that this altered pathway, in which GCSAML might have a role, could be cause or consequence of the complex FM phenotype. As shown in Figure 2, the HPA axis is the primary stress response system, and its activation results in downstream production of cortisol and a dampening of the immune response [57]. FM syndrome was found to be associated with hypocortisolism [58], and low cortisol levels may be associated to immune system hyper-reactivity with subsequent activation of inflammatory markers. Peripheral inflammatory mediators have been shown to directly induce transcriptional modulation in the brain [59]. Completing this loop, the CNS, in particular the brain stem catecholaminergic centers, may in turn regulate the HPA axis [60]. In addition, the HPA axis has been implicated in the pathophysiology of depression [61], in turn associated with peripheral inflammatory markers. Unfortunately, it is not possible to establish any causal relationship among the evidenced pathways, and future longitudinal designs are encouraged to clarify the contribution of the factors involved in FM.

Limitations and Future Research Directions
A major strength of this study is the inclusion of biological siblings unaffected by FM as controls, but certain limitations should be highlighted. First, it is still inconclusive if the identified DNA methylation differences are mechanisms of the disease or result from a response towards environmental stimuli. Early environmental stressors can cause CpGs hypermethylation, altering the HPA axis responses to stress: considering epigenetic factors with no correlation with personal life experiences can be deeply misleading. The second limitation is related to the analyzed population: it included 84 participants and may thus be too small to detect differential DNA methylation and strong associations with the participants characteristics. Further investigation into differential methylation between FM patients and healthy controls remains necessary. In addition, only women were included, and thus the results may not be generalized to male patients with FM. Third, potential transcriptional changes related to the altered methylation were not investigated because they required higher starting material, and this should be considered in future studies. Moreover, DNA methylation in peripheral blood is not necessarily directly reflective of central pain mechanisms, but it could serve as a peripheral epigenetic biomarker, as similar levels of DNA methylation were observed in blood and brain tissues at multiple pain regions in previous studies. Replication studies using specific brain and dorsal root ganglia tissues should help to further clarify the role of DNA methylation in FM.
inflammatory mediators have been shown to directly induce transcriptional modulation in the brain [59]. Completing this loop, the CNS, in particular the brain stem catecholaminergic centers, may in turn regulate the HPA axis [60]. In addition, the HPA axis has been implicated in the pathophysiology of depression [61], in turn associated with peripheral inflammatory markers. Unfortunately, it is not possible to establish any causal relationship among the evidenced pathways, and future longitudinal designs are encouraged to clarify the contribution of the factors involved in FM.

Figure 2.
Proposal of an etiopathogenesis model for FM: the HPA axis is the primary stress response system, and it is regulated by the CNS. In addition, depression and poor sleep/fatigue development, both characterized by elevated peripheral inflammatory markers, are influenced by the HPA axis.

Conclusions
The results of the present study recall the bidirectional communication between the brain and the immune system, and they are consistent with clinical data showing a complex involvement of depression in FM pathogenesis. FM seems to be the result of a complex interplay between stress system alterations that might trigger depression and pain pathway dysfunctions. In addition, we identified GCSAML and GRM2 as interesting targets that need to be considered in future research to unravel their role in FM and provide useful biomarkers to improve diagnosis and treatment for it.
Supplementary Materials: The following are available online at https://www.mdpi.com/article/ 10.3390/jcm10214992/s1, Table S1. Data collected for each female participant (C1-C42: healthy controls; P1-P42: Fibromyalgia patients): age, diagnosis code, Fibromyalgia Impact Questionnaire (FIQ), Widespread Pain Index (WPI), Symptom Severity scale (SSS), Fibromyalgia Severity (FS) score (the sum of the WPI and SSS), Visual Analog Scale (VAS) for main FM symptoms, Pittsburgh Sleep Quality Inventory (PSQI) and Beck Depression Inventory (BDI) scores. Average (AVER) and standard deviation (STDEV) are reported for the collected data., Table S2: Differentially methylated Cytosines (DMCs) test (Metilene output), including both significant and not significant p values: Metilene output reveals DNA methylation mean scores for patients (42 FM women) and controls (42 related healthy sisters) measured in FM symptoms related regions. Figure S1: Study design., Figure S2: The sequence analyzed related to the GCSAML gene: the region is on the first exon in a CpG region of the gene that has multiple transcripts.