Druggable Genome-Wide Mendelian Randomization Identifies Distinct Therapeutic Targets for Low Back Pain, Intervertebral Disk Degeneration, and Sciatica

Introduction

Low back pain (LBP) refers to pain occurring in the area from the twelfth rib down to the iliac crests and is one of the most common causes of loss of working capacity, constituting a significant global health issue.¹ The causes of LBP are complex and generally categorized into specific and non-specific types. Although the vast majority of cases are classified as non-specific, intervertebral disc degeneration (IVDD), as a common degenerative spinal condition, is a key contributor to its onset.² When the intervertebral disc degenerates, the nucleus pulposus protrudes posteriorly from the site of rupture, leading to compression of the surrounding nerves. When the protruding nucleus stimulates the nerve roots in the lower back, it can cause LBP. Additionally, if the L4 or L5 lumbar nerves or the S1 sacral nerve are affected, sciatica may occur. LBP, IVDD, and sciatica can affect individuals of all ages, placing significant strain on global healthcare systems.^3,4 The current treatment for LBP, IVDD, and sciatica faces challenges due to limited understanding of their etiology and pathogenesis. Most medications are primarily aimed at symptom relief, which fails to achieve disease cure or progression delay.⁵

Fortunately, with the rapid advancement of human genetics and multi-omics research, large-scale genome-wide association studies (GWAS) have identified thousands of expression quantitative trait loci (eQTLs) and protein quantitative trait loci (pQTLs) for genes and proteins. Integrating transcriptomics and proteomics for multidimensional exploration offers new avenues for discovering potential biomarkers and therapeutic targets for complex diseases at the gene transcription and protein expression levels.⁶ In addition, research on the druggable genome, which refers to the subset of genes encoding proteins that are known or predicted to interact with drug-like compounds, offers a systematic framework for therapeutic target identification, although its application still faces methodological and translational challenges.⁷ Mendelian randomization (MR) is a genetic epidemiology method that helps reduce bias caused by confounding variables and reverse causation, providing evidence for elucidating causal relationships.⁸ Conducting MR analysis on eQTL/pQTL to identify disease-associated genes or proteins, and subsequently integrating this with the druggable genome, has led to the identification of potential drug targets for challenging diseases such as sarcopenia and Alzheimer’s disease.^9,10

Previous MR studies have primarily focused on non-pharmacological interventions such as diet, lifestyle, and physical activity, while investigations into potential drug development targets remain lacking.^11,12 Therefore, we performed a systematic druggable genome-wide MR to identify therapeutic targets for LBP, IVDD, and sciatica. First, we conduct two-sample MR analysis in conjunction with previous studies on the druggable genome to identify candidate genes. Next, we further validate and confirm potential therapeutic targets through colocalization analysis, SMR analysis, and the Steiger filtering test. Finally, we use phenome-wide MR to assess the potential side effects of the identified genes.

Materials and Methods

Study Design

First, we took the union of data from Drug–Gene Interaction Database (DGIdb) and the results of a previous study,¹³ yielding 5,583 druggable genes. Second, we conducted MR analysis using cis-eQTL and cis-pQTL data from the eQTLGen Consortium and the UKBB Proteomics Project to identify potential targets for LBP, IVDD, and sciatica, intersecting these with previously identified druggable genes. Bayesian colocalization analysis, SMR analysis, and the Steiger filtering test were then applied to further identify the therapeutic targets. Finally, a Phe-wide MR study was conducted to evaluate the side effects or other indications of the selected therapeutic targets. The study design is shown in Figure 1.

Figure 1 Flowchart of the study design. Abbreviations: DGIdb, Drug–Gene Interaction Database; SMR, summary-data-based Mendelian randomization; UKBB, UK Biobank; eQTL, expression quantitative trait loci; pQTL, protein quantitative trait loci; LBP, low back pain; IVDD, intervertebral disk degeneration.

Data Sources

Identification of Druggable Genes

Druggable genes were obtained from the DGIdb V.4.0, which integrates drug–gene interaction information from literature reports, publications, and more than 30 external datasets. We downloaded the “Categories Data”, including all genes classified under druggable categories and mapped to Entrez gene identifiers.¹⁴ In addition, we incorporated data from a recent comprehensive review on gene druggability to further expand our list of candidate therapeutic targets. Clinical development drugs were collected from pharmaceutical company pipelines, ChEMBL database (release 17), and clinicaltrials.gov. The study selected candidate genes by integrating information on known drug–gene interactions, bioactive compounds, and GWAS associations, ensuring comprehensive coverage of potential therapeutic targets.¹³ The combined dataset of druggable genes used in this study is summarized in Table 1.

Table 1 Information on Data Sources for Analysis

Data Sources for eQTL and pQTL

The eQTL and pQTL data were sourced from the eQTLGen Consortium and the UK Biobank Pharma Proteomics Project (Table 1). We selected cis-eQTL and cis-pQTL as genetic variants that are more relevant to drug development, defined as being located within a 1 Mb range of the gene coding sequence.¹⁰ The eQTLGen Consortium includes 31,684 participants and 19,250 gene transcripts, from which 16,987 genes with cis-eQTLs were identified.¹⁵ The UK Biobank Pharma Proteomics Project comprises 54,219 participants and 2,923 proteins, from which 1,982 proteins with cis-pQTLs were identified.¹⁶

Data Sources for LBP, IVDD, and Sciatica

The summary GWAS data for LBP, IVDD, and sciatica were sourced from the latest R11 version released by the FinnGen Consortium.¹⁷ The diagnostic criteria for LBP are primarily based on ICD-10 M54.5, comprising 37,567 cases and 322,314 controls. For IVDD, the criteria are mainly defined by ICD-10 M51, including 46,205 cases and 322,314 controls. The diagnostic standard for sciatica is based on ICD-10 M54.3, encompassing 23,458 cases and 322,314 controls. Detailed inclusion and exclusion criteria are provided in Table 1.

MR Analysis

A number of quality control measures were taken to satisfy the three assumptions of the genetic instrumental variables (IVs) in the MR analysis.¹⁸ First, to identify IVs directly associated with the exposure, we set FDR<0.05 for cis-eQTLs and P<5×10⁻⁸ for cis-pQTLs. Second, to retain independent IVs, we utilized the 1000 Genomes Project reference panel to remove linked cis-eQTLs with r²<0.1 and within 100 kb,^19,20 as well as linked cis-pQTLs with r²<0.001 and within 10,000 kb.²¹ Third, weak IVs with F-statistics<10 were excluded.²² The summary data used in this study were all derived from individuals of European ancestry, with exposure and outcome data sourced from different countries, ensuring no overlap between the datasets.

Inverse variance weighted (IVW) method was the primary approach used in this study, as it is considered the most efficient and statistically robust method when multiple IVs are present. For exposures with a single IV, the Wald ratio method was employed. Additionally, various methods were used to assess causality, including the simple mode, weighted mode, weighted median, and MR-Egger. In sensitivity analyses, we used the Cochrane Q test to assess heterogeneity and the MR Egger intercept test to detect pleiotropy. We performed Bonferroni correction for multiple testing, considering the MR analysis results significant for cis-eQTLs when P<2.94×10⁻⁶ (0.05/16,987) and for cis-pQTLs when P<2.52×10⁻⁵ (0.05/1,982). All statistical analyses were conducted using R software version 4.3.1, with data processing and analysis performed using the TwoSampleMR package.

Colocalization Analysis

To investigate whether the potential therapeutic targets (trait 1) identified through MR analysis are associated with the same genetic variations in LBP, IVDD, and sciatica (trait 2), we performed Bayesian colocalization analysis using the “coloc” package in R software.²³ Based on the Bayesian model, we assessed the posterior probabilities of five hypotheses (H0, H1, H2, H3, and H4) to determine whether a single variant is shared between trait 1 and trait 2.²⁴ The single nucleotide polymorphism (SNP) range considered was pQTL/eQTL±500 kb. The prior probabilities used for detecting SNP associations were set as follows: p1=1×10⁻⁴ for SNP-trait 1 association, p2=1×10⁻⁴ for SNP-trait 2 association, and p12=1×10⁻⁵ for SNP associations with both trait 1 and trait 2. Each of the five models generates a posterior probability, with their sum equal to 100%. A higher posterior probability indicates greater likelihood of the model hypothesis being true. The posterior probability for H4 (PPH4) is considered to have strong colocalization support if it exceeds 80%.²⁵

SMR Analysis and Steiger Filtering Test

A summary data-based Mendelian randomization (SMR) analysis was further employed as a complementary method to verify causal relationships between therapeutic targets with LBP, IVDD, and sciatica.²⁶ Analyses were conducted using SMR software (SMR v1.3.1), with significance set at P < 0.05.

We applied the Steiger filtering test to assess the causal direction between exposure and outcome. When the r² value of the exposure variable exceeds that of the outcome variable, the SNP effect direction is classified as “TRUE”, indicating a lower likelihood of reverse causation.²⁷

Phenome-Wide MR

To further assess whether the identified therapeutic targets may impact the onset of other diseases, potentially resulting in side effects or additional indications, we performed a phenome-wide MR (Phe-MR analysis) evaluation using the AstraZeneca portal (https:// azphewas.com/). The study utilized 14,706 phenotypes from UK Biobank participants, with P<2×10⁻⁹ considered indicative of a significant association.²⁸

Ethics Approval and Consent to Participate

This study used publicly available GWAS summary statistics. All original studies had ethical approval and informed consent, and in accordance with International Ethical Guidelines for Health-related Research Involving Humans (CIOMS, 2017) guidelines²⁹ and the Measures for Ethical Review of Life Science and Medical Research Involving Human Subjects (China, 2023; Article 32, Items 1 and 2),³⁰ no additional ethical review was required.

Results

Druggable Genome

From the DGIdb database and a previous study, 3,953 and 4,463 druggable genes were obtained, respectively (Supplementary Table S1-S2). A total of 5,883 unique druggable genes with Human Genome Organisation Gene Nomenclature Committee approved names were ultimately identified (Supplementary Table S3).

MR Analysis Identified 36 Candidate Druggable Genes for LBP, IVDD, and Sciatica

In MR analysis, the Wald ratio or IVW method was primarily used. For eQTL, after Bonferroni correction (P<2.94×10⁻⁶), 26 potential genes were found to be associated with LBP (Figure 2A and Supplementary Table S4), 45 with IVDD (Figure 2A and Supplementary Table S5), and 26 with sciatica (Figure 2A and Supplementary Table S6). For pQTL, after Bonferroni correction (P<2.52×10⁻⁵), 2 potential proteins were identified to be associated with LBP (Figure 2B and Supplementary Table S7), 6 with IVDD (Figure 2B and Supplementary Table S8), and 2 with sciatica (Figure 2B and Supplementary Table S9). Subsequently, candidate genes identified through MR analysis were intersected with the previously obtained 5,583 druggable genes. Ultimately, 10 candidate druggable genes (BPTF, GRK4, NFKBIL1, P2RY13, SLC5A6, MANBA, STX1A, FGL2, BTN1A1, APOE) were found to be associated with LBP, 18 genes (AKR1C2, CAPN10, CYP27B1, DLK1, ILKAP, KHK, MEG3, SLC14A1, SLC29A1, GFPT1, STX1A, CGREF1, DAG1, BTN1A1, EIF2AK3, SMAD3, SERPINA1, IGFBP3) with IVDD, and 8 genes (GPX1, MEG3, NT5C, SUMO2, BMP6, STX1A, MST1, DAG1) with sciatic. Results from additional MR methods were consistent with those of the primary approach. In sensitivity analyses, no heterogeneity or horizontal pleiotropy was observed for any gene. Detailed results are presented in Figure 3 and Supplementary Table S10.

Figure 2 The Volcano plot of MR analyses results. (A) eQTL MR analyses results. The horizontal dashed line represents the corrected threshold of 2.94×10−6. Only the top 10 genes with the largest effect size changes are shown. (B) pQTL MR analyses results. The horizontal dashed line represents the corrected threshold of 2.52×10−5. The x-axis represents the Beta value, the y-axis represents the -log10(P value) of the MR result. Abbreviations: MR, Mendelian randomization; eQTL, expression quantitative trait loci; pQTL, protein quantitative trait loci; LBP, low back pain; IVDD, intervertebral disk degeneration.

Figure 3 Candidate druggable genes for LBP, IVDD, and sciatica. Abbreviations: LBP, low back pain; IVDD, intervertebral disk degeneration; IVW, Inverse variance weighted; SNP, single nucleotide polymorphism.

Colocalization Analysis, SMR Analysis and Steiger Filtering Determined 9 Therapeutic Targets for LBP, IVDD, and Sciatica

We conducted a Bayesian colocalization analysis to determine whether candidate druggable genes share the same genetic variations with LBP, IVDD, and sciatica (Table 2). The results indicated high colocalization support for the association between P2RY13 (PPH4=96.3%) and LBP. CAPN10 (PPH4=91.3%), AKR1C2 (PPH4=98.5%), BTN1A1 (PPH4=83.7%), and EIF2AK3 (PPH4=96.8%) showed high colocalization support with IVDD. NT5C (PPH4=96.4%), GPX1 (PPH4=85.8%), SUMO2 (PPH4=95.9%), and DAG1 (PPH4=86.9%) demonstrated high colocalization support with sciatica.

Table 2 Summary Results for MR Analysis, Colocalization, SMR, and Steiger Filtering Test

We performed SMR analysis to further validate the causal relationship between druggable genes with high colocalization support and the outcomes (Table 2). The results demonstrated that all genes passed the test (P_SMR<0.05). To verify the direction of causality, we also conducted a Steiger filtering test, which revealed that the effect direction of all druggable genes was “TRUE”, indicating no reverse causality (Table 2).

Based on comprehensive examination and validation, we identified nine potential druggable targets. Specifically, increased transcription levels of P2RY13 (OR=1.09, 95% CI: 1.05 to 1.12) were associated with a higher risk of LBP. Increased gene prediction levels of CAPN10 (OR=1.22, 95% CI: 1.13 to 1.32) were associated with a higher risk of IVDD, while elevated gene prediction levels of AKR1C2 (OR=0.94, 95% CI: 0.92 to 0.96), BTN1A1 (OR=0.85, 95% CI: 0.80 to 0.91), and EIF2AK3 (OR=0.75, 95% CI: 0.67 to 0.85) were associated with a reduced risk of IVDD. High levels of NT5C (OR=1.15, 95% CI: 1.10 to 1.21) and GPX1 (OR=1.53, 95% CI: 1.31 to 1.80) were associated with an increased risk of sciatica, whereas elevated gene prediction levels of SUMO2 (OR=0.81, 95% CI: 0.74 to 0.87) and DAG1 (OR=0.46, 95% CI: 0.33 to 0.63) were linked to a decreased risk of sciatica (Table 2).

Phenome-Wide MR Results

Phenome-wide MR analysis indicated that inhibiting NT5C expression may have adverse cardiometabolic effects (Supplementary Table S11). No other therapeutic targets were found to have potential side effects or additional indications.

Discussion

In this eQTL and pQTL-based study, we used MR analysis to identify 10 candidate druggable genes associated with LBP, 18 candidate druggable genes with IVDD, and 8 candidate druggable genes with sciatica. Through a series of validations including Bayesian colocalization, SMR analysis, and Steiger filtering, we identified one therapeutic target (P2RY13) for LBP, four therapeutic targets (CAPN10, AKR1C2, BTN1A1, EIF2AK3) for IVDD, and four therapeutic targets (NT5C, GPX1, SUMO2, DAG1) for sciatica with compelling evidence. Phenome-wide MR analysis revealed a potential adverse cardiometabolic effect associated with NT5C.

Although LBP has become one of the most common causes of disability, its specific etiology is challenging to determine.³¹ Recent evidence indicates that adding educational interventions to physiotherapy provides no additional benefit for chronic non-specific LBP,³² highlighting the limitations of current non-pharmacological strategies and emphasizing the need to explore novel therapeutic targets. This study identified P2RY13 (purinergic receptor P2Y13) as a potential therapeutic target. The product of the P2RY13 gene is the P2Y13 receptor, which belongs to the P2Y receptor subgroup of G protein-coupled receptors. It plays a significant role in regulating cholesterol and glucose metabolism, bone homeostasis, and functions within the nervous system, such as neurotransmitter release, pain transmission, and neuroprotection.³³ Therefore, the P2Y13 receptor may be involved in the stimulation of various elements, such as soft tissues, vertebrae, facet joints, sacroiliac joints, intervertebral discs, and neurovascular structures, by noxious, neuropathic, and harmful or nonspecific factors, thereby mediating the onset and progression of LBP.

The pathogenesis of IVDD is complex, influenced by various factors including environmental and genetic components, and involves intricate biological and biomechanical processes.³⁴ Through multiple validation methods, this study identified BTN1A1, EIF2AK3, CAPN10, and AKR1C2 as potential therapeutic targets for IVDD. BTN1A1 (butyrophilin 1A1) belongs to the butyrophilin (BTN) family, which primarily regulates T cell immune responses.³⁵ The intervertebral disc is an immune-privileged organ that becomes a target for immune system attack when injured, triggering an immune response. Activated T cells can regulate the inflammatory response, thereby mediating the occurrence and progression of IVDD.³⁶ Animal studies have indicated that BTN1A1 has a regulatory role in immune inflammation, capable of inhibiting T cell activation.³⁷ The BTN1A1 protein may mediate a more direct effect of immune inflammation on IVDD by regulating T cell activation. EIF2AK3 (eukaryotic translation initiation factor 2 alpha kinase 3) is a transmembrane protein located in the endoplasmic reticulum membrane, primarily responsible for mediating endoplasmic reticulum stress (ERS) signaling.³⁸ Numerous risk factors, such as genetics, aging, abnormal mechanical load, obesity, and diabetes, contribute to the pathogenesis and progression of IVDD by inducing ERS, which mediates processes like nucleus pulposus degeneration, annulus fibrosus rupture, and cartilaginous endplate degeneration.³⁹ Studies have shown that inflammatory factors, such as IL-1β and TNF-α, mediate damage to nucleus pulposus cells by activating the EIF2AK3 and IRE1-α pathways. Blocking these two pathways can restore the functionality of nucleus pulposus cells, suggesting that endoplasmic reticulum stress, specifically through EIF2AK3 and IRE1-α pathways, may represent potential molecular targets for the treatment of IVDD.⁴⁰ CAPN10 (calpain 10) is a member of the calcium-dependent intracellular cysteine protease family, playing a significant role in glucose metabolism and other cellular processes.⁴¹ Emerging perspectives suggest that immunometabolic changes play a crucial role in IVDD. Conditions such as type 2 diabetes and obesity synergistically activate a pro-inflammatory microenvironment, promoting insulin resistance and hyperglycemia. These changes can activate various catabolic pathways, leading to extracellular matrix degradation, reduced cell viability, and metabolic dysfunction at the cellular level, thereby contributing to the development of IVDD.⁴² Therefore, CAPN10 may be involved in the degeneration of intervertebral disc structure and function by influencing metabolic homeostasis. The AKR1C2 (aldo-keto reductase family 1 member C2) gene encodes a member of the aldehyde/ketone reductase superfamily. Research has identified and validated iron metabolism-related biomarkers in IVDD and suggests that AKR1C2 may slow the progression of IVDD by degrading lipid peroxides to counteract ferroptosis.⁴³

Sciatica is a clinical syndrome caused by inflammation, injury, or compression of the sciatic nerve, which is composed of lumbar and sacral nerve roots (L4-S1). It is characterized by radiating leg pain, with or without associated neurological deficits related to the affected nerve. The majority of cases are secondary to lesions in the local and surrounding structures of the sciatic nerve, such as lumbar disc herniation and lumbar osteoarthritis.⁴⁴ This study identified 4 potential drug targets (NT5C, GPX1, SUMO2, DAG1) for the treatment of sciatica. The NT5C (5’, 3’-nucleotidase, cytosolic) gene encodes a nucleotidase that belongs to the enzyme family responsible for dephosphorylating non-cyclic nucleoside monophosphates into nucleosides and inorganic phosphate, playing a role in maintaining cellular nucleotide and nucleoside levels. NT5C may mediate glucose uptake and lipid oxidation by regulating the intracellular AMP/ATP ratio, contributing to the degeneration of intervertebral disc structure and function, and consequently promoting the occurrence and progression of sciatica.⁴⁵ The GPX1 (glutathione peroxidase 1) gene encodes a protein that catalyzes the reduction of organic hydroperoxides and hydrogen peroxide using glutathione, thereby protecting cells from oxidative damage.⁴⁶ Oxidative stress plays a significant role in intervertebral disc degeneration, and the expression of GPX1 may help delay the progression of disc degeneration by alleviating oxidative stress-induced cellular damage, potentially preventing the onset of sciatica. The SUMO2 (small ubiquitin-like modifier 2) gene encodes a protein that is a member of the SUMO protein family. This protein can modify other proteins in a manner similar to ubiquitin, and it is associated with various cellular activities.⁴⁷ SUMOylation and deSUMOylation regulate ion channels that mediate pain signaling, potentially contributing to the symptoms of sciatica.⁴⁸ DAG1 (dystroglycan 1) is a key component of the dystrophin glycoprotein complex (DGC), widely expressed in muscle and both the central and peripheral nervous systems. It mediates the interaction between the extracellular environment and the actin cytoskeleton.⁴⁹ Research has shown that dystroglycan plays a crucial role in the formation and/or maintenance of peripheral nerve myelin, as well as in stabilizing Na+ channels. DAG1 may be involved in the development and progression of sciatica by mediating extracellular interactions between neuronal cell subpopulations and the extracellular matrix, as well as intracellular interactions.⁵⁰

Our study has several strengths. First, we utilized large-scale eQTL and pQTL data to identify cis-SNPs more closely associated with gene regulatory mechanisms and signaling pathways. This enabled us to perform Mendelian randomization (MR) analyses at both gene transcription and protein expression levels. Additionally, the use of the latest and largest GWAS datasets for LBP, IVDD, and sciatica enhances the credibility of our findings. Second, we integrated previously validated drug target data to increase the success rate of drug development. Third, our findings were rigorously validated through Bayesian colocalization analysis, SMR analysis, and the Steiger filtering test, ensuring the stability of the results. Finally, we confirmed the safety of therapeutic targets through Phenome-wide MR analysis, providing valuable insights for future drug development.

Our study also has limitations. First, the number of IVs for eQTLs and pQTLs was relatively small due to stringent selection criteria, with some analyses based on only a single SNP. This limited our ability to perform sensitivity analyses, such as heterogeneity and pleiotropy assessments, which may affect the stability of the MR results, despite attempts at validation using other methods. Second, all data were derived from European populations, limiting the generalizability of our findings to other ethnic groups. Third, the Phenome-wide MR analysis was based solely on the UK Biobank, so additional side effects may remain undetected. Fourth, the restricted number of IVs in some MR analyses may limit the confidence in the results. Fifth, our phenome-wide MR analyses were restricted to UK Biobank data, which may not capture the full spectrum of potential adverse effects. Lastly, the modest effect sizes of the identified targets are consistent with the polygenic nature of LBP, IVDD, and sciatica, and may partly reflect weak instruments. Therefore, these MR findings should be considered hypothesis-generating. Despite multiple validation approaches, residual horizontal pleiotropy and unmeasured confounding cannot be entirely excluded.

Conclusion

Our study identifies one promising therapeutic target (P2RY13) for LBP, four (CAPN10, AKR1C2, BTN1A1, EIF2AK3) for IVDD, and four (NT5C, GPX1, SUMO2, DAG1) for sciatica. These findings prioritize candidate targets for future experimental validation and provide a framework to guide drug discovery efforts for these common and debilitating conditions. Rigorous clinical trials will be required to confirm their therapeutic potential and translate these genetic insights into effective treatments.