Are Epigenetic Factors Implicated in Chronic Widespread Pain?

Background Chronic widespread musculoskeletal pain (CWP) is the cardinal symptom of fibromyalgia and affects about 12% of the general population. Familial aggregation of CWP has been repeatedly demonstrated with estimated heritabilities of around 50%, indicating a genetic susceptibility. The objective of the study was to explore genome-wide disease-differentially methylated positions (DMPs) for chronic widespread pain (CWP) in a sample of unrelated individuals and a subsample of discordant monozygotic (MZ) twins. Methodology/Principle Findings A total of N = 281 twin individuals from the TwinsUK registry, including N = 33 MZ twins discordant for self-reported CWP, were part of the discovery sample. The replication sample included 729 men and 756 women from a subsample of the KORA S4 survey–an independent population-based cohort from Southern Germany. Epigenome-wide analysis of DNA methylation was conducted using the Illumina Infinium HumanMethylation 450 DNA BeadChip in both the discovery and replication sample. Of our 40 main loci that were carried forward for replication, three CPGs reached significant p-values in the replication sample, including malate dehydrogenase 2 (MDH2; p-value 0.017), tetranectin (CLEC3B; p-value 0.039), and heat shock protein beta-6 (HSPB6; p-value 0.016). The associations between the collagen type I, alpha 2 chain (COL1A2) and monoamine oxidase B (MAOB) observed in the discovery sample–both of which have been previously reported to be biological candidates for pain–could not be replicated. Conclusion/Significance Our results may serve as a starting point to encourage further investigation in large and independent population-based cohorts of DNA methylation and other epigenetic changes as possible disease mechanisms in CWP. Ultimately, understanding the key mechanisms underlying CWP may lead to new treatments and inform clinical practice.


Introduction
Chronic widespread musculoskeletal pain (CWP) is the cardinal symptom of fibromyalgia (FM), emphasizing axial pain as a constant feature, as well as the presence of pain in the upper and lower quadrants of abdomen, and the right and left sides of the body [1]. Musculoskeletal pain affects about 12% of the general population, with the proportion of women being generally twice as high compared to men [2]. CWP is a diagnosis that challenges the social and cultural boundaries of current perception of illness and health in the Western world and the lack of a biomedical diagnosis denies the sufferer's legitimacy of the "sick role". CWP is not only a huge burden for sufferers, but also difficult to treat and costly for society overall. Despite recent research efforts, the CWP etiology and mechanisms remain unclear, with evidence from a range of epidemiologic studies suggesting an interplay of socio-demographic, psychosocial, affective, physiological, as well as genetic factors [3][4][5][6][7][8].
Familial aggregation of CWP has been repeatedly demonstrated with estimated heritabilities of around 50% [3,5]. Although several studies have been performed to search for biomarkers, they have been mostly underpowered, restricted by the heterogeneity of the phenotype, constrained by study design, and consequently remain unreplicated [9]. Nevertheless, several biochemical factors (e.g., androgens) [10] and candidate genes-in particular genes involved in neurotransmission and part of the dopaminergic and serotonergic pathways and/or important for the hypothalamic-pituitary-adrenal (HPA) axis-have been suggested as related to etiology of the condition [9]. Furthermore, the most recent and so far largest genome-wide association (GWAS) meta-analysis has found evidence for chaperonin containing TCP1 subunit 5 (CCT5) and family with sequence similarity 173 member B (FAM173B) as promising candidates in pain regulation of CWP [7].
Despite evidence of a genetic basis, environmental components seem to be equally important in the pathogenesis of CWP, suggesting a possible role of gene-environment interaction in the development of the condition [3,5]. Thus, exploration of the DNA epigenetic patterns that regulate gene expression profiles may offer a new model to understand the dynamic interactions between stochastic, environmental, genetic and epigenetic factors that influence an individual's musculoskeletal pain perception and expression. As such, specific epigenetic markers may provide quantifiable measures of lifetime environmental exposures and thus could offer a new biological framework for the cognitive-behavioral causes underlying CWP.
Recent research attempts implementing this methodology have provided promising novel insights into the etiological mechanisms underlying CWP/FM. In a genome-wide methylation study, Menzies and colleagues were able to observe significant differences in methylation patterns between 10 women with FM and 42 healthy controls, with the differentially methylated sites including brain derived neurotrophic factor (BDNF), N(alpha)-acetyltransferase 60, NatF catalytic subunit (NAT15), or protein kinase C alpha (PRKCA)-all of which show biological relevance to FM [11]. While these initial results are promising, they remain unreplicated, restricted to a small sample and therefore extremely underpowered. Furthermore, and to the best of our knowledge, no studies have looked at epigenetic markers for CWP. Given previous support for a role of methylation changes in FM-related pain, we conducted a detailed epigenome-wide analysis study (EWAS) of CWP-differentially methylated positions (DMPs) in peripheral blood DNA using the Illumina Infinium HumanMethylation 450 DNA BeadChip. Analyses were conducted on a sample of N = 281 unrelated individuals with available CWP data and in a subsample of 66 monozygotic (MZ) twins (33 twin pairs) discordant for CWP from the UK. Epigenetic studies of disease-discordant MZ twins, who are completely matched for genetics, age, sex, cohort effects, maternal influences and-to a certain extentcommon environment, are considerably more powerful in detecting disease-related epigenetic differences and therefore the ideal for similar studies [12]. The 20 most significant loci from both the complete sample and the twin pairs subsample were carried forward for replication in an iindependent sample of 729 men and 756 women.

Participants and Study Design
The discovery sample was drawn from the TwinsUK registry, which has been shown to be representative and comparable to the general population in terms of behavior, lifestyle factors and diseases [13]. For detailed information on the twin cohort, see [14,15]. Zygosity was confirmed by genotyping [16]. Collection of CWP phenotypes and genetic/epigenetic information (from blood samples) was carried out during clinical visits and via postal questionnaires. The twin subjects were unaware of the research interests of the present study. In the end, matching phenotype and epigenome-wide data was available for a total of N = 281 subjects. All twins provided written informed consent and the study was approved by St. Thomas' Hospital Research Ethics Committee.
The independent replication sample included 729 men and 756 women from a subsample of the KORA S4 survey, which was performed between 1999 and 2000. The KORA study includes a series of independent population-based surveys (S1-S4) recruiting participants from the region of Augsburg in Southern Germany and was initially set up as a cross-sectional survey to explore development and course of chronic diseases in a sample of adult individuals [17]. The KORA study was conducted according to the principles expressed in the Declaration of Helsinki. Written informed consent had been given by each participant. The study was reviewed and approved by the local ethical committee (Bayerische Landesärztekammer). Individuals from both the discovery and replication sample were of Caucasian background to limit biases resulting from ethnic disparities. Individuals from both the discovery and replication samples were of Northern European background to limit genetic heterogeneity.
Two different analyses were conducted in the discovery sample. First, methylation levels were compared in the full discovery sample, comprising 281 individuals. In addition, methylation status was compared in a subsample of 33 MZ twin pairs discordant for CWP. The study of epigenetics and disease using discordant MZ twins offers the opportunity to control for many potential confounders encountered in general population studies, such as differences in genetic background, early-life environmental exposure, age, gender, and cohort effects. It therefore represents a more advantaegous model and offers more power using a smaller sample to detect potential differences between disease-affected and non-affected individuals. The 20 most significant loci from both analyses (the discordant MZ subsample and the full discovery sample) were then carried forward for replication in the independent KORA sample (N = 1485).

CWP Phenotype
Different instruments were used in the replication and discovery cohort to assess and define the main phenotype CWP. In the discovery sample, the London Fibromyalgia Epidemiology Symptom Screening Questionnaire (LFESSQ) was used to screen for self-reported CWP [18]. This 6-item questionnaire was originally designed to screen for FM in the general population and in specific patient groups and includes 4 items relating to widespread pain, and 2 items relating to fatigue. For the assessment of CWP, the four items pertaining to the "pain subscale" and asking about pain left and right of body and above and below diaphragm lasting at least 7 days in the previous 3 months were considered. In order to be classified as having CWP, participants had to respond "yes" to all four pain items with either both a right-and left-side positive response or a positive response for pain at both sides. The utility of this phenotype assessment is supported by the contribution these twins have made to previous studies, such as the genome-wide association meta-analysis conducted by Peters and colleagues [7] or the twin studies by Burri and colleagues [5].
In KORA S4, pain was assessed as part of a self-report questionnaire based on the question "To what extent did pain hinder you in your daily tasks at home and at work during the past 4 weeks?," for which the participants could select from following answers: "no pain" (0), "not at all" (1), "slightly" (2), "moderately" (3), "quite" (4), "very" (5), "do not know." The pain variable was examined as a dichotomous variable, identifying controls (CWPq < 3) and cases (CWPq ! 3).

Illumina 450k DNA methylation analysis
Discovery Sample. Data and methylation analyses were conducted by common standard procedures in the discovery and the replication sample. In the discovery sample, DNA methylation analyses were performed in the R environment (http://cran.r-project.org) using raw intensity data files (IDAT files) as an input [19]. In all samples, the fraction of failed probes per sample did not exceed a quality contol cut-off (< 0.05). Multi-dimensional scaling and hierarchical clustering analyses were used to assess samples for outliers and no outliers were identified. Probes exhibiting detection p-values >0.01 were excluded from further analysis. Detection p-values predict whether the target probes signals are distinguishable from a set of background probe. For statistical analysis M-values (log2-ratios of methylated and unmethylated probes intensities) were used due to their presumed higher statistical validity [20]. Analyses of DNA methylation in 33 MZ twin pairs who screened discordant for CWP were carried out with the RnBeads package [21]. After removal of probes with detection p-values >0.01 in any of the samples, 483,017 probes were retained. Normalization was performed using the manufacturer-recommended scaling to internal controls as implemented in methylumi (method = "illumina"). Background correction was carried out with the Lumi option also as implemented in methylumi (bgcorr.method = "methylumi.lumi"). Paired-analysis of differential DNA methylation for the 33 MZ twin pairs was performed using paired Student's t-test. Age, array number, array position and processing batch were included as covariates in the analysis. The false discovery rate (FDR) after correction for multiple testing was set to 0.05. A priority-ranked list of differentially methylated CpGs was obtained, with rank based on a combination of statistical significance and absolute and relative effect size. The reason to select a rank approach in the discordant twin pairs analysis was to account also for effect sizes, which were the commonly very small between discordant twins.
DNA methylation patterns were further compared in a sample of unpaired twin individuals. One twin from each twin pair was chosen randomly, resulting in a sample of N = 281 unrelated individuals, out of which 200 screened negative and 81 screened positive for CWP. The data import, quality control and normalization steps of the minfi Bioconductor package were used [22,23]. Probes with detection p-values >0.01 in any of the samples were removed. Normalization and background correction were performed following the Illumina recommended procedures (bg.correct = TRUE, normalize = "controls") Q-Q plots of observed versus expected −log10 (p-values) in the unrelated CWP sample before and after normalization are shown in Fig 1. Differential methylation between the CWP positively and negatively screened individuals was estimated with the CpGassoc package [24]. Age, array number and array position were included as covariates in the analysis and processing batch was included as a fixed effect. The FDR (after Benjamini-Hochberg correction for multiple comparisons) was set at 0.05.
Replication Sample. In the replication sample, array-based DNA methylation data was also measured using the Infinium HumanMethylation450 BeadChip Array. The bisulfite conversion and genome-wide methylation assessment were performed as described in Zeilinger et al. [25]. Normalization of the methylation data was conducted following the pipeline of Lehne et al. [26] beginning with exclusion of 65 rs-probes (SNP probes) and background correction using the minfi R package. Probes were set to NA if the detection p-value was !0.01or if the number of beads was 3. Samples were excluded if the detection rate was 0.95. Quantile normalization, as given by Lehne et al. [26] then resulted in intensity values used to calculate the β-values. Following normalization, a site-by-site detection rate of 95% was applied to the KORA S4 cohort. To eliminate technical effects, the final methylation values for each CpG site were calculated as the residuals of a linear regression with the methylation value as response and the 20 control probe principal components (PCs) as covariates. The control probe PCs were calculated on all green and red intensities of all (non-negative, autosomal) control probes together, after background correction. The linear model for methylation (adjusted for the technical bias) and the binary CWP was estimated with the statistical software R, using age and sex as covariates.

Database for Annotation, Visualization and Integration Discovery (DAVID) analysis
Genes to which any of the top 20 differentially methylated CpG sites were annotated were entered in the Database for Annotation, Visualization and Integration Discovery (DAVID), a high-throughput and integrated data-mining environment. Functional annotation clustering was carried out and clusters with enrichment scores > 1.3 (a previously suggested threshold) are reported [27].

Results
The full final discovery sample of unrelated female twins (N = 281) consisted of 81 (28.8%) individuals having CWP (mean age ± SD was 62.7 ± 9.9 years). For the CWP-discordant MZ analyses, 33 full MZ twin pairs were available (mean age ± SD for the twin pairs was 61.8 ± 10.7). The final replication sample consisted of N = 1535. The overall mean age of the sample was 54.41 for men and 53.83 for women.

Differentially methylated positions in CWP discordant MZ twins
A SNP heatmap of the 65 genotyping probes on the Illumina Infinium Human Methylation 450 BeadChip which measures single nucleotide polymorphisms, was used to assess genotyperelated similarities between the samples. It showed clustering of each twin pair and no convincing additional clusters among the subjects (Fig 2). A list of the top ranked 20 DMPs between the CWP discordant twin pairs can be found in Table 1. Among the top 20 ranked DMPs we identified one mapped to monoamine oxidase B (MAOB; p = 0.08)-a gene previously implicated in neuropathic and postoperative pain [28,29]. Functional annotation clustering (FAC) of genes with hypermethylated CpG sites in CWP-affected MZ twins using the DAVID database further identified a cluster related to DNA binding, regulation of transcription, nucleus as well as developmental proteins (Enrichment score 1.81, genes: BANP, NKX6-2, T). FAC of genes with hypomethylated CpG sites in CWP-affected twins did not show any significantly enriched clusters (Table 2).

Differentially methylated sites in unrelated twin individuals
To identify DMPs in individuals with CWP compared to healthy controls in the unrelated sample, we fitted a fixed effects model regressing methylation levels at each probe on CWP  Table 3. These 20 DMPs included CpG islands (50%); north shores (15%); south shores (10%); as well as south shelves (5%) and sites that were not annotated into the previously noted categories (20%). Overall, the top 20 DMPs were localized around17 different genes ( Table 3). The highest statistical significance between the groups was found for cg11251499, which is annotated to the 1 st exon/5 0 -UTR of RE1-silencing transcription factor (REST), named also neuronrestrictive silencer factor (NRSF) (p-value 2.23 × 10 −6 )-a transcriptional repressor, which represses the expression of neuronal genes in non-neuronal cells [30][31][32]. REST has previously been implicated in pain perception [30]. Among the 20 most strongly associated DMRs, another previously reported biological candidate for pain-collagen type I alpha 2 chain (COL1A2)-could be detected. FAC of genes with hypomethylated or hypermethylated CpG sites in CWP positively screened individuals in the unrelated sample did not show significant enrichment.

Independent replication analyses of the 40 most strongly associated DMRs
The same Infinium HumanMethylation450 BeadChip Array was used for the replication of in an independent sample of n = 1485 individuals. Of these 40 DMRs, three reached EWAS significance in the indepenedent replication analyses. The three CpGs reaching significant p-values in the replication sample were malate dehydrogenase 2 (MDH2; p-value 0.017), tetranectin (CLEC3B; p-value 0.039), and heat shock protein beta-6 (HSPB6; p-value 0.016; Table 4). The significant associations between the previously reported biological candidates for pain identified in the discovery cohort-including REST, MAOB, COL1A2-could not be successfully replicated (Table 4).

Discussion
Several interesting candidates could be identified in our discovery analyses, including REST, MAOB, COL1A2. Although none of these candidates were successfully replicated in our independent replication sample, further investigation of these candidates is warranted given previous reports of their involvement and importance in pain regulation. Although small in number, study findings in animal models and humans, for example, have repeatedly underlined the analgesic properties of reversible MAO-A inhibitors in the treatment of nociceptive and psychogenic pain [33][34][35]. The majority of the evidence stems from pharmacological trials investigating the therapeutic effects of MAO inhibitors on patients reporting different pain disorders. Although they are not officially approved for this indication, antidepressantsincluding tricyclic antidepressants (TCAs), selective serotonin reuptake inhibitors (SSRIs), and MAOIs-have been used in the treatment of FM. A study by Pillidar and colleagues [34], for example, found moclobemide-a reversible inhibitor of monoamine oxidase type A-to be an effective option in the management of perceived pain, as well as of symptoms of somatization. Similarly, a review by Tort and colleagues [36] suggested a moderate effect of MAOIs on pain and a small effect on tender points. Supporting this conclusion, a meta-analysis conducted by Häuser and colleagues found medium effect sizes for MAOIs for pain reduction [37]. Further evidence for an involvement of MAO in pain comes from genetic studies, such as the experimental study conducted by Treister et al. [38]. Using a Transmission Disequilibrium Test (TDT) on 192 individuals, the authors found significant associations between cold pain tolerance and MAO-A polymorphisms, suggesting an involvement of MAO in pain sensitivity. Similarly, the significant associations between the A/G polymorphism of the MAOB gene and postoperative pain intensity found by Sery et al., also point towards a potential role of MAOB in the perception of pain intensity [28].
Results from our larger sample of unrelated individuals also revealed a high, albeit non-significant, association between COL1A2 and CWP. The COL1A2 gene provides instructions for the creation of type I collagen, which is known to strengthen and support many tissues in the body (including cartilage, bone, tendon, skin, etc.). A number of previous studies have shown a link between collagen and FM. The 27 female FM patients included in a muscle biopsy study by Gronemann et al. [39], for example, showed significantly lower amount of intramuscular collagen compared to healthy controls. The authors argued that the reduced collagen may lower the threshold for muscle micro-injury and consequently resulting in non-specific signs of muscle pathology. Further exploring the distribution and amount of collagen in skin from non-tender-point areas in N = 27 female FM patients, Ribel-Madsen et al. detected differences between the amino acid composition of skin proteins, i.e., fewer amount and less density, in FM patients compared with controls [40]. Yet another study conducted on 20 FM patients even found a marginal symptomatic improvement of chronic FM symptoms by taking collagen hydrolysat-a prescription-free food supplement [41].
Functional annotation clustering of our CWP-discordant MZ results further revealed enrichment in genes with hypermethylated CpGs in the CWP-affected twins for DNA binding, regulation of transcription, nucleus as well as developmental proteins (i.e., BANP, NKX6-2). These findings are in line with a previous study by Menzies and colleagues which identified enrichment in differentially methylated sites to clusters involved in neuron differentiation, development (nervous and skeletal/organ system) and chromatin compaction in women with FM [11].
Of our 40 DMRs that were carried forward for replication, three CpGs reached significant p-values in the replication sample and included MDH2, CLEC3B, and HSPB6. While, to the best of our knowledge, none of them have previously been related to pain perception or expression, a recent comparative analysis of gene expression in tendinopathy found evidence of CLEC3B (and COL1A2) upregulation in tendon cells in response to cyclic strain [42]. Widespread tenderness of muscles and the areas around tendon insertions (as well as increased prevalence of tendinitis) is a common clinical presentation in patients with FM and CWP and an involvement of CLEC3B in these musculoskeletal and connective tissue disorders is likely [43].

Limitations
Several study limitations have to be considered such as the relatively small sample size of the discovery sample, especially the subsample of discordant MZ twins. This might be one reason why statistical significance at FDR < 0.05 was not reached. Furthermore, comorbidities such as anxiety and depression were not taken into account, as this would have resulted in very low statistical power. Due to the relatively high mean age of our sample, extrapolation of the detected differences to other study populations is not necessarily possible. In the replication and discovery cohort, different instruments were used to assess and define the main phenotype CWP. This could be one reason why replication of the main findings fell short. However, we believe that it can equally be seen as a strength rather than a weakness that we were able to partially replicate our findings although phenotype assessment (and study population) were different. Furthermore, the different mean ages of the discovery and replication sample, as well as the fact that the cohorts were of different cultural background (UK vs Germany)-albeit homogenous in terms of ethnic background-could be other reasons why the main findings did not replicate. Finally, we used blood samples for the detection of DNA methylation patterns. It is well known that epigenetic patterns are tissue-specific, therefore the study of central and peripheral nervous system tissues would be of great interest, but are very hard to obtain in humans [44].

Conclusion
In this large EWAS study of CWP we found evidence for the involvement of epigenetic factors. These preliminary results serve as a starting point to encourage further replication in large and independent population-based cohorts and subsequent biological investigation to elucidate possible disease mechanisms. Ultimately, understanding the key mechanisms underlying CWP may lead to new treatments and inform clinical practice and developments in psychiatric nosology.