Common hidden genetic risks linking fibromyalgia and mental health traits revealed by advanced genetic analysis

All GWASs included in this study had been approved by a relevant review board. No additional ethics approval was required as this study was conducted based on summary-level statistics.

The study comprised several key steps: [1] a phenotype-specific GWAS meta-analysis was conducted on the phenotypes of fibromyalgia, insomnia, depression, and anxiety, respectively, to increase statistical power; [2] Genomic SEM was employed to identify the common genetic factor underlying these conditions; [3] downstream analyses were performed to pinpoint lead single nucleotide polymorphisms (SNPs), as well as the most likely implicated genes, pathways, and tissues associated with this common factor; [4] the Latent Causal Variable (LCV) method was applied to identify modifiable risk factors for the common factor and to determine which diseases might be influenced by it; [5] brain-wide Mendelian randomization (MR) was utilized to identify brain regions linked to the common factor; [6] proteome-wide MR was conducted to identify potential blood biomarkers, and [7] multi-layer brain molecular quantitative trait locus (QTL) analyses was conducted to investigate the underlying mechanisms of mvFibroPsych. The flowchart of the study design is presented in Fig 1.

We conducted phenotype-specific GWAS meta-analyses for fibromyalgia, insomnia, depression, and anxiety, respectively. GWAS data for fibromyalgia were sourced from the BioMe Biobank [11], while data for insomnia and depression were obtained from the UK Biobank [15,16], and anxiety data from the iPSYCH consortium [17]. Additionally, corresponding GWAS datasets for these phenotypes were retrieved from FinnGen R11 to further strengthen the analysis [18]. Comprehensive quality control procedures were implemented across all datasets, including pre-imputation filtering, imputation, and post-imputation assessments, with adjustments for age, sex, and principal components to control for population stratification. All four phenotypes—fibromyalgia, insomnia, depression, and anxiety—were modeled as binary (case–control) traits in the meta-analysis. Fibromyalgia cases were defined based on the ICD-10 code M79.7. Insomnia, depression, and anxiety were defined according to cohort-specific criteria, including self-reported symptoms, clinical interviews, or diagnostic codes, as described in the original GWAS publications from which the summary statistics were obtained. Relevant references and sample size details are provided in S1 Table. We computed the effective sample size (Neff) for each cohort using the standard case–control formula [13]. For meta-analyses of multiple cohorts, we computed Neff for each component dataset separately and summed them to generate the total effective sample size per phenotype. The meta-analysis was performed using the fixed-effect inverse-variance-weighted approach in METAL [19], ensuring robust and reliable combined estimates across studies. All subsequent analyses were conducted based on the summary statistics from this meta-analysis.

The GWAS of the latent common factor for fibromyalgia, insomnia, depression, and anxiety (mvFibroPsych) was conducted by jointly modeling the cross-trait liability of these conditions using a common factor model within the framework of Genomic SEM [12]. We selected the common factor model based on prior evidence of substantial genetic correlations among the included traits, and in line with prior GenomicSEM studies that have modeled shared latent genetic liabilities [12,20–22]. The LDSC was undertaken to investigate the genetic correlations across the four traits. The Genomic SEM integrates genetic correlations and SNP heritability derived from GWAS summary statistics of individual traits, even when the sample overlap is variable or unknown. This method enables the modeling of multivariate genetic associations across phenotypes and can identify variants that influence shared genetic liability across multiple traits [12]. Prior to the multivariable GWAS, allele alignment across fibromyalgia, insomnia, depression, and anxiety was performed using the HAPMAP3 reference panel. Quality control was applied by selecting SNPs with a minor allele frequency greater than 0.01 and an INFO score above 0.9. The multivariable LDSC was then estimated across fibromyalgia, insomnia, depression, and anxiety using the 1000G European reference panel.

The multivariate GWAS was implemented using the commonfactorGWAS() function in GenomicSEM. The latent factor model was estimated using diagonally weighted least squares (DWLS). The regression is performed at the summary-statistics level, and no additional covariates are included at this stage, as covariate adjustment (e.g., age, sex, PCs) had already been applied in the original GWAS of each trait. The analysis was restricted to individuals of European ancestry to ensure LD homogeneity. The effective sample size for the mvFibroPsych GWAS was calculated following the GenomicSEM recommendations. Finally, we integrated the LDSC output with the fibromyalgia, insomnia, depression, and anxiety summary statistics to perform the multivariable common factor GWAS. The GWAS output included SNP-level estimates of beta, standard error, Z-score, and two-sided P-value for the association with the latent factor. In addition, we computed Q_SNP statistics and associated Q_p-values to test for heterogeneity in SNP effects across the contributing traits. The Q_SNP statistic evaluates whether the effect of a SNP is fully mediated by the common factor or whether trait-specific residual effects remain. SNPs with significant Q_SNP values indicate potential heterogeneity and were reported separately.

To evaluate model stability and trait-specific contribution, we performed sensitivity analyses using multiple combinations of traits. These included (i) a three-trait model excluding fibromyalgia, (ii) three-trait models pairing fibromyalgia with two psychiatric traits, and (iii) user-specified pairwise models with factor loading of fibromyalgia was fixed to 1. Model fit and factor loading significance were used to assess robustness and overlap structure. All models were run using the GenomicSEM v0.0.5 R package with default settings. Model fit was evaluated using Comparative Fit Index (CFI) and Standardized Root Mean Square Residual (SRMR). All multivariate GWAS results were independently replicated by multiple members of the research team using separate computational pipelines and software installations.

The post-GWAS analysis of the mvFibroPsych was conducted using the Functional Mapping and Annotation (FUMA) platform (https://fuma.ctglab.nl/home↗) to perform SNP annotation, gene prioritization, and enrichment analyses [23]. SNPs were annotated based on their functional consequences using reference databases, including Annotate Variation (ANNOVAR) and The Combined Annotation-Dependent Depletion (CADD) scores, to identify likely causal variants [24,25]. Gene prioritization was carried out by mapping significant SNPs to genes through positional mapping, expression quantitative trait loci (eQTL) mapping, and chromatin interaction data. Additionally, gene-based association analysis was performed using the Multi-marker Analysis of Genomic Annotation (MAGMA) to identify genes associated with the mvFibroPsych. Tissue enrichment and pathway enrichment analyses were then undertaken to determine the biological pathways and tissues most relevant to the identified genes, utilizing databases such as GTEx for tissue-specific expression and MsigDB for pathway annotations. Standard parameters and thresholds were applied throughout these analyses to ensure the robustness of the results.

The LCV method was applied to infer causal relationships between the common factor and other traits phenome-wide. Given its ability to leverage genome-wide genetic correlation patterns without requiring strong instruments, LCV was specifically selected for identifying potentially causal directions between mvFibroPsych and complex traits, especially where effect size estimation was not the primary aim [26]. The LCV model introduces a latent variable (L) that is assumed to exert causal effects on both traits through their genetic correlation. This model evaluates whether one trait potentially causes the other (vertical pleiotropy) by comparing the strength of association between L and each trait. An absolute genetic causality proportion (GCP) value exceeding 0.6 and a p-value below the FDR-corrected (Benjamini–Hochberg false discovery rate) threshold is considered evidence of a causal link between the two traits. A positive GCP suggests that trait A is likely to influence trait B, while a negative GCP implies the reverse [26]. The phenome-wide causal relationships between the common genetic factor identified in our analysis and various traits were examined using GWAS data from the UK Biobank, specifically the second wave of data released by the Neale Lab (www.nealelab.is/uk-biobank/↗), which primarily includes individuals of European ancestry. The covariates adjusted in the UK Biobank GWAS data have been extensively described elsewhere [27].

The two-sample Mendelian randomization framework was employed to quantify the direction and magnitude of putative causal effects between Brain imaging-derived phenotypes (IDPs) and mvFibroPsych. GWAS summary statistics for IDPs were obtained from a cohort of 33,224 individuals of European ancestry, as released in the 2020 UK Biobank dataset [28,29]. IDPs were obtained from structural MRI (sMRI) and diffusion MRI (dMRI) modalities, with sMRI capturing brain anatomical variations and dMRI assessing structural connectivity between brain regions. We filtered the initial 3,935 IDPs using a stepwise approach to ensure reliable and reproducible results, following the methodology by Yang et al [30]. First, we addressed redundancy by removing 1,201 IDPs measured in the same brain regions with different tools, retaining only those assessed by commonly used methods such as FreeSurfer, FIRST, tract-based spatial statistics, and probabilistic tractography. Next, we excluded 358 IDPs derived from small brain areas with low contrast in MR images. After applying these filters, we selected 227 sMRI and 360 dMRI IDPs for further analysis.

For assessing the impact of IDPs on mvFibroPsych, we applied a more lenient set of parameters for instrumental variable (IV) selection (p-value threshold of 5 × 10^-6, an r² threshold of 0.001, and a 1 Mb window size). This approach was necessitated by the relatively smaller sample size of the IDPs, which limited the number of SNPs reaching the conventional genome-wide significance threshold of 5 × 10^-8. Conversely, for evaluating the influence of mvFibroPsych on IDPs, we utilized more stringent IV selection criteria (p-value threshold of 5 × 10^-8, an r² threshold of 0.001, and a 1 Mb window size) to ensure robust and reliable findings. LD estimation was based on the 1000 Genomes European data (phase 3). Outlier IVs identified by RadialMR were excluded from the analysis, and any IVs with associations to potential confounders (p < 5 × 10⁻⁸ based on the GWAS Catalog) were also removed. The primary analysis was conducted using the Inverse-Variance Weighted (IVW) method [31]. The MRPRESSO package was employed to address potential horizontal pleiotropy [32]. Sensitivity analyses were performed to evaluate the robustness of the results [33]. Additionally, the MR Steiger method was utilized to assess the true direction of causality [34].

The Summary data-based Mendelian Randomization (SMR) method was employed to investigate the relationship between blood protein levels and mvFibroPsych, to identify potential blood biomarkers [35]. This analysis utilized cis-pQTLs identified from the UK Biobank Pharma Proteomics Project (UKB-PPP), which included plasma samples from 54,219 individuals of European ancestry and measured 2,923 plasma proteins using the Olink Explore platform [36]. Cis-pQTLs were defined as SNPs located within 1 Mb of the transcription start site (TSS) of the gene encoding the corresponding protein. Only index cis-pQTLs associated with protein levels at a genome-wide significance threshold (P < 5 × 10^-8) were included in the SMR analysis. To differentiate between true pleiotropy and linkage, where distinct but linked causal variants could influence protein levels and the disease phenotype, the Heterogeneity in Dependent Instruments (HEIDI) test was conducted. A protein with a p-SMR below the FDR correction threshold and a p-HEIDI exceeding 0.01 was considered to have a true causal relationship with mvFibroPsych, not driven by linkage disequilibrium. Additionally, alternative datasets from the deCODE (N = 35,559), the FENLAND (N = 10,708), and the INTERVAL (N = 3301) datasets were applied as the validation cohort for this analysis [37].

To investigate the influence of gene expression, DNA methylation, and alternative splicing on mvFibroPsych, we applied the SMR method, incorporating eQTL (estimated effective N = 2,443), mQTL (estimated effective N = 1160), and sQTL (estimated effective N = 2,443) data from the BrainMeta data sources [38,39]. First, we used cis-eQTLs to assess the impact of brain-wide gene expressions on mvFibroPsych, where cis-eQTLs were defined as SNPs located within 1 Mb of the TSS of the corresponding gene and significantly associated with gene expression at the genome-wide significance threshold (P < 5 × 10^-8). Subsequently, mQTL and sQTL data were utilized to explore whether DNA methylation or alternative splicing of these genes also exhibited causal effects on mvFibroPsych. The SMR method was used to estimate the causal relationship between these molecular traits and mvFibroPsych, and the HEIDI test was performed to distinguish true pleiotropy from linkage. A significant SMR result (p-SMR below the FDR correction threshold) combined with a HEIDI p-value exceeding 0.01 indicated a likely causal association between the molecular trait and mvFibroPsych, free from confounding by linkage disequilibrium.

After meta-analysis, the estimated effective sample sizes (Neff) were 17,827 for fibromyalgia, 115,173 for insomnia, 391,264 for depression, and 191,496 for anxiety disorder (S1 Table). In particular, meta-analysis modestly increased the number of genome-wide significant loci, supporting improved statistical power. A positive genetic correlation was observed among the four univariate input GWASs with r_g ranging from 0.55 to 0.84, (P-value ranging from 1.23 × 10^-94 to 1.29 × 10^-10) (Fig 2A and S2 Table). Model fit for the common factor shows a good fit to the implied genetic covariance matrix between fibromyalgia, insomnia, depression, and anxiety with a comparative fit index (CFI) = 0.999 and a standardized root mean square residual (SRMR) = 0.015, providing evidence for a shared genetic factor. The standard factor loadings for the four input phenotypes are presented in Fig 2B, with factor loading varying from 0.63 to 0.94, aligned well with prior literature [20].

In the sensitivity analyses, we systematically assessed the robustness of the latent factor model by constructing GenomicSEM models using subsets of the original four phenotypes. Across all three-trait configurations—including fibromyalgia with any two psychiatric traits—the models demonstrated consistently acceptable fit (CFI > 0.95, SRMR < 0.08) (), and fibromyalgia retained statistically significant factor loadings, indicating its stable contribution to the shared genetic architecture. In pairwise models where the factor loading of fibromyalgia was fixed to 1 and the loading of the psychiatric trait was freely estimated, the psychiatric traits also exhibited statistically significant loadings onto the latent factor (). These findings confirm that the observed multivariate structure is not solely driven by highly powered psychiatric traits, and that fibromyalgia contributes non-redundant shared variance, despite its smaller sample size relative to the other phenotypes. S3 Table S3 Table

The common factor model built in the Genomic SEM was applied to incorporate individual variants, enabling the generation of a multivariate GWAS that estimated 7,627,423 associations at the SNP level for the shared factor mvFibroPsych. A total of 49 lead SNPs across 43 independent loci were identified at genome-wide significance (P < 5 × 10 ⁻ ⁸) (Fig 2C and S4 Table). Among these, 32 loci had not been previously reported in any of the eight input GWAS datasets, underscoring the enhanced power and resolution offered by the GenomicSEM approach. To assess the relative contribution of each trait to the genome-wide significant loci, we examined the Z-statistics of lead SNPs across the four input phenotypes. The results revealed that the strongest effects were not consistently driven by depression; in many cases, fibromyalgia, insomnia and anxiety had comparable or greater Z-values (S5 Table). This indicates that the top SNPs identified in the mvFibroPsych GWAS capture shared variance across multiple traits, rather than being solely attributable to the best-powered phenotype. To evaluate whether genome-wide significant loci from the mvFibroPsych GWAS act homogeneously across all component phenotypes, we computed Cochran’s Q statistic for each SNP. Among the significant loci, only one SNP showed evidence of cross-trait heterogeneity (rs1245129, Q_P-value less than the Bonferroni corrected threshold), suggesting trait-specific effects at this locus. Notably, several genome-wide significant loci for the mvFibroPsych latent factor showed modest or non-significant marginal effects in the individual univariate GWASs. For example, the lead variant rs4865477 reached genome-wide significance in the latent factor GWAS (P = 2.31 × 10 ⁻ ²⁴) despite all input Z-scores for the individual traits falling within ±1.96. Importantly, rs4865477 did not show significant Q statistics, indicating that its effect is consistent with a pure common-factor model and reflects shared genetic liability across traits rather than trait-specific heterogeneity. Full results, including Q statistics and heterogeneity-adjusted p-values, are provided in S5 Table. The genomic control (λ_GC) was estimated at 1.207, with an LDSC intercept of 0.812 (se = 0.0074). The effective sample size for the mvFibroPsych GWAS was estimated to be 5,927,502. To evaluate whether the observed inflation in test statistics reflected true polygenicity or residual confounding, we calculated the attenuation ratio using the formula (LDSC intercept − 1)/ (λ_GC − 1), which yielded an attenuation ratio of approximately -0.908. This negative value is atypical in single cohort GWAS but has been observed in multivariate analyses where phenotype sparsity, cohort heterogeneity, or conservative correction procedures may drive intercept deflation [40–42]. We also confirmed that the heritability Z-scores and LDSC cross-trait intercepts showed no evidence of uncontrolled population structure.

The SNP-to-gene mapping was performed using FUMA (v1.6.1). Among all SNPs in LD (r² ≥ 0.6) with independent genome-wide significant SNPs, 54.9% were mapped to genes by positional proximity, 48.1% via eQTL mapping, and 0% via chromatin interaction. Mapping strategies were applied in parallel, and some SNPs were annotated by multiple methods (S6 Table). Based on position mapping and eQTL mapping, gene prioritization identified 342 protein-coding genes associated with mvFibroPsych (S7 Table). The MAGMA gene-set analyses revealed that mvFibroPsych has an enrichment in the pathways involved in synaptic function (S8 Table). The MAGMA tissue expression analysis using the GTEx V8 53 general tissue datasets shows that the brain and the pituitary are the most associated tissues (Fig 2D).

Of the 1,266 phenotypes examined by LCV, a total of 133 phenotypes demonstrated a |GCP| greater than 0.6 and remained significant after FDR correction. Among these, 131 phenotypes showed evidence of genetically causal effects on mvFibroPsych (negative GCP), while mvFibroPsych itself was identified as a potential causal factor for 2 phenotypes (positive GCP) (Fig 3 and S9 Table). The mvFibroPsych showed a positive causal association with “Substances taken for anxiety: Drugs or alcohol (more than once)” (r_g = 0.855, GCP = 0.889, FDR = 3.96 × 10⁻¹⁵). This strong positive genetic correlation indicates that individuals with higher genetic liability captured by mvFibroPsych are more likely to engage in substance use as a coping mechanism for anxiety, aligning with the affective dysregulation underlying the latent construct. Positive causal associations were also observed in diseases of “Back pain for 3+ months” (r_g = 0.385, GCP = 0.749, FDR = 5.74 × 10⁻³). Traits with GCP values below –0.6 showed strong evidence of exerting upstream causal influences on mvFibroPsych. Notably, excessive, frequent and irregular menstruation (ICD10: N92), deep vein thrombosis (DVT), syncope and collapse (ICD10: R55), synovitis and tenosynovitis (ICD10: M65), and urinary tract/kidney infection demonstrated pronounced effects, suggesting that these disorders may act as upstream contributors predisposing individuals to the latent mvFibroPsych factor (r_g > 0, GCP < -0.6, FDR < 0.05) (Fig 3 and S9 Table).

The F-statistics for the selected IVs were all greater than 10, indicating a low likelihood of weak instrument bias. In the analysis of the causal direction from brain imaging-derived phenotypes (IDPs) to mvFibroPsych, 11 out of 227 sMRI phenotypes and 14 out of 360 dMRI phenotypes demonstrated a suggestive causal association with mvFibroPsych (p < 0.05) (S10 and S11 Tables and Fig 4). However, after FDR correction, only one dMRI phenotype—fractional anisotropy (FA) in the splenium of the corpus callosum—remained significantly associated with mvFibroPsych. Specifically, per 10 standard deviation (SD) increase in FA in the splenium of the corpus callosum was associated with an 8.1% reduction in the odds of mvFibroPsych (OR = 0.919, 95% CI 0.880–0.960, FDR = 0.048). FA measures the degree of anisotropic diffusion —that is, the extent to which water molecules diffuse more readily along the direction of white matter fibers than perpendicular to them. FA values typically decrease in structurally compromised tissue, such as damaged or demyelinated white matter tracts. The observed negative causal association between FA in the splenium of the corpus callosum and mvFibroPsych suggests that the integrity of this brain region may play a protective role against the development of mvFibroPsych. The MR Steiger test confirmed the directionality of the causal relationship. Furthermore, sensitivity analyses, including MR-Egger, MRPRESSO, and leave-one-out tests, demonstrated the robustness of the MR estimates.

No significant associations were found in the bidirectional Mendelian Randomization analysis examining the potential causal effects of mvFibroPsych on brain IDPs (and). This lack of evidence suggests that brain changes are more likely to be a causal factor in mvFibroPsych rather than a consequence of the condition, supporting the hypothesis that neuroanatomical alterations may predispose individuals to mvFibroPsych rather than result from it. S12 S13 Tables

Of the 2,923 proteins examined in the UKB-PPP datasets, 1,987 had sufficient instrumental variables to yield reliable effect estimates, and 25 proteins met the FDR-corrected significance threshold (P < 5.29 × 10^-4, FDR < 0.05). Among these, five proteins—UBE2L6, VWC2, BTN3A2, FES, and HSPA1A—exhibited a P_HEIDI value greater than 0.01, suggesting that their association with mvFibroPsych is not due to linkage disequilibrium but rather reflects a true causal relationship (S14 Table). Notably, all five proteins showed a positive causal association with mvFibroPsych, underscoring their potential relevance as biomarkers for this condition. Remarkably, VWC2 demonstrated a significant positive association with the occurrence of mvFibroPsych, consistently passing the FDR correction threshold in the UKB-PPP and all the validation datasets including the deCODE, the FENLAND, and the INTERVAL datasets (Fig 5 and S15-17 Tables). The effect directions were consistent, and the P_HEIDI value exceeded 0.01, indicating a robust causal relationship, prioritizing the potential importance of VWC2 as a biomarker for mvFibroPsych.

Using the SMR-HEIDI method, we identified 20 genes whose brain-wide expression had significant associations with mvFibroPsych (FDR < 0.05, p-HEIDI>0.01) (S18 Table). Among these, the top three genes were SLC12A5 (beta = 0.0044, FDR = 1.83 × 10^-3, p-HEIDI > 0.01), BPTF (beta = -0.0248, FDR = 6.45 × 10^-3, p-HEIDI > 0.01), and FURIN(beta = -0.0039, FDR = 7.78 × 10^-3, p-HEIDI > 0.01), indicating strong and statistically significant associations between their expression levels and mvFibroPsych. The methylation QTL (mQTL) analysis further identified four genes—AMZ1 (cg23627948, beta = 0.001, FDR = 2.1 × 10^-2, p-HEIDI > 0.01), CD40 (cg09053081, beta = 0.0048, FDR = 2.2 × 10^-2, p-HEIDI > 0.01), GSDME (cg27436603, beta = -0.0087, FDR = 2.0 × 10^-2, p-HEIDI > 0.01), and TMEM258 (cg18171955, beta = -0.0033, FDR = 3.6 × 10^-2, p-HEIDI > 0.01)—where methylation was significantly associated with mvFibroPsych (S19 Table). In the splicing QTL (sQTL) analysis, five genes—AREL1, ARHGAP19, BPTF, CD40, and OTOA—exhibited significant splicing variants associated with mvFibroPsych (S20 Table). Notably, CD40 once again demonstrated significant associations at splicing levels (splicing beta = 0.0103, FDR = 1.8 × 10^-2, p-HEIDI > 0.01), underscoring its potential role as a key regulatory factor in the pathology of mvFibroPsych (Fig 6). Interestingly, the functional roles of these genes align with known pathways implicated in neuroinflammation, synaptic plasticity, and cellular apoptosis, which are hypothesized to be relevant to the pathogenesis of fibromyalgia and associated psychological conditions. The association of genes such as BPTF and ARHGAP19, which are involved in chromatin remodeling and cell signaling, respectively, provides further evidence that mvFibroPsych may be driven by dysregulated molecular processes that affect both neurodevelopment and neurodegeneration [43,44].