An omics investigation into chronic widespread musculoskeletal pain reveals epiandrosterone sulfate as a potential biomarker.

Chronic widespread musculoskeletal pain (CWP) is common, having a population prevalence of 10%. This study aimed to define the biological basis of the CWP/body mass association by using a systems biology approach. Adult female twins (n = 2444) from the TwinsUK registry who had extensive clinical, anthropometric, and "omic" data were included. Nontargeted metabolomics screening including 324 metabolites was carried out for CWP and body composition using dual-energy X-ray absorptiometry. The biological basis of these associations was explored through a genome-wide association study and replicated in an independent population sample (Cooperative Health Research in the Region of Augsburg [KORA] study, n = 2483). A causal role for the genetic variants identified was sought in CWP using a Mendelian randomisation study design. Fat mass/height2 was the body composition variable most strongly associated with CWP (TwinsUK: P = 2.4 × 10(-15) and KORA: P = 1.59 × 10(-10)). Of 324 metabolites examined, epiandrosterone sulfate (EAS) was highly associated with both CWP (P = 1.05 × 10(-09) in TwinsUK and P = 3.70 × 10(-06) in KORA) and fat mass/height2. Genome-wide association study of EAS identified imputed single nucleotide polymorphism rs1581492 at 7q22.1 to be strikingly associated with EAS levels (P ≤ 2.49 × 10(-78)), and this result was replicated in KORA (P = 2.12 × 10(-9)). Mendelian randomization by rs1581492 genotype showed that EAS is unlikely to be causally related to CWP. Using an agnostic omics approach to focus on the association of CWP with body mass index, we have confirmed a steroid hormone association and identified a genetic variant upstream of the CYP genes, which likely controls this response. This study suggests that steroid hormone abnormalities result from pain rather than causing it, and EAS may provide a biomarker that identifies subgroups at risk of CWP.

1. Introduction

Chronic widespread musculoskeletal pain (CWP)—a key feature of fibromyalgia—is common in the general population with a prevalence of 5% to 15%.[21,28] The condition is associated with frequent physician consultations and high public health cost.[31] Of the many aetiological factors that have been proposed in CWP, the association with elevated body mass index (BMI) is one of the strongest and most consistently reported and has been demonstrated both in cross-sectional[14,19] and longitudinal studies.[18] The observation that raised BMI predicts the onset of pain suggests that BMI may be casually related to CWP and not simply a reflection of weight gain resulting from, for example, pain-induced immobility. We have shown previously that there is a common genetic basis for multisite musculoskeletal pain[33] and have recently contributed to an international genome-wide association study (GWAS) meta-analysis for CWP.[22] The power of “agnostic” metabolomics as an intermediate phenotype in the association analysis of complex traits has been demonstrated.[29]

In this study, a systems biology approach was used to dissect the biological basis of the relationship between BMI and CWP. We used omics data (genomics and metabolomics) to perform agnostic testing of 2 large population samples of Northern European origin. First, we defined the relationship between BMI and CWP using whole-body dual-energy X-ray absorptiometry to define the component of BMI (fat vs muscle) mediating the relationship. Then, we examined a panel of 364 metabolites to determine associations with both CWP and body composition. Identified metabolites were further explored through GWAS in a strategy that has proved useful in other common complex traits[29] including age-related traits[16] and diabetes.[15] The combined omic approach represents a novel and powerful attempt to apply new technological methods[10] to the elucidation of biological pathways underlying CWP.

2. Methods

The study samples included the TwinsUK as the discovery set with replication performed using data available in a sample from the Cooperative Health Research in the Region of Augsburg (KORA) study.[8]

2.1. Study samples—phenotypes

The TwinsUK discovery sample available for this study included 1269 monozygotic and 1175 dizygotic twins.[17,27] The register of healthy adult female twins has been developed over 22 years, and participants are representative of the general UK population.[2] Twins are sent regular questionnaires and may be invited to a clinical visit and sample collection. Where possible, twins were not made aware of the precise hypothesis being tested before inclusion in a particular study.

Twins included in this study had been sent the London Fibromyalgia Epidemiology Symptom Screening Questionnaire (LFESSQ) for self-completion without reference to the co-twin.[32] Twins with pain on both left and right sides of the body, above and below the diaphragm, and with a duration of 7 days or more within the preceding 3 months were considered as cases. At the clinical visit, height was measured in metres and weight in kilograms, and BMI was calculated (kg/m2). Participants also underwent a whole-body dual-energy X-ray absorptiometry (Hologic Discovery W; Hologic, Bedford, MA) following the manufacturer's recommendations.[3] This method quantified the contributions of lean body mass and fat body mass.

Participants from KORA included 1231 men and 1252 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. The KORA study is a cross-sectional survey of the development and course of chronic diseases in a randomly collected sample of adult individuals.[8] In KORA, 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 twice, first, as a dichotomous variable with controls (CWPq < 3) vs cases (CWPq ≥ 3) and then as a continuous trait. For the assessment of lean mass and fat mass, bioelectrical impedance analysis (BIA) measurements of resistance (R), reactance (Xc), and the phase angle (α) were taken using a bioelectrical impedance analyser (BIA 2000-S; Data Input GmbH, Frankfurt, Germany). Based on these measurements, lean and fat mass were calculated by Segal's equations.[25,26]

2.2. Metabolomics

Nontargeted ultrahigh-performance liquid chromatography and mass spectrometry was performed on fasting plasma samples of TwinsUK participants (n = 5003) and fasting serum samples of KORA S4 participants (n = 1614) using the Metabolon's platform.[24,29] The metabolomic data set contained amino acids, peptides, acylcarnitines, sphingomyelins, glycerophospholipids, lipids, carbohydrates, nucleotides, vitamins, steroids, and xenobiotics. Raw data were median normalized for daily fluctuations of the method and then inverse normalized. In the discovery set, missing values were imputed using the minimum measures of each metabolomic feature for the particular run days. To avoid spurious false-positive associations because of small sample size, metabolic traits having >20% missing values were excluded.

2.3. Genomics

TwinsUK subjects had been genotyped for association markers using a combination of Illumina arrays (Human Hap300 and the Human Hap610Q) as previously reported.[12] For this analysis, single nucleotide polymorphisms (SNPs) were excluded if call rate <97% (SNPs with minor allele frequency, MAF ≥5%) or <99% (for 1% ≤ MAF <5%), Hardy–Weinberg p values <10−6, and MAF <1%. Subjects were removed if genotyping failed in >2% SNP. The overall genotyping efficiency was 98.7%. Imputation of genotypes was carried out using the software IMPUTE version 2[13] using HapMap II as the reference panel. Population substructure was examined using the STRUCTURE program[23] and correspondingly controlled for spurious associations. In KORA, Affymetrix Axiom chip had been used for genotyping. HapMap build 37 served as population reference, and the criteria of call rate >98% and P (Hardy–Weinberg) >5 × 10−6 were applied as filters for SNP quality. Genotyped SNPs were imputed with IMPUTE v2.3.0 using the 1000G set as reference panel.

2.4. Analytical approach

The analysis was conducted in several steps. First, potential risk factors for CWP including age and anthropometric measurements (weight, height, BMI, and body composition variables) were examined in univariate analysis (Student t test). Second, we sought metabolites significantly associated with CWP, preserving power by including only those metabolites available in ≥2000 subjects in TwinsUK. To this aim, all 324 metabolites available were tested separately in a series of t tests comparing each metabolite among the affected vs nonaffected individuals. Metabolites showing metabolome-wide significant association with CWP (false discovery rate <5%[5]) were tested for association with BMI and body composition variables and then underwent GWAS using GenABEL software,[4] adjusting for familial relationship.

To determine the relative contributions of factors influencing CWP, we conducted binary logistic regression, taking into account twin relatedness. Variables that were introduced stepwise into the regression model included age, BMI/body composition, metabolites, the most highly associated genotyped SNP, and co-twin status. This analysis included 4 metabolites that had been identified at the previous stage as statistically significant and independently associated with CWP. We also note that GWAS revealed only a few highly significant “candidate” SNPs, which were included in the binary logistic regression as potential covariates.

To assess the causal influence of the identified metabolites on CWP, we performed Mendelian randomisation analysis using the instrumental variable method.[20] Multivariate regression analysis of the selected metabolite(s) with simultaneous adjustment for the top SNP and other significant covariates was used. In KORA, where the sample comprised unrelated individuals and pain phenotype was scored as a semiquantitative variable (CWPq), we used multivariable linear regression analyses to determine the association of CWP with all potential covariates, including age, body composition variables, and circulating levels of metabolites. In addition, the CWP was further classified as a dichotomous variable—controls (CWPq < 3) vs cases (CWPq ≥ 3).

3. Results

3.1. Body mass phenotype

The discovery sample comprised 2444 TwinsUK participants and the replication sample comprised 2483 participants from KORA (of which 1614 had metabolomic data). Table 1 shows a comparison of the age and body composition risk factors by case–control status for the 2 groups. Controls were found to be significantly younger (Twins UK: 53.5 vs 58.3, P = 7.8 × 10−17; KORA: 56.9 vs 59.6 years, P = 2.3 × 10−10), taller, and leaner than cases of CWP. Chronic widespread musculoskeletal pain cases demonstrated greater adiposity than controls in all body composition variables of which relative fat mass (fat body mass/H2—referred to as fat mass index or FMI hereafter) was most significantly associated with CWP in both samples (TwinsUK: P = 2.4 × 10−15, KORA: P = 1.7 × 10−08) and was used as the main variable in subsequent analyses (Table 1).

Table 1

Descriptive statistics for the discovery and replication samples from TwinsUK and KORA, respectively.

3.2. Metabolomics

Of the 324 metabolites having complete data in TwinsUK, 6 showed significant association with CWP (after adjusting for age, false discovery rate <0.05; Table 2). Epiandrosterone sulfate (EAS) was the most strongly associated before (P = 6.81 × 10−29) and after adjustment for age (P = 1.05 × 10−09), demonstrating a strong inverse association, that is, lower levels of EAS associated with higher risk of CWP. In univariate analyses, 5 of the 6 metabolites found associated with CWP were also significantly correlated with FMI, the exception being unidentified metabolite X-1440 (Table 2). The significant associations of all 6 metabolites with CWP survived adjustment for other covariates (age, FMI, twin relatedness) in multiple logistic regression (Table 2). The corresponding OR ranged between 0.676 (P = 1.05 × 10−09) for EAS and 0.838 (P = 4.84 × 10−03) for nonadecanoate (19:0).

Table 2

Association of metabolite levels with CWP and with relative fat mass.

In KORA, similar associations were observed (Table 2). Testing pain as a continuous phenotype, with simultaneous adjustment for age, sex, and FMI, confirmed the association of pain with 4 of the 5 metabolites, including EAS (P = 3.70 × 10−06) and DHEAS (P = 1.86 × 10−06). Of note, on subgroup analysis by gender in KORA, the associations were consistently significant and in the same direction in females and males (data not shown).

3.3. Genomics of metabolites

The discovery GWA study of EAS using 2.5 million genotyped and imputed SNPs in TwinsUK revealed a highly significant association at a single peak on chromosome 7q22.1 between 98.85 and 99.05 Mbp (P ≤ 2.49 × 10−78; Fig. 1). The lead SNP rs1581492 was imputed and was in strong linkage disequilibrium with several other genotyped and imputed highly associated SNPs (for the 6 genotyped SNPs D′ = 0.89-1.00 (Fig. 2)). The proportion of variance in EAS explained by the variants varied between 7.6% and 8.1%. In silico replication in KORA provided confirmation of the signal with the overlapping SNPs showing an almost identical pattern, with the lead SNP rs1581492 having P = 2.12 × 10−09 (Fig. 1B).

Figure 1

Manhattan plot of GWAS of metabolite M33973 in TwinsUK and KORA. The upper figure shows the GWAS results for EAS in TwinsUK*, with KORA results below. *Including only unaffected individuals in the GWAS gives almost identical results, although the corresponding P values were attenuated.

Figure 2

Regional plot of association results and recombination rates for the chromosome 7q22.1. −log10 P values (y axis) of the SNPs are shown according to their chromosomal positions (x axis) with lead SNP shown as a purple diamond. The colour intensity of each symbol depicting an SNP reflects the extent of LD with the rs1581492, coloured red (r2 > 0.8) through to blue (r2 < 0.2). Genetic recombination rates (cM/Mb), estimated using HapMap CEU samples, are shown with a light blue line. Physical positions are based on build 36 (NCBI) of the human genome. Also shown are the relative positions of genes mapping to the region of association. Genes have been redrawn to show the relative positions and, therefore, the maps are not to physical scale.

For nonadecanoate (19:0), we found no genome-wide significant association signals. The results of GWAS of 3-(4-hydroxyphenyl)lactate revealed a single genome-wide significant peak (P ≥ 10−9), which mapped to chromosome 17 between 77.67 and 77.77 Mbp. Of 9 SNPs genotyped in this region, 8 were in almost perfect LD (D′ ∼ 1, r2 > 0.99). Although the nearest gene is CCDC57, the SLC16A3 locus also lies in that region and is more likely to harbour the causal variant, as this gene is known to encode a protein that transports lactate and derivatives thereof.

Finally, all covariates significantly associated with CWP were included in a single multiple logistic regression model simultaneously to test their relative contribution to CWP risk (Table 3). In particular, we tested the independent association of genotyped SNP, rs10235235, which was virtually in perfect LD with our top but imputed SNP rs1581492. The analysis revealed that besides age and FMI, EAS, SNP rs10235235, and co-twin affection status (residual genetic effect) made major contributions to overall model significance (χ2(10) = 340.8, P << 0.001) and the prediction of CWP. When all 6 metabolites were tested simultaneously, only 3 remained significantly and independently associated with CWP. Of these, EAS was highly associated with CWP (P = 5.9 × 10−10) in a combined model including age, co-twin status, FMI, and SNPs (Table 3).

Table 3

Relative contribution of metabolomic and genomic predictors to risk of CWP in TwinsUK.

3.4. Mendelian randomisation

Taken together, the results suggested that EAS levels influence the development of CWP. To test this hypothesis, we performed a Mendelian randomization study of the lead genotyped SNP. If circulating levels of EAS genetically determined by rs10235235 are causally related to CWP, individuals carrying the rs10235235 C allele would be expected to have reduced EAS levels and a higher prevalence of CWP. In fact, our results showed that CWP prevalence among C allele carriers was lower compared with noncarriers (0.17 vs 0.21, P = 0.046). Conducting the 2 stage instrumental variable analysis using genotype-predicted EAS levels showed that predicted levels do not significantly influence the risk of CWP (β = 0.458 ± 0.239, P = 0.055), whereas EAS residuals, after adjustment for genotype effect and other covariates, were significantly associated with CWP (β = −0.512 ± 0.066, P = 2.01 × 10−14). These finding are consistent with EAS levels falling as a consequence of pain, rather than predisposing to it, determined by the genotype at 7q22.1. However, the presence of an unobserved confounding variable influencing the outcome of this analysis cannot be excluded.

4. Discussion

Fibromyalgia is a highly prevalent health problem in the EU and USA comprising CWP, fatigue, and sleep disturbance.[6] Chronic widespread musculoskeletal pain is recognised to coexist with other common pain states, and they are thought to share a genetic underlying predisposition.[30] We used novel omics technologies to dissect the biological mechanisms underlying the CWP–BMI relationship. Our analysis, replicated in an independent sample, has shown that adiposity or FMI is the major mediating body composition factor. The associations we observed were consistently significant in the 2 data sets, despite differences in their phenotyping of both CWP and body composition. Metabolomic data analysis revealed joint independent associations between steroid hormones, FMI, and CWP in both TwinsUK and KORA cohorts with an inverse association of CWP risk with EAS—that is, the risk of CWP increasing with falling levels of EAS.

We are not the first to identify steroid hormone metabolism abnormalities in chronic pain.[1,7] However, our approach had the advantage of being an agnostic search of metabolites in a large population sample and having an independent population cohort for replication. Although CWP was associated with cortisone levels on metabolomic screen (P = 0.00069), the association was only marginally statistically significant when adjustment was made for FMI (P = 0.052, data not shown). Taken together, our data provide compelling evidence of a strong association between androgen hormone metabolism and CWP, suggesting that chronic pain leads to a reduction in hormone levels.

The use of metabolites as an intermediate phenotype in GWAS is providing a tractable approach to understanding better the genetic variants, and hence the pathways, involved in common complex traits.[29] GWAS of EAS levels revealed multiple highly associated SNPs on chromosome 7q22.1 in TwinsUK, a finding that was replicated in KORA, providing robust evidence of a true association. The lead genotyped SNP rs10235235 explained 8% of the total variance of EAS, and it was our expectation that this represented a novel genomic locus predisposing to CWP. Note, however, that the inclusion of the top imputed SNP rs1581492, instead rs10235235, gives virtually the same result. Chromosome 7q22.1 is a gene-rich region and includes both the zinc finger gene ZNF789 and cytochrome p450 gene CYP3A5, which lies 0.3 Mbp upstream of the SNP (Fig. 2). The latter is known to be involved in intracellular drug metabolism[11] and synthesis and breakdown of a variety of lipids including cholesterol and steroid hormones and so very likely has an influence on the androgen steroid metabolism pathway containing both dehydroepiandrosterone and its breakdown product, EAS. Mendelian randomization analysis, however, did not confirm that EAS lies in the causal pathway for CWP. Explanations for this include lack of power to detect a real effect: power for this analysis was estimated at 56%. Alternatively, there may be pleiotropic effects of SNP rs10235235 on CWP and EAS or else reverse causation: that CWP leads to reduced EAS levels. Finally, it is possible that other factors (genetic loci or epigenetic influences) could be responsible for this consistently observed correlation. Either way, falling EAS levels may represent a sensitive marker of chronic pain manifestation and has the potential greatly to assist in the clinical management of CWP.

There are a number of limitations of the study. Questions used to define the pain phenotype in the TwinsUK and KORA collections differed, as did the precise definition of CWP that was applied. This is a field of study where standardisation of phenotype is greatly needed; however, both samples have contributed to the successful GWAS meta-analysis of CWP,[22] which was limited by similarly varied diagnostic criteria from among the many contributors. The advantage of using population cohorts is the availability of large sample sizes—which are essential in omic studies. With less well-defined conditions, the trade-off with phenotype consistency would serve to bias findings towards the null—and lessen the chance of positive findings or successful replication. TwinsUK has CWP prevalence high enough to provide informative case/control numbers, and the questions used were taken from a validated questionnaire (LFESSQ),[32] whereas the CWP data for KORA consistently replicated the findings regardless of how the phenotype was defined (semiquantitative or dichotomous). Our findings are pertinent to chronic pain in the community as both TwinsUK and KORA are population samples. TwinsUK has been shown to be similar for common traits and outcomes to age-matched singleton women.[2] These results might also extend to other chronic pain syndromes such as irritable bowel syndrome and chronic pelvic pain, and this needs to be investigated. Other limitations include the predominance of females in the TwinsUK sample, something that might be regarded as an advantage here as CWP is more prevalent in women. It is noteworthy that in KORA, a similar relationship between body composition and metabolites in men and women was observed—despite the use of a different method to determine body fat. These results suggest that our findings are robust and pertinent to both sexes.

To date, genetic studies of pain phenotypes have explored an array of candidate genes taken from putative and diverse biochemical pathways.[9] By narrowing our focus to the association between CWP and BMI/body composition, we have shed light on possible neuroendocrine and central hormonal mechanisms that are shared between the 2 traits. The findings show that CWP influences EAS in a highly genotype-dependent manner, the specific nature of which remains to be established. It is possible that this effect is mediated by the cytochrome P450 enzyme 3A5, the gene, which is mapped to the genomic region highly associated with EAS variation. Further exploration of the chromosomal region is under way, and the use of EAS as a clinical biomarker in CWP and other chronic pain states will need to be assessed in independent samples.

Conflicts of interest statement

The authors have no conflicts of interest to declare.