In order to comprehensively assess the causal relationship between sleep disturbance and LBP, we selected five traits (insomnia, sleep duration, short sleep duration, long sleep duration, and daytime sleepiness) that can reflect sleep disturbance as genetic variants (Jia et al., 2022; Ni et al., 2022).

Insomnia is a common clinical condition characterized by difficulty initiating or maintaining sleep, accompanied by symptoms such as irritability or fatigue during wakefulness (Buysse, 2013). In the present study, the genetic variants of insomnia were obtained from the genome-wide association study (GWAS) among 336,965 individuals of European ancestry from the UK Biobank.

Other sleep traits such as sleep duration, short sleep duration, and long sleep duration, were described in detail in another GWAS summary statistics involving 128,266 subjects (Jones et al., 2016). Sleep duration was recorded based on self-reported sleep time and subjects who reported more than 18 h of sleep within 24 h were excluded from this study. Then, adjustment for age, gender, and study center was made to obtain the model residuals, and inverse normalizing was performed to ensure a normally distributed phenotype. Descriptions of short sleep duration and long sleep duration were defined based on average sleep duration. Short sleep duration refers people who slept less than 6 h per day on average and 28,980 subjects were included in this group. Long sleep duration reported an average of more than 9 h of sleep per day and 10,102 subjects were included in this group. GWAS summary statistics of daytime sleepiness were derived from the UK Biobank (Wang et al., 2019). A total of 452,071 participants of European genetic ancestry self-reported the frequency of daytime sleepiness using the question: ‘How likely are you to doze off or fall asleep during the daytime when you don’t mean to? (e.g., when working, reading or driving),’ with the answer categories ‘never’ (N = 347,285), ‘sometimes’ (N = 92,794), ‘often’ (N = 11,963), or ‘all of the time’ (N = 29). All participants who had taken any sleep-related medication were excluded. The latest GWAS summary statistics describing LBP (finn-b-M13_LOWBACKPAIN) was obtained by visiting the website¹. This GWAS dataset consisting of 13,178 cases and 164,682 controls from European ancestry was identified in 2021.

Researchers screened the genetic variants that met the conditions based on strict quality control from the GWAS summary statistics of various sleep traits including insomnia, sleep duration, short sleep duration, long sleep duration, and daytime sleepiness. It was worth to emphasize that when performing MR analysis using genetic variants (usually single nucleotide polymorphisms [SNPs]) as instrumental variables (IVs), the IVs also need to satisfy three core assumptions: (1) Genetic variants are strongly associated with exposure factors; (2) Genetic variants are associated with the outcome only through the exposure of interest; and (3) Genetic variants are not associated with other confounders affecting the outcome (Boef et al., 2015). First, the single nucleotide polymorphisms (SNPs) associated with five sleep traits with genome-wide significance (p < 5e-8) were extracted. To obtain more IVs associated with the exposure of interest, relaxed thresholds, setting the maximum threshold to 5e-6, were used. This approach to threshold relaxation has been reported in other studies (Chen et al., 2020; Kwok and Schooling, 2021). Because the existence of linkage disequilibrium (LD) may cause corresponding bias, controlling LD before subsequent analysis was necessary. In this study, independent SNPs were selected by setting r² < 0.001 and window size = 10,000 kb.

To further understand whether selected genetic variants were associated with potential confounders and the outcome, researchers visited PhenoScanner², a website which provides details about genetic variants and phenotype information. We focused on physical and psychosocial factors including obesity, smoking, and mood disorders (e.g., anxiety and depression) associated with LBP. IVs that were significantly associated with the above confounders were eliminated before proceeding. However, because genetic variants were not readily available, this process was not performed under extremely stringent conditions.

Statistical analyses were conducted using the R programming language (version 4.0.5). MR analysis was performed based on the “TwoSampleMR” package (version 0.5.6), and the “MRPRESSO” package (version 1.0) was used to apply MRPRESSO analysis to identified the outliers.

Various MR analysis approaches including MR-Egger, weighted median, inverse-variance weighted (IVW), penalized weighted median and maximum likelihood (ML) were performed to estimate the casual effects between sleep traits (insomnia, sleep duration, short sleep duration, long sleep duration, and daytime sleepiness) and LBP. The unbiased estimate of the causal effect can be obtained by IVW regression due to the non-existence of horizontal pleiotropy (Davies et al., 2019). The results of IVW analysis were considered as the major outcome. IVW was able to combine the effect of individual SNP on the outcome and obtained the ratio estimates (β). The ratio estimates were converted to acquire the corresponding odds ratios (ORs) and 95% confidence intervals (95%CIs). MR-Egger, which performing the directional pleiotropy, test of causal effect and the estimate of the causal effect, was an analytical method for MR using pooled genetic summary data (Burgess and Thompson, 2017). The pleiotropy of genetic variants may lead to the failure of the core assumption of IVs, but the causal effect can still be more accurately calculated when up to 50% of the information comes from invalid IVs in the method of weighted median estimator (Bowden et al., 2016). Both of the MR-Egger regression and weighted median could perform to improve the evaluation of IVW due to the robust estimates they could offer. Other methods such as ML (Hartwig et al., 2017) and penalized weighted median were mainly used to assess the robustness of MR results.

The horizontal pleiotropy of genetic variables is important because the results of MR analysis can be significantly affected and cause instability in the effect estimates. The test of correlation horizontal pleiotropy was mainly evaluated based on MR-Egger intercept and MR-PRESSO analysis; the former estimated the possibility of horizontal pleiotropy by calculating the term of intercept obtained after linear regression analysis. F statistics were calculated to evaluate the strength of the IVs. MR-PRESSO analysis identifies the outliers that may possess the characteristic of horizontal pleiotropy. The number of distributions in MR-PRESSO analysis was set to 1,000, and the robustness of MR analysis results was evaluated by comparing whether the casual relationship was affected before and after the outlier elimination. IVW and MR-Egger regression were applied to test the heterogeneity and Q statistics were calculated to quantitatively evaluate the heterogeneity. If the heterogeneity existed (p < 0.05), then the results of random effect IVW were dominant, otherwise it referred to the results of fixed effect IVW.

To further evaluate whether there was a single SNP with a large contribution of pleiotropy that may cause deviation to the results, leave-one-out analysis was conducted to eliminate SNPs one by one and then re-estimate the causal effect. The forest plot was used to evaluate the effect estimation between the genetic variants and LBP, and MR-Egger regression and IVW to calculate the combined effects. If a correlation between parts of the SNPs and LBP was found, a specific rule determined whether a related SNP should be eliminated. More specifically, if p_exposure > p_outcome, then the SNP was not included in the MR analysis. Funnel plot was used to evaluate the publication bias and applied to assess the potential directional pleiotropy in this study.

Supplementary Tables 1–3 shows the process of IV filtering in detail. SNPs were removed in the following situations. First, in the process of extracting SNPs from the outcome GWAS dataset, parts of the SNPs not found in the outcome dataset were removed. The results obtained from analysis of the four groups were as follows: insomnia-LBP (rs11804386), sleep duration-LBP (rs282086 and rs12537376), short sleep duration-LBP (rs573615914, rs546786239, rs8008258, and rs74500417), and long sleep duration-LBP (rs2387776, rs73196898, rs9915132, and rs4006399). Second, ambiguous SNPs with non-concordant alleles or palindromic SNPs with ambiguous strand were removed (Wu et al., 2020). In the analysis of insomnia-LBP, short sleep duration-LBP, long sleep duration-LBP, and daytime sleepiness-LBP, rs10280045, and rs2644128; rs9474974; rs2973993; rs3803763, rs58460356, rs61696052, rs6557066, rs72831782, and rs9475029, respectively, were eliminated. Third, genetic variants associated with the outcome and confounders were eliminated using PhenoScanner. Only rs73196898 was excluded in the MR analysis of long sleep duration and LBP in the present study.

Different methods were used to evaluate the causal relationship between five sleep characteristics and LBP. To further investigate the effect estimation of LBP on sleep traits, reverse casual inference was also implemented, and Table 1 and Supplementary Table 4 show the MR analysis results in detail. Table 2 shows the details of GWAS summary statistics of LBP and sleep traits.

The causal effect of insomnia on LBP was the focus in the present study. The results showed a causal relationship between insomnia and LBP (IVW: OR = 1.954, 95% CI: 1.119–3.411, p = 0.019; MR-Egger: OR = 6.269, 95% CI: 1.388–28.313, p = 0.025; weighted median: OR = 2.095, 95% CI: 1.024–4.286, p = 0.043; maximum likelihood: OR = 1.998, 95% CI: 1.262–3.166, p = 0.003; penalized weighted median: OR = 1.958, 95% CI; 0.947–4.046, p = 0.070). Weak IVs due to the F statistics value was a low possibility (F = 13.928). The results of reverse causal inference indicated LBP can also have a causal effect on insomnia (IVW: OR = 1.015, 95% CI: 1.004–1.026, p = 0.006; MR-Egger: OR = 1.013, 95% CI: 0.988–1.039, p = 0.317; weighted median: OR = 1.012, 95% CI: 0.996–1.028, p = 0.153; maximum likelihood: OR = 1.016, 95% CI: 1.004–1.027, p = 0.007; penalized weighted median: OR = 1.012, 95% CI: 0.995–1.029, p = 0.162). The results of heterogeneity test are shown in Table 1 and Supplementary Table 5 and the output of heterogeneity analysis of insomnia-LBP indicated possible heterogeneity (MR-Egger: Q statistics = 36.01452, p = 0.07137978; IVW: Q statistics = 39.81181, p = 0.04070194). However, the reverse causal inference result appeared stable (MR Egger: Q statistics = 25.74813, p = 0.105658; IVW: Q statistics = 25.77981, p = 0.1364423). MR-Egger intercept was used to evaluate the horizontal pleiotropy and the results are shown in Table 1 and Supplementary Table 6 (intercept _{(insomnia–LBP)} = −0.01610656, p_{(insomnia–LBP)} = 0.11 70137; intercept_{(LBP–insomnia)} = −0.01610656, p_{(LBP–insomnia)} = 0.1170137). MR-PRESSO analysis showed no obvious outliers.

A causal relationship between sleep duration and LBP was not found (IVW: OR = 1.44, 95% CI: 0.851–2.436, p = 0.175; MR-Egger: OR = 0.449, 95% CI: 0.038–5.263, p = 0.589; weighted median: OR = 1.48, 95% CI: 0.785–2.791, p = 0.226; maximum likelihood: OR = 1.446, 95% CI: 0.85–2.46, p = 0.174; penalized weighted median: OR = 1.48, 95% CI: 0.794–2.757, p = 0.217). A similar result was observed in reverse MR analysis (Supplementary Table 4). Abnormality was not found in the MR-PRESSO analysis, heterogeneity test, or MR-Egger intercept (Supplementary Tables 5–7).

The IVW analysis showed that genetic predisposition to short sleep duration did not increase the risk of LBP (OR = 1.45, 95% CI: 0.575–3.656, p = 0.431) and a similar output was found in reverse MR analysis (OR = 1.004, 95% CI: 0.993–1.016, p = 0.443). Heterogeneity test results and MR-PRESSO analysis in bidirectional causal inference were stable (Supplementary Tables 5, 7) and terms of intercept were −0.01105083 and 0.0004451583, respectively (Supplementary Table 6).

When estimating the causal effect of long sleep duration on LBP, results indicated the heterogeneity cannot be ignored (MR-Egger: Q statistics = 27.01214, p = 0.05789008; IVW: Q statistics = 28.88768, p = 0.04976916; Supplementary Table 5). Therefore, the output of random effect IVW analysis was considered to be more reliable. Next, the genetic predisposition to long sleep duration was shown to not affect the risk of LBP (OR = 1.566, 95% CI: 0.887–2.765, p = 0.122) and the results of reverse causal inference were similar (OR = 0.997, 95% CI: 0.989–1.006, p = 0.545) with terms of intercept −0.008811378 and −0.000102443, respectively (Supplementary Table 6). MR-PRESSO analysis showed no significant outliers (Supplementary Table 7).

The results of the sensitivity analysis are summarized in Supplementary Tables 5, 6; both MR-Egger and IVW analysis showed heterogeneity (MR-Egger: Q statistics = 41.43407, p = 0.03739338; IVW: Q statistics = 41.70184, p = 0.04624392) and random effect IVW analysis showed no causal relationship between daytime sleepiness and LBP. However, the MR-PRESSO analysis identified an outlier (rs4765936). When comparing the causal effect before and after removing this outlier, bias was not found in the adjusted result (raw: causal estimate = −0.3669115, SE = 0.493133, T-stat = −0.7440417, p = 0.4619627; outlier-corrected: causal estimate = −0.1740137, SE = 0.4572936, T-stat = −0.3805295, p = 0.705991; Supplementary Table 7). Reverse causal inference indicated that genetic predisposition to LBP could affect the risk of daytime sleepiness (IVW: OR = 1.011, 95% CI: 1.004–1.017, p = 0.00114; maximum likelihood: OR = 1.011, 95% CI: 1.004–1.017, p = 0.00144). MR-Egger intercept did not show horizontal pleiotropy [Intercept _{(daytime sleepiness–LBP)} = 0.008773135; Intercept _{(LBP–daytime sleepiness)} = 0.001033245]. Outliers were not identified in MR-PRESSO analysis and more detailed results are shown in Supplementary Tables 5–7.

Leave-one-out analysis and funnel plot results are shown in Supplementary Figures 1, 2. Due to the existence of individual, potentially influential SNP, the results should be interpreted with caution. The funnel plot showed no probable directional pleiotropy. Figure 1 shows the individual putative causal effect between insomnia and LBP. The intercepts calculated in every method of MR analysis were close to zero which indicated the possibility of horizontal pleiotropy was low. A significant positive correlation was observed between insomnia and LBP. Figure 2 shows the causal effect estimates between each SNP and the outcome, and the combination of the effect estimates based on IVW and MR-Egger regression. Supplementary Table 3 shows the p_exposure of every SNP was less than p_outcome, thus, additional SNPs were not removed in this procedure. Figure 3 shows the design flow chart for the MR study.