Linkage Disequilibrium Score Regression (LDSR)

LDSR can estimate cross-trait genetic correlations using GWAS summary statistics data [12]. We conducted LDSR analyses with the LD Hub interface [13] and the LDSC package (https://github.com/bulik/ldsc). We restricted our LDSR analyses to well-imputed SNPs (imputation quality score > 0.8) and SNPs with Minor Allele Frequency (MAF) ≥ 5% in the study populations. In addition, because outliers can unduly influence the regression, we also excluded SNPs with extremely large effect sizes (X12 > 80) [14]. The MS GWAS (47,429 MS cases and 68,374 non-MS healthy controls) summary statistics was obtained from the recent MS GWAS meta-analysis of the International Multiple Sclerosis Genetics Consortium (IMSGC) [15]. As for other musculoskeletal traits, we used GWAS summary statistics from the UK Biobank with the LD hub pipeline to estimate genetic heritability and assess correlations with MS.

Mendelian Randomization (MR)

MR is a widely used method using measured variation in genes (genotypes, SNPs) of exposure to examine the causal effect of exposures on disease outcomes in observational studies [16, 17].

Most GWAS summary statistics data were downloaded from the “MRC IEU OpenGWAS database” repository (https://gwas.mrcieu.ac.uk/) or each GWAS consortia.

The information on diagnostic criteria, phenotype definition, sample size, name of consortium, and dataset IDs from the database are available in Supplementary Table 1. The overall study design for MR is displayed in Fig. 1.

Fig. 1

Conceptual study designs and MR analyses with univariable, multivariable, and bidirectional MR models

Exposure data

Instrumental variables (IVs) for MS were obtained from the recent GWAS meta-analysis of the IMSGC [15]. The detailed analysis and results description are provided in the previous IMSGC publication [15]. The estimated MS heritability is ~ 48%. Among the genome-wide significantly associated SNPs (p-values ≤ 5 × 10–8), 200 SNPs are located in the autosomal and non-MHC genomic region. Among these 200 SNPs, we further excluded 91 SNPs in LD with the lead SNPs. A total of 109 genome-wide associated SNPs were selected as IVs for MR analysis. The Effect allele, reference allele, effect allele frequency, and effect size information on these 109 SNPs are summarized in Supplementary Table 2. The minor alleles of these 109 SNPs are all larger than 10% in study populations. The effect size (OR) ranges from 1.06 to 1.23.

Outcome data

Fracture and fall GWAS data used in the analysis were selected from the UK Biobank data, fractured/broken in the last 5 years (ID: ukb-b-13346, n = 44,502 cases and 415,887 controls) [18] and falls last year (ID: ukb-b-2535, n = 461,725). Frailty GWAS data were obtained from a study by Atkins et al. in 2019 (n = 175,226) [19].

Predisposing factors of musculoskeletal outcomes

Known bone and muscle strength-related factors were used as covariates in multivariable MR (MVMR). They are heel-estimated BMD (eBMD) [19], appendicular lean mass, whole body fat-free mass, right handgrip strength, left grip strength [20] vitamin D [21]. These covariates were considered in MVMR to control for their potential confounding effects on musculoskeletal three outcomes.

Statistical analysisUnivariable MR between MS and three musculoskeletal traits as well as predisposing factors

Univariable MR was applied to assess the causal relationship between MS and 3 musculoskeletal outcomes (fracture, falls, and frailty) and 5 predisposing factors. We additionally performed MR between 5 predisposing factors and 3 musculoskeletal outcomes to examine potential underlying pathways between them. Multi-instrument MR strategy was applied. This approach involves conducting single instrumental variable MR analyses for each exposure or risk factor independently, and then the MR results from these multiple instrumental variables were combined through a fixed-effect meta-analysis. This meta-analysis allowed us to estimate the overall causal relationships between the exposures (MS and predisposing factors), and musculoskeletal outcomes.

Multivariable MR (MVMR) between MS and three musculoskeletal outcomes adjusting for predisposing factors

We also applied regression-based MVMR to account for many variants having pleiotropic effects that are associated with multiple predisposing factors. This method provides coefficients or effect sizes that represent the direct causal effects of multiple sclerosis (MS) on three distinct musculoskeletal outcomes while holding other predisposing factors constant. The estimation of causal effects in MVMR involves performing regression analyses of the associations between genetic variants and the outcome of interest on the associations between the same genetic variants and the selected predisposing factors. The intercept in these regression analyses is set to zero, and the weights used are the inverse variances of the associations with the outcome [22]. Predisposing factors that demonstrated significant associations in the univariable MR analysis with each musculoskeletal trait were included in the MVMR analysis.

Bidirectional MR among the three musculoskeletal outcomes

Bidirectional MR was performed among fracture, falls, and frailty to examine their mutual effects and directionality considering their highly correlated nature. In bidirectional MR, instruments for both exposure and outcome are used to evaluate whether the “exposure” variable causes the “outcome” or whether the “outcome” variable causes the “exposure”[23].

Heterogeneity analysis

We employed Cochran’s Q statistic as a diagnostic tool to assess the heterogeneity of IVs. Heterogeneity can indicate a potential violation of the necessary assumptions for MR or the presence of pleiotropic effects, where a single genetic variant influences multiple traits. Cochran’s Q statistic derived from the inverse variance weighted (IVW) estimate should follow a χ2 distribution with degrees of freedom equal to the number of SNPs minus 1. When there is excessive heterogeneity, it suggests that either the modeling assumptions for MR have been violated or that some of the genetic variants are not adhering to the MR assumptions [24]. This is termed ‘horizontal pleiotropy’[25]. When Q statistics showed significant heterogeneity, we applied the MR-PRESSO [26] method to reassess the MR effects after filtering out outlier SNPs that may be responsible for the observed heterogeneity. After removing these outliers, we compared the MR results obtained from the MR-PRESSO adjusted analysis with the fixed-effects MR results derived from the IVW approach.

Sensitivity analysis

The sensitivity analysis including IVW, weighted median, weighted mode, and MR-Egger tests was performed in univariable MR analysis. For MVMR, IVW and MR-Egger tests were evaluated. We also applied the MR-Egger regression test [27] to verify horizontal pleiotropy in each MR. Among the multiple methods, we set the results from IVW and outlier removal results from MR-PRESSO as primary findings.

All the MR analyses were conducted in R (version 4.0.0) and specialized packages, including TwoSampleMR [28], MendelianRandomization [29], MRPRESSO [30], and the MR-Base platform [31]. The causal effect size (beta or OR) with standard error and p-values were presented as appropriately. We applied a false discovery rate (FDR, q) to account for multiple testing corrections. The criteria used to select final instrument variables in each MR were GWAS significance: p < 5e−8, LD clumping: r2 > 0.001 within a 10 MB window, and proxy SNPs in 1000 Genomes EUR, r2 = 0.8 [32].

Colocalization

We performed colocalization analysis to identify genetic variants that were associated with both multiple sclerosis (MS) and falls as potential pleiotropic genetic effects. The R COLOC package (https://chr1swallace.github.io/coloc/index.html) was used [33]. The R COLOC package calculates approximate Bayes Factors, which help quantify the evidence for colocalized genetic variants between two traits. We focused on examining a 500 kb region around each site (250 kb on either side) [34] specifically targeting nine SNPs that displayed a significant association in the MR analysis between MS and falls. The posterior probabilities of H0 (no causal variant), H1 (causal variant for trait 1 only), H2 (causal variant for trait 2 only), H3 (two distinct causal variants), and H4 (one common causal variant) are calculated. Among their output results, PP (Post-probability > 95%) H4 was used to determine the high probability of colocalized genetic variants between 2 traits. Conversely low PP4 (< 50%) means we cannot identify which individual SNP is jointly causal with confidence [35].