Investigating The Bidirectional Relationship Between Low Back Pain, Hip Pain, Knee Pain, And Gait Abnormalities: A Mendelian Randomization Study

Summary

IVW analyses suggest low back and hip pain may be associated with gait abnormalities, whereas knee pain and reverse associations remain inconclusive.

Abstract

Gait abnormalities are associated with altered lower-limb biomechanics and functional impairment in musculoskeletal disorders. Although low back pain, hip pain, and knee pain have been linked to impaired gait in observational studies, their causal relationships remain unclear. We conducted a bidirectional Mendelian randomization (MR) study to assess potential causal associations between site-specific pain and gait abnormalities. Independent single nucleotide polymorphisms from large-scale genome-wide association studies were used as genetic instruments for low back pain, hip pain, and knee pain. Gait abnormalities, defined in the source GWAS as self-reported difficulty walking, were used as the outcome in forward MR analyses and as the exposure in reverse MR analyses. The inverse-variance weighted (IVW) method was the primary analysis, supported by sensitivity analyses for heterogeneity, horizontal pleiotropy, and robustness. In forward analyses, genetically predicted low back pain (OR = 1.53, 95% CI = 1.058–2.219, P = 0.024) and hip pain (OR = 1.55, 95% CI = 1.067–2.252, P = 0.021) were associated with increased risk of gait abnormalities in IVW analyses. However, the four additional MR methods were not statistically significant, although estimates were generally directionally consistent. Genetically predicted knee pain showed no significant association (OR = 2.15, 95% CI = 0.234–19.789, P = 0.498), with substantial uncertainty reflected by the wide confidence interval. Reverse MR analyses found no evidence that genetic liability to gait abnormalities increased the risk of low back, hip, or knee pain, although only four gait-abnormality SNPs were available. No substantial heterogeneity or horizontal pleiotropy was detected. Overall, this study provides limited IVW-based genetic evidence that low back pain and hip pain may be associated with gait abnormalities. Findings should be interpreted cautiously and validated using larger GWAS datasets with refined phenotypes.

Introduction

Gait abnormalities can stem from motor or sensory impairments, and their clinical features are contingent upon the location and nature of the underlying pathology. Specific abnormal gait patterns may provide important diagnostic clues for particular diseases. Notably, in musculoskeletal disorders, gait abnormalities often present with distinctive biomechanical characteristics and are closely related to pain, joint dysfunction, and reduced mobility^1,2,3. Studies have shown that long-term gait abnormalities can aggravate the condition of patients with osteoarthritis⁴. Moreover, persistent abnormal gait may increase instability of the locomotor joints and alter joint loading; repetitive joint activities under abnormal mechanics may further exacerbate joint injury, thereby inducing pain, promoting gradual muscle atrophy, and resulting in more pronounced muscle strength loss. These factors may interact to form a vicious cycle, contributing to the progressive worsening of pain and functional impairment in patients with musculoskeletal disorders⁵. Therefore, clarifying the direction and potential causality of associations between site-specific musculoskeletal pain and gait abnormalities may help interpret observational pain–gait relationships and inform future research on gait dysfunction and rehabilitation.

However, despite accumulating observational evidence linking site-specific pain—such as low back pain, hip pain, and knee pain—to impaired gait performance, conventional observational studies are vulnerable to confounding, measurement error, and reverse causality, making it difficult to determine the direction and causality of these associations. In addition, the relationship between pain and gait may involve central nervous system adaptation and modulation of nociceptive information received from the periphery, which further complicates the interpretation of observed pain–gait associations. Mendelian randomization (MR) is an epidemiological approach designed to overcome key limitations of observational studies and has been widely applied in causal inference research. MR uses independent single nucleotide polymorphisms (SNPs) as instrumental variables to infer potential causal relationships between exposures and outcomes⁶. For example, rather than assigning individuals to develop low back pain, MR can examine whether genetic liability to low back pain is associated with gait-related outcomes in large Genome-wide association study (GWAS) datasets. Because genetic variants are randomly allocated at conception, meaning at the formation of the embryo, MR can reduce bias from confounders such as age, body mass index, physical activity, and comorbid musculoskeletal conditions. This approach is suitable for exploring population-level causal relationships when well-powered GWAS data and valid genetic instruments are available. However, the phenotypes used in this study were broad self-reported GWAS traits and may not fully capture pain severity, symptom duration, clinician-confirmed diagnoses, or objective gait parameters. By leveraging the random allocation of genetic variants at conception, MR can effectively mitigate bias from confounders and reverse causality and thus yield more plausible causal evidence.

The objective of this study was to perform a bidirectional MR analysis using GWAS summary data on self-reported low back pain, hip pain, and knee pain, together with GWAS data on gait abnormalities from the IEU Open GWAS database. In this study, the pain phenotypes represented broad self-reported site-specific pain traits, and gait abnormalities represented self-reported difficulty walking rather than clinician-confirmed diagnoses, laboratory-measured gait parameters, or standardized functional scores. This bidirectional design aimed to evaluate the direction of potential causal associations between site-specific musculoskeletal pain and gait abnormalities using GWAS summary statistics⁷.

Protocol

This study used publicly available genome-wide association study (GWAS) summary-level data from the IEU OpenGWAS database and GWAS Catalog. No individual-level participant data were accessed, and no new human participants were recruited. Ethical approval and written informed consent had been obtained in the original contributing GWAS studies. Therefore, additional institutional review board approval was not required for the present secondary analysis of publicly available summary statistics.

Study principle and design

The basic steps of MR study include obtaining GWAS summary data, screening and evaluating SNPs, conducting statistical analysis, and implementing quality control measures. As shown in Figure 1, the accuracy of the MR analysis is based on the fulfillment of three crucial assumptions: (1) the relevance assumption: IVs need to be relevant to the exposure phenotypes, namely low back pain, hip pain, and knee pain⁸; (2) the independence assumption: IVs are irrelevant to confounding factors affecting "exposure-outcome"; (3) the exclusivity assumption: IVs affect the outcome only through exposure and not through other pathways⁹. Exposure factors and outcome variables will be interchanged in each analysis to determine if reverse causality exists between the two.

To improve reproducibility, the complete computational workflow used for local GWAS summary-statistic import, SNP extraction, LD clumping, outcome SNP extraction, allele harmonization, MR analysis, sensitivity analysis, diagnostic plotting, and output saving is provided as Supplementary Code 1. The main text describes the key protocol steps, whereas Supplementary File 1 provides the corresponding command-level and function-level R implementation.

Data sources

Data related to low back pain, hip pain, and knee pain are obtained from the IEU openGWAS database, and the website is: https://gwas.mrcieu.ac.uk/. Specific questions used to define these conditions can be found in Supplementary File 2. Low back pain, hip pain, and knee pain were defined according to self-reported site-specific pain items in the original GWAS datasets. The low back pain dataset GWAS ID number is: ebi-a-GCST90018797 , with a sample size of 468269, which includes 22,413 patients with low back pain, and 445,856 controls in the population, with a total of SNP number of 24174741. The hip pain dataset GWAS ID number is: ebi-a-GCST90013968, with a sample size of 407746 and a total SNP count of 11039206. The knee pain dataset GWAS ID number is: ukb-b-16254, with a sample size of 461857, which includes 98704 patients with knee pain, as well as a control population of 363153, with a total number of SNPs of 9851867. GWAS data for gait abnormalities is also obtained from the IEU openGWAS database, Gait abnormalities dataset GWAS ID number is: finn-b-R18 _ABNORMALITI_GAIT_MOBIL, with a sample size of 210717, including 1348 cases of gait abnormalities and 209369 control populations, and a total of 1638044 SNPs loci were identified. Details are shown in Table 1.

For reproducible data acquisition, researchers should search each GWAS ID listed above in the IEU OpenGWAS database, download the corresponding GWAS summary-statistic file where available, and save the file as a local text or comma-separated file before importing it into R. In the computational workflow used in this study, local GWAS summary-statistic files were imported into R using read.table() or read.csv(), depending on the source-file format. The local-file-based workflow was adopted to ensure that the same downloaded GWAS summary-statistic files could be repeatedly processed using fixed column mappings and identical analysis commands. The exact R commands for importing and processing the local files are provided in Supplementary File 1.

Screening and entry of SNPs

SNPs significantly associated with each exposure phenotype were extracted from the GWAS summary datasets using a genome-wide significance threshold of P < 5 × 10⁻⁸. If fewer than three independent SNPs were available at this threshold, the threshold was relaxed to P < 5 × 10⁻⁶ to obtain sufficient genetic instruments for MR analysis. Linkage disequilibrium was assessed using r² < 0.001 within a 10,000 kb window to ensure the independence of instrumental variables.

At the script level, exposure GWAS summary statistics were first imported into R and filtered according to the prespecified P-value threshold. The selected exposure SNPs were then saved as an exposure file and formatted according to the requirements of the TwoSampleMR package. The required column mappings included the SNP identifier, effect estimate, standard error, effect allele, other allele, and P value. In the reproducible R workflow, exposure data were imported using read_exposure_data(), and LD clumping was performed by setting clump = TRUE. The corresponding R command structure was: read_exposure_data(filename = "exposure.csv", sep = ",", snp_col = "rsids", beta_col = "beta", se_col = "sebeta", effect_allele_col = "alt", other_allele_col = "ref", pval_col = "pval", clump = TRUE). Column names were adapted according to the actual downloaded GWAS files.

Instrument strength was assessed using the F statistic:

F = [(N − k − 1)/k] × [R²/(1 − R²)],

where N represents the sample size, k represents the number of instrumental variables, and R² represents the proportion of exposure variance explained by the SNPs. R² was calculated as

Σ[2 × MAF × (1 − MAF) × β²/(SE² × N)],

where MAF is the minor allele frequency, β is the allele effect estimate, and SE is the standard error.

SNPs with F statistics > 10 were retained. SNPs associated with potential confounders, including congenital anatomical abnormalities and body mass index, were identified using PhenoScanner V2 and removed.

Outcome SNPs were extracted by merging the clumped exposure SNPs with the corresponding local outcome GWAS summary-statistic file according to the SNP identifier. If the outcome GWAS reported −log10(P), the P value was converted using P = 10^(−LP). The extracted outcome SNPs were then imported using read_outcome_data(). Exposure and outcome datasets were harmonized using harmonise_data() to align effect alleles and remove SNPs unsuitable for MR analysis. Palindromic SNPs with ambiguous allele orientation were excluded, and SNPs with mr_keep = TRUE were retained for subsequent MR analyses¹¹. The detailed script for outcome SNP extraction, P-value conversion, and harmonization is provided in Supplementary File 1.

Statistical analysis

The causal effects between exposure and outcome variables were estimated in this study using five methods: the IVW method, the MR-Egger regression method, the weighted median method, the weighted mode method, and the simple mode method¹². The IVW method is considered the standard method for MR analysis, which assumes that all instrumental variables are valid, and is based on the principle of combining the Wald ratio estimates, and a random effects model is used if there is heterogeneity and a fixed effects model is used if there is no heterogeneity. The MR-Egger regression method detects potential multicollinearity and takes into account the presence of an intercept term in the regression¹³. The weighted median method requires more than 50% of the instrumental variables to be valid SNPs. The weighted mode method requires smaller sample sizes than other methods and ensures less bias and lower type I error rates. The simple mode method allows SNPs with similar effects to be grouped according to whether the causal effects being estimated are similar or not¹⁴.

At the function level, MR analyses were performed using the mr() function in the TwoSampleMR package. The method list used in the analysis was specified as c("mr_ivw", "mr_egger_regression", "mr_weighted_median", "mr_simple_mode", "mr_weighted_mode"). Odds ratios and 95% confidence intervals were generated using generate_odds_ratios(). The full R implementation, including the exact method list and output-saving commands, is provided in Supplementary File 1.

Quality control and sensitivity analysis

In order to test the stability and reliability of the MR results, first, Cochran Q test was used to assess the heterogeneity of SNPs, if Cochran Q test was statistically significant then it indicated that there was a significant heterogeneity in the analyzed results¹⁵. Second, the MR-Egger intercept test was used to assess potential horizontal pleiotropy. A statistically significant MR-Egger intercept indicates the presence of directional horizontal pleiotropy in the MR analysis. Third, Mendelian randomization pleiotropy residual sum and outlier (MR-PRESSO) was applied to search for the presence of outlier SNPs in the results, and if they existed, they were eliminated and reanalyzed¹⁶. Fourth, the "leave-one-out" method was applied to test the robustness of the results, and the combined effect of the remaining SNPs was calculated by excluding SNPs one by one to assess the effect of individual SNPs on the association between exposure and outcome variables¹⁷.

At the function level, heterogeneity was assessed using mr_heterogeneity(), directional pleiotropy was assessed using mr_pleiotropy_test(), single-SNP effects were assessed using mr_singlesnp(), and leave-one-out analysis was performed using mr_leaveoneout(). Diagnostic plots were generated using mr_scatter_plot(), mr_funnel_plot(), and mr_leaveoneout_plot(). MR-PRESSO was performed when the number of available instrumental variables was sufficient; if MR-PRESSO could not be performed because of an insufficient number of SNPs, this was recorded and reported. The MR analysis and quality control procedures were performed using R version 4.3.2 and the TwoSampleMR package version 0.5.8, with a significance level of α = 0.05. The complete commands are provided in Supplementary File 1.

Output saving, quality checking, and reporting

All intermediate and final outputs were saved for reproducibility. These included the selected exposure SNPs, clumped exposure instruments, extracted outcome SNPs, harmonized datasets, retained SNP lists, MR estimates, odds ratios, heterogeneity results, MR-Egger intercept results, MR-PRESSO results, single-SNP results, leave-one-out results, diagnostic plots, and R session information. Before reporting, retained SNPs were checked for genome-wide significance, LD independence, allele harmonization status, ambiguous palindromic variants, instrument strength, heterogeneity, directional horizontal pleiotropy, and outlier SNPs. The primary MR results were reported as odds ratios with 95% confidence intervals and P values. Sensitivity-analysis results and diagnostic plots were reported in the supplementary tables and figures. The output-saving and file-checking commands are provided in Supplementary File 1.

Results

GWAS data for three site-specific musculoskeletal pain phenotypes, including low back pain, hip pain, and knee pain, and one gait-abnormality phenotype were analyzed in this study. Detailed information on IVs for each exposure factor can be found in Supplementary File 3.

MR analysis of site-specific musculoskeletal pain on gait abnormalities

In this study, SNPs with linkage disequilibrium and palindromic structure were excluded and ...

Discussion

Gait abnormalities can alter lower-limb joint loading and movement strategies, potentially accelerate 821+_osteoarthritic processes and contributing to functional limitations that markedly impair quality of life. This clinical relevance has drawn increasing attention to the interplay between gait dysfunction and musculoskeletal pain. Observational studies suggest that site-specific musculoskeletal pain may be associated with impaired gait; for example, a cross-sectional study of community-dwelling older adults aged 65 ye...

Materials

List of materials used in this article

Name	Company	Catalog Number	Comments
IEU OpenGWAS database	MRC Integrative Epidemiology Unit, University of Bristol	N/A	Source of GWAS summary statistics for exposure and outcome phenotypes.
MR-PRESSO R package	Marie Verbanck / R package	N/A	Used for outlier detection and assessment of horizontal pleiotropy.
PhenoScanner V2	University of Cambridge	V2	Web resource used to identify SNPs associated with potential confounders.
R	R Foundation for Statistical Computing	Version 4.3.2	Statistical computing environment used for MR analysis and quality control.
TwoSampleMR R package	MRC Integrative Epidemiology Unit, University of Bristol	Version 0.5.8	Used for IVW, MR-Egger, weighted median, weighted mode, and simple mode MR analyses.