Introduction

Intervertebral disc degeneration (IVDD) is a prevalent chronic degenerative disease, with its incidence significantly increasing with age. IVDD manifests as symptoms such as lower back pain, lower limb numbness, and cauda equina syndrome, which severely impact patients’ quality of life.¹ Although current treatment approaches for IVDD, such as physical therapy, medication,² and surgery,³ can alleviate symptoms to a certain extent, there is still a lack of fundamental therapies that effectively halt or reverse the progression of IVDD.

Studies have shown that genetic factors play a crucial role in the pathogenesis of IVDD.⁴ To identify potential susceptibility gene loci for IVDD, researchers have conducted numerous candidate gene studies and genome-wide association analyses (GWAS), identifying several genetic variants significantly associated with IVDD.⁵ However, most of these studies are limited to exploring the statistical association between gene polymorphisms and IVDD phenotypes, making it difficult to determine the causal relationship between the relevant genes and IVDD, ie, whether they are the cause or consequence of IVDD. Moreover, as a complex degenerative disease, the development and progression of IVDD involve pathological changes in various cells and tissue components of the intervertebral disc, such as nucleus pulposus cell apoptosis, inflammatory response, matrix degradation, and metabolic disorders.⁶ While GWAS have identified potential genetic associations, the causality and mechanisms linking these genes to IVDD pathology remain unclear. Therefore, systematically elucidating the genetic pathogenesis of IVDD and revealing the molecular pathways through which key driving genes induce IVDD by regulating the cellular function and microenvironmental homeostasis of intervertebral disc tissue are crucial for deepening our understanding of the complex pathogenesis of IVDD and exploring new prevention and treatment strategies.

In recent years, MR has emerged as a powerful causal inference method that can effectively avoid the influence of reverse causality and unmeasured confounding factors, providing new insights for inferring potential disease causes.⁷ eQTL and pQTL are important data for multi-omics integration, which can reveal genetic variants that explain differences in gene and protein expression levels. The measurement of these markers can reflect an individual’s health status and may provide new insights into the impact of diseases.⁸ Furthermore, mediation MR analysis can also detect whether an intermediate phenotype (eg, inflammatory factors, metabolites) mediates the causal chain between an exposure factor (eg, gene expression) and a disease, thereby revealing the pathogenic mechanisms of the disease.⁹ Thus, integrating MR analysis with multi-omics data can reduce the impact of false positives or false negatives from single-omics data, and is expected to screen out key driving factors of IVDD at the genome-wide level, revealing the cellular and molecular mechanisms mediating its occurrence and development, and providing important clues for elucidating the genetic pathogenesis of IVDD and developing new intervention strategies.

This study aims to elucidate the genetic mechanisms underlying IVDD through a comprehensive genome-wide multi-omics integration approach. We will combine transcriptomic data from lumbar intervertebral disc tissues, single-cell sequencing, and MR analysis to identify key genes involved in IVDD. Initially, MR analysis using blood eQTL will screen candidate genes associated with IVDD, which will be validated through differential expression analysis of publicly available IVDD gene expression profiles. Single-cell sequencing will further enable cell type identification, subgroup clustering, and differential gene expression analysis in tissues with varying degrees of degeneration, revealing the cellular and molecular mechanisms influenced by these genes. Additionally, mediation MR analysis will assess whether inflammatory factors and serum metabolites mediate the relationship between the identified genes and IVDD development. Finally, in vitro experiments involving gene overexpression and silencing will validate the regulatory roles of these key genes in IVDD. This integrated strategy is expected to identify novel molecular markers and therapeutic targets, enhancing risk prediction, early diagnosis, and the precise prevention and treatment of IVDD.

Methods

Data Sources

This study utilized multiple large-scale genomic datasets. Firstly, the whole blood eQTL data we employed was obtained from the eQTLGen consortium, which includes 31,684 blood and peripheral blood mononuclear cell (PBMC) samples from 37 eQTLGen consortium cohorts, covering 19,942 genes.¹⁰ Secondly, the GWAS summary data for 91 plasma inflammatory markers was derived from 14,824 subjects, measured using the Olink platform.¹¹ Additionally, we incorporated the results of a genome-wide association analysis conducted by Chen et al on 1,091 blood metabolites and 309 metabolite ratios.¹² To ensure the validity of the genetic instrumental variables used, we applied stringent selection criteria, only including loci that met the following standards: (a) exhibited genome-wide significance (P < 1×10⁻⁸) in the GWAS; (b) had linkage disequilibrium (LD) r² < 0.001 between loci, with a distance greater than 10,000 kb; (c) F-statistic, F = R²(n-2)/(1- R²), used to measure the strength of each genetic instrument, where n is the effective sample size used for SNP association in the GWAS, retaining strong instrumental variables with F-statistic greater than 10.¹³ Finally, the GWAS summary data for IVDD was obtained from the latest R10 version of the FinnGen consortium (https://www.finngen.fi/fi), comprising 41669 IVDD cases and 294770 control individuals (Supplementary Files 1).

MR Analysis

In this study, we employed the Two-Sample Mendelian Randomization method, utilizing summary data from GWAS, to evaluate the causal relationship between individual gene expression levels in whole blood and IVDD.¹⁴ First, we performed stringent screening of eQTL for each gene, only including genetic instrumental variables that met the following criteria: (1) significantly associated with gene expression levels at the genome-wide level (P<5×10⁻⁸); (2) linkage disequilibrium (LD) r²<0.001 between loci to ensure the independence of instrumental variables; (3) F-statistic>10 to exclude the influence of weak instrumental variables.

For each screened gene, we conducted independent univariate two-sample MR analyses using the “TwoSampleMR” package¹⁵ in R software (version 4.1.3), with its eQTL as the instrumental variable (exposure) and IVDD as the outcome measure. We selected the Inverse Variance Weighted (IVW) method as the primary approach to estimate the causal effect. Simultaneously, we combined multiple sensitivity analysis methods, including MR Egger regression, weighted median method, simple mode method, and weighted mode method, to assess the robustness of the results.^15–18 Additionally, we evaluated the heterogeneity among the genetic instrumental variables for each gene using Cochrane’s Q test. To further exclude the influence of horizontal pleiotropy and outliers, we selected MR-Egger intercept test, MR-PRESSO (Mendelian Randomization Pleiotropy RESidual Sum and Outlier) global test,¹⁹ and Leave-One-Out Analysis as sensitivity analysis methods. The MR-PRESSO test was performed using the “MRPRESSO” R package. The significance threshold for all statistical tests was set at a P-value<0.05.

Transcriptome Sequencing Data Analysis

We downloaded two transcriptome datasets related to IVDD from the Gene Expression Omnibus (GEO) (https://www.ncbi.nlm.nih.gov/geo/) database:²⁰ GSE15227 and GSE23130. Among them, GSE15227 includes 12 mildly degenerated and 3 severely degenerated intervertebral disc samples, while GSE23130 contains 6 mildly degenerated and 17 severely degenerated intervertebral disc samples. Subsequently, we performed annotation, batch effect removal, and data integration on the two datasets. After completing data preprocessing, we used the limma software package to conduct differential expression gene (DEG) analysis on the integrated dataset, comparing the gene expression differences between the severely degenerated group and the mildly degenerated group. We set the screening criteria for DEGs as an adjusted P-value < 0.05 and logFC > 0.58.

Single-Cell RNA Sequencing Data Analysis

We obtained a single-cell RNA sequencing (scRNA-seq) dataset (GSE244889) related to IVDD from the GEO database.²¹ This dataset includes one mildly degenerated and one severely degenerated intervertebral disc sample. To analyze these scRNA-seq data, we employed the Seurat software package²² for quality control and filtering of the raw data. The filtered data underwent normalization and standardization to eliminate technical biases in expression levels between different cells. Subsequently, we performed principal component analysis (PCA) on the filtered data. To visually demonstrate the heterogeneity among different cells, we used UMAP for dimensionality reduction and visualization of the cells. We then applied a graph-based clustering method to divide the cells into different clusters. To identify the specific cell types of each cluster, we utilized the SingleR software package²³ for cell type annotation. We then conducted inter-group differential analysis to verify whether the differentially expressed genes obtained from transcriptomics and MR analyses were also significantly highly expressed at the single-cell level. Subsequently, we performed pseudotime analysis on the single-cell data using the Monocle 3 software package.²⁴ By constructing cell developmental trajectories, we discovered that the chondrocyte clusters with significantly high expression of the aforementioned genes might represent a critical developmental state during the process of IVDD.

Reverse MR Analysis

To conduct subsequent mediation MR analysis and identify potential inflammatory proteins and metabolites, we performed reverse MR analysis between the genes screened by multi-omics and IVDD. We excluded genes with potential causal relationships, using the IVW method as the primary approach, and combined MR Egger, weighted median, simple mode, and weighted mode as supplementary measures to assess the causal effects. Results with p<0.05 were considered statistically significant.

Mediation Analysis

To investigate the potential inflammatory protein and metabolite mediators between the aforementioned genes and IVDD, we employed a two-step approach for mediation analysis. First, we used all eligible serum metabolites and inflammatory factors as exposure factors and IVDD as the outcome measure. We performed univariate two-sample Mendelian randomization analysis using the “TwoSampleMR” software package. We selected the IVW method as the primary approach to estimate causal effects and combined sensitivity analysis methods such as MR Egger, weighted median, simple mode, and weighted mode for supplementation. Additionally, we assessed the heterogeneity among genetic instrumental variables (Cochrane’s Q test) and the potential influence of horizontal pleiotropy (MR-Egger intercept test, MR-PRESSO global test) and outliers (leave-one-out analysis) on the results.

In the second step, we used the inflammatory proteins or serum metabolites with significant causal effects from the first step analysis as outcomes and the candidate genes obtained from multi-omics data screening as exposures. We again performed univariate two-sample Mendelian randomization analysis using the IVW method to identify potential mediators. For the effect size of each mediator, we adopted the product of coefficients method.²⁵ This method involves first estimating the effect of genetic risk on the mediator and then multiplying this effect by the effect of the mediator on IVDD adjusted for genetic risk,²⁶ thereby obtaining the indirect effect (ie, the effect of genetic risk on IVDD through the mediator). When estimating the effect of the mediator on the outcome, we excluded SNPs that significantly influenced the mediator by genetic risk to avoid double counting. The indirect effect divided by the total effect estimates the proportion of the total effect of genetic risk on IVDD explained by each mediator. The standard error was calculated using the delta method (based on the effect estimates from the two-sample Mendelian randomization analysis). We analyzed all mediation effects to screen for inflammatory proteins and serum metabolites that may serve as mediators.

Cell Culture

Human primary nucleus pulposus cells (NPCs) were obtained from ScienCell Research Laboratories (Catalog#4800, USA). Cells were seeded at a density of 2×10⁵ cells/mL in T25 culture flasks and maintained in a humidified incubator at 37°C with 5% CO₂. The complete culture medium consisted of Dulbecco’s Modified Eagle Medium (DMEM; Gibco, Invitrogen, USA) supplemented with 10% fetal bovine serum (FBS; Gibco, Invitrogen, USA) and 1% penicillin-streptomycin (Pen-Strep; Gibco, Invitrogen, USA). The medium was refreshed every 2–3 days to ensure optimal nutrient availability and waste removal. Upon reaching 90–95% confluence, cells were passaged by aspirating the medium, rinsing with phosphate-buffered saline (PBS; Gibco, Invitrogen, USA), and adding 0.25% trypsin-EDTA (Catalog#J8150; Solarbio, China) to detach the cells at 37°C for 2–3 minutes with gentle agitation. The trypsin reaction was neutralized by adding an equal volume of complete culture medium, followed by centrifugation at 1000 rpm for 5 minutes. The resulting cell pellet was resuspended in fresh medium and subcultured at a 1:3 split ratio to maintain continuous cell growth for subsequent experiments.

Transfection Assay

To investigate the roles of CLEC11A and ARTN in NPCs, we conducted small interfering RNA (siRNA) knockdown and plasmid overexpression experiments. NPCs were initially seeded in six-well plates at a density of 2×10⁵ cells per well and cultured until reaching 80–90% confluence. Prior to transfection, the cells underwent serum starvation for 30 minutes to synchronize their cell cycle and enhance transfection efficiency. In the siRNA knockdown experiments, the experimental groups were transfected with 50 nM of either CLEC11A siRNA or ARTN siRNA (GenePharma, Suzhou, China), while the control groups were transfected with 50 nM of scrambled siRNA. For the overexpression experiments, the experimental groups were transfected with 2.5 μg of either CLEC11A or ARTN overexpression plasmids (Unibio, Chongqing, China). These plasmids were constructed by cloning the full-length complementary DNA (cDNA) of human CLEC11A or ARTN into the pcDNA3.1(+) vector; the control groups received 2.5 μg of the empty vector. All transfections were performed using Lipofectamine 3000 reagent (Invitrogen, USA) according to the manufacturer’s instructions. Specifically, siRNA or plasmid DNA was mixed with Lipofectamine 3000 in serum-free medium and incubated at room temperature for 15 minutes to form liposome-nucleic acid complexes. This mixture was then added to the cells and incubated at 37°C for 6 hours. Following transfection, the medium was replaced with fresh growth medium. To evaluate transfection efficiency and the effects of CLEC11A and ARTN on NPCs, cells were harvested at various time points post-transfection: 24 hours for RNA analysis and 48 hours for protein detection.

Quantitative Real-Time Polymerase Chain Reaction (qRT-PCR)

To analyze the mRNA expression levels of inflammation markers associated with IVDD, quantitative reverse transcription polymerase chain reaction (qRT-PCR) was employed. Total RNA was extracted from NPCs using TRIzol Reagent™ (AG21102; Accurate Biology, Hunan, China) according to the manufacturer’s instructions. The concentration and purity of the extracted RNA were assessed using a NanoDrop spectrophotometer (Thermo Fisher Scientific), with A260/A280 ratios between 1.8 and 2.0 indicating high purity. cDNA synthesis was performed in a 20 μL reaction volume following established protocols. For qRT-PCR, the SYBR Premix Ex TaqTM II kit (MR05101M; Monad, Suzhou, China) was utilized in a CFX Connect Touch Real-Time System (1855204; Bio-Rad, China). The reaction mixture consisted of cDNA template, gene-specific primers designed using Primer3 software, and SYBR Green Master Mix. Each sample underwent qRT-PCR in triplicate. Thermal cycling conditions were as follows: an initial denaturation step at 95°C for 30 seconds, followed by 40 cycles of denaturation at 95°C for 5 seconds, and annealing/extension at 60°C for 30 seconds. Melting curve analysis was performed to verify the specificity of the qRT-PCR products. Relative gene expression levels were calculated using the 2^−ΔΔCt method, normalizing against glyceraldehyde-3-phosphate dehydrogenase as the internal reference gene. Fold changes in mRNA expression were determined by comparing the normalized target gene expression levels in experimental groups to those in control groups.

Western Blot Analysis

Following two washes with PBS, cells were lysed using a mixture of IP lysis buffer (P0013; Beyotime), proteasome inhibitor (P1045; Beyotime), and phosphatase inhibitor (1081; Beyotime). The lysates were then centrifuged at 12,000 rpm for 10 minutes to remove cellular debris, and the supernatant was collected for protein sample preparation. Equal amounts of protein samples were subsequently loaded onto a 10% Tris-glycine sodium dodecyl sulfate-polyacrylamide gel for electrophoretic separation. Upon completion of the electrophoresis, the separated proteins were transferred to a polyvinylidene fluoride (PVDF) membrane. To minimize non-specific binding, the membrane was blocked using a buffer containing 5% skim milk for 2 hours. Following this, the PVDF membrane was incubated overnight at 4°C with the following primary antibodies: anti-actin (23660-1-AP; Proteintech), anti-ARTN (ab178434; Abcam), anti-IL-1β (AF5103; Affinity), and anti-MMP13 (AF5355; Affinity). The membrane was then incubated with a secondary antibody (A0208; Beyotime) at 37°C for 2 hours. Finally, the membrane was analyzed using an enhanced chemiluminescence kit (P0018AS; Beyotime). Quantitative analysis was performed using ImageJ software (Media Cybernetics, MD, USA).

Statistical Analysis

Experiments were performed in triplicate to ensure accuracy. Data analysis was conducted using GraphPad Prism 9.5, with results presented as mean ± standard deviation (SD). Statistical significance was evaluated using the Student’s t-test, one-way ANOVA followed by post hoc least significant difference test, and two-way ANOVA with multiple t-tests. A significance level of P < 0.05 was considered statistically significant.

Results

Genes Associated With IVDD in Blood

We processed 19,942 whole blood gene data from the GWAS database and obtained 5,430 exposure genes after removing linkage disequilibrium and weak instrumental variables. We performed univariate two-sample Mendelian randomization analysis on these genes and IVDD, and excluded exposures with pleiotropy and heterogeneity, while retaining exposure genes with consistent beta value directions. Ultimately, we identified 240 exposure factors strongly associated with IVDD (Supplementary Files 2), among which 115 were inductive exposures and 125 were protective exposures. To further understand the genes related to IVDD at the transcriptome level, we conducted differential expression analysis (logFC>0.58, p<0.05) on datasets GSE15227 and GSE23130 using the limma R package. The results showed that a total of 765 genes were associated with the induction of IVDD, while 217 genes exhibited protective effects (Figure 1A and B). Through Venn diagram intersection analysis, we discovered 7 genes significantly associated with IVDD (Figure 1C and D). Among them, the expression levels of AEN (OR:1.086, 95% CI: 1.014 ~ 1.165, P = 0.019), CLEC11A (OR:1.072, 95% CI: 1.016 ~ 1.132, P = 0.011), HMGN1 (OR:1.086, 95% CI: 1.030 ~ 1.144, P = 0.002), TAF7 (OR:1.060, 95% CI: 1.003 ~ 1.120, P = 0.037), and TREM1 (OR:1.037, 95% CI: 1.003 ~ 1.072, P = 0.033) were significantly associated with an increased risk of IVDD. In contrast, the expression levels of DNAJA4 (OR:0.926, 95% CI: 0.876 ~ 0.978, P = 0.006) and LRRC25 (OR:0.958, 95% CI: 0.924 ~ 0.993, P = 0.020) were significantly associated with a decreased risk of IVDD (Figure 2A and B).

Figure 1 Transcriptome sequencing analysis of IVDD. (A) Heatmap of differential analysis. (B) Volcano plot of differential analysis, logFC > 0.58, P < 0.05. (C) Venn diagram showing the intersection of upregulated differentially expressed genes and risk genes with OR > 1 in the MR results. (D) Venn diagram showing the intersection of downregulated differentially expressed genes and protective genes with OR < 1 in the MR results.

Figure 2 MR analysis. (A) Causal relationships between key genes and IVDD. nsnp: total number of instrumental variables used for analysis. OR > 1: exposure increases the risk of the outcome. OR < 1: exposure reduces the risk of the outcome. (B) Reverse MR analysis of the effect of key genes on the risk of IVDD. (C) Scatter plots of the MR analysis of key genes and the risk of IVDD. Bold numbers in the figures represent statistically significant results (p < 0.05).

To further understand whether these important genes have a reverse causal relationship with IVDD, we conducted reverse MR analysis. The results revealed that DNAJA4 (OR:0.753, 95% CI: 0.584 ~ 0.971, P = 0.029) had a reverse causal relationship with IVDD (Figure 2C), suggesting that IVDD may affect the expression of DNAJA4. Therefore, in subsequent mediation analysis, we excluded DNAJA4 to avoid potential interference from confounding factors. The remaining 6 genes were considered potentially important genes inducing IVDD.

Single-Cell Analysis Reveals IVDD

We analyzed the single-cell dataset GSE244889, which includes one mildly degenerated intervertebral disc sample and one severely degenerated intervertebral disc sample. After strict quality control, we performed UMAP for visualization. Subsequently, we divided the cells into 8 clusters (Figure 3A) and annotated them as identifiable cell types using the SingleR R package (Figure 3B). The main cell types included Chondrocytes, Monocytes, Endothelial_cells, Erythroblasts, and B_cells. To further determine the important role of Chondrocytes in degeneration,²⁷ we divided Chondrocytes into four clusters: Chondrocytes 1, Chondrocytes 2, Chondrocytes 3, and Chondrocytes 4. We then analyzed the expression of the 6 genes associated with the risk of IVDD in each cluster. We found that AEN and LRRC25 had low expression in all clusters, HMGN1 was mainly highly expressed in Monocytes, CLEC11A was mainly highly expressed in Chondrocytes 1, Chondrocytes 2, and Chondrocytes 3, TREM1 was mainly highly expressed in Chondrocytes 1, and TAF7 was mainly highly expressed in Endothelial_cells (Figure 3C and D). Subsequently, we performed differential analysis between the two samples in each cluster (Figure 4A and B). The results showed that in Chondrocytes 1, CLEC11A was significantly highly expressed in severe IVDD compared to mild IVDD (P=0.00057); in Chondrocytes 2, CLEC11A was significantly highly expressed in severe IVDD compared to mild IVDD (P=7.94E-22); in Chondrocytes 3, TREM1 was significantly highly expressed in severe IVDD compared to mild IVDD (P=3.92E-10); in Monocytes, HMGN1 was significantly highly expressed in severe IVDD compared to mild IVDD (P=0.001). These results preliminarily suggest that CLEC11A, TREM1, and HMGN1 may be high-risk genes for the occurrence of IVDD.

Figure 3 Clustering analysis of single-cell data and expression of key genes in IVDD. (A) UMAP plot segregating cells into 8 clusters. (B) SingleR annotating cells into 5 cell types. (C and D) Expression levels of 6 key genes across different clusters.

Figure 4 Differential analysis between clusters and pseudotemporal analysis using Monocle 3 in single-cell data. (A and B) In clusters 0 and 1, CLEC11A is significantly upregulated in the severe IVDD group compared to the mild IVDD group. In cluster 2, TREM1 is significantly upregulated in the severe IVDD group compared to the mild IVDD group. In cluster 3, HMGN1 is significantly upregulated in the severe IVDD group compared to the mild IVDD group. (C and D) Pseudotemporal analysis showing the progression of Chondrocytes 1, Chondrocytes 2, Chondrocytes 3, and Chondrocytes 4 during the IVDD process. (E) Sample distribution of mild and severe IVDD on the UMAP plot.

Since Chondrocytes are the key cell type causing IVDD, we performed pseudotime analysis on Chondrocytes to determine the degree of IVDD by inferring the developmental trajectory of cells. The results showed that Chondrocytes 2 was located at the initial stage of degeneration, Chondrocytes 1 was located at the middle stage of degeneration, and Chondrocytes 3 and Chondrocytes 4 were located at the terminal stage of degeneration (Figure 4C and D). The sample distribution map (Figure 4E) results showed that the main intersection region of the mildly degenerated intervertebral disc sample and the severely degenerated intervertebral disc sample was located in Chondrocytes 2, indicating that Chondrocytes 2 might be an important cluster causing the occurrence of IVDD, which is consistent with the pseudotime analysis results.

Mediating Role of Inflammatory Factors in the Relationship Between TREM1/CLEC11A/HMGN1 and IVDD

We first conducted a comprehensive MR study to reveal the causal relationship between 91 inflammatory factors and IVDD. After removing linkage disequilibrium and weak instrumental variables, we performed univariate two-sample Mendelian randomization analysis and excluded inflammatory factors with pleiotropy in the results. Finally, we determined that two inflammatory proteins, ARTN (OR:1.069, 95% CI: 1.001 ~ 1.141, P = 0.045) and CCL8 (OR:1.077, 95% CI: 1.027 ~ 1.128, P = 0.002), had significant causal relationships with the occurrence of IVDD (Figure 5A–C), and the inflammatory protein OPG (OR:0.907, 95% CI: 0.837 ~ 0.984, P = 0.019) had a significant causal relationship with protection against IVDD (Figure 5A and D). Subsequently, we used TREM1/CLEC11A/HMGN1 as exposures and the three inflammatory factors as outcomes, respectively, and performed univariate two-sample Mendelian randomization analysis. The results showed that only CLEC11A had a significant causal relationship with the inflammatory factor ARTN (OR:1.078, 95% CI: 1.004 ~ 1.158, P = 0.038, Beta1=0.075) (Figure 6A–D). Subsequent heterogeneity tests (P=0.580) and pleiotropy tests (P=0.587) indicated no heterogeneity or pleiotropy. We then conducted MR analysis on the inflammatory protein ARTN and IVDD to calculate Beta2. To satisfy a key assumption of mediation MR analysis, ie, the genetic instrumental variables used to estimate beta1 and beta2 should be independent, we first excluded the significant SNPs used to calculate Beta1. The results showed that the inflammatory factor ARTN had a significant causal relationship with IVDD (OR:1.069, 95% CI: 1.001 ~ 1.141, P = 0.045, Beta2=0.066), and subsequent heterogeneity tests (P= 0.257) and pleiotropy tests (P= 0.960) indicated no heterogeneity or pleiotropy. To calculate the total effect value beta all, we performed MR analysis on CLEC11A and IVDD and conducted heterogeneity and pleiotropy tests. The results showed that CLEC11A had a significant causal relationship with IVDD (OR:1.072, 95% CI: 1.015 ~ 1.131, P = 0.011, Beta all=0.069), and heterogeneity tests (P=0.523) and pleiotropy tests (P=0.462) indicated no heterogeneity or pleiotropy. The mediation analysis results showed a mediated effect value of 0.00506 (0.00138, 0.00875) and a mediated proportion of 7.26% (1.98%, 12.5%) (Table 1). The results suggest that CLEC11A may promote the risk of IVDD by promoting the expression of the inflammatory protein ARTN.

Table 1 Mediation Effects and Proportion Mediated by Inflammatory Proteins

Figure 5 MR analysis of 91 inflammatory proteins and the risk of IVDD. (A) Forest plots of three positive inflammatory proteins: ARTN, CCL8, and OPG. (B–D) Scatter plots of three positive inflammatory proteins: ARTN, CCL8, and OPG.

Figure 6 MR analysis of inflammatory proteins as mediators. (A) Scatter plot of the MR analysis between CLEC11A and the mediating inflammatory protein ARTN. (B) Forest plot of the MR analysis between CLEC11A and the mediating inflammatory protein ARTN. (C) Scatter plot of the MR analysis between the adjusted mediating inflammatory protein ARTN and the risk of IVDD. (D) Forest plot of the MR analysis between the adjusted mediating inflammatory protein ARTN and the risk of IVDD.

Mediating Role of Metabolites in the Relationship Between TREM1/CLEC11A/HMGN1 and IVDD

Through a comprehensive MR study, we revealed the causal relationship between 1,400 serum metabolites and IVDD. After removing linkage disequilibrium and weak instrumental variables, we performed univariate two-sample Mendelian randomization analysis and excluded metabolites with pleiotropy in the results. Finally, we determined that 64 serum metabolites had a causal relationship with IVDD (Supplementary Files 3). Subsequently, we used TREM1/CLEC11A/HMGN1 as exposures and the 64 serum metabolites as outcomes, respectively, and performed univariate two-sample Mendelian randomization analysis. The results showed that CLEC11A had significant causal relationships with three serum metabolites: X-12731 levels (OR:0.906, 95% CI: 0.852 ~ 0.996, P = 0.043, Beta1=−0.097), X-18901 levels (OR:1.090, 95% CI: 1.007 ~ 1.179, P = 0.034, Beta1=0.085), and 1-palmitoyl-2-arachidonoyl-gpc (16:0/20:4n6) levels (OR:1.097, 95% CI: 1.013 ~ 1.188, P = 0.023, Beta1=0.092). TREM1 had significant causal relationships with three serum metabolites: 5alpha-androstan-3beta,17alpha-diol disulfate levels (OR:0.948, 95% CI: 0.899 ~ 0.999, P = 0.048, Beta1=−0.053), 1-(1-enyl-stearoyl)-GPE (p-18:0) levels (OR:1.060, 95% CI: 1.012 ~ 1.111, P = 0.014, Beta1=0.058), and Hypotaurine levels (OR:1.046, 95% CI: 1.000 ~ 1.093, P = 0.049, Beta1=0.044). HMGN1 had significant causal relationships with two serum metabolites: 1-stearoyl-2-arachidonoyl-gpc (18:0/20:4) levels (OR:0.921, 95% CI: 0.853 ~ 0.994, P = 0.034, Beta1=−0.082) and Glucose to mannitol to sorbitol ratio (OR:1.091, 95% CI: 1.011 ~ 1.178, P = 0.024, Beta1=0.087) (Figure 7A and B). Subsequent heterogeneity tests and pleiotropy tests indicated no heterogeneity or pleiotropy. We then performed MR analysis on the 8 serum metabolites and IVDD to calculate Beta2, first excluding the significant SNPs used to calculate Beta1. The results showed that the 8 serum metabolites had significant causal relationships with IVDD: X-12731 levels (OR:0.942, 95% CI: 0.889 ~ 0.999, P = 0.046, Beta2= −0.059), X-18901 levels (OR:1.078, 95% CI: 1.021 ~ 1.138, P = 0.007, Beta2= 0.075), 1-palmitoyl-2-arachidonoyl-gpc (16:0/20:4n6) levels (OR:0.970, 95% CI: 0.942 ~ 1.000, P = 0.048, Beta2=−0.030), 5alpha-androstan-3beta,17alpha-diol disulfate levels (OR:1.072, 95% CI: 1.020 ~ 1.126, P = 0.006, Beta2= 0.069), 1-(1-enyl-stearoyl)-GPE (p-18:0) levels (OR:0.926, 95% CI: 0.869 ~ 0.987, P = 0.018, Beta2=−0.076), Hypotaurine levels (OR:0.936, 95% CI: 0.887 ~ 0.989, P = 0.018, Beta2=−0.065), 1-stearoyl-2-arachidonoyl-gpc (18:0/20:4) levels (OR:0.976, 95% CI: 0.953 ~ 0.999, P = 0.045, Beta2=−0.024), and Glucose to mannitol to sorbitol ratio (OR:0.905, 95% CI: 0.832 ~ 0.984, P = 0.019, Beta2=−0.100). Subsequent heterogeneity tests and pleiotropy tests indicated no heterogeneity or pleiotropy (Figure 8A and B). To calculate the total effect value beta all, we performed MR analysis on TREM1/CLEC11A/HMGN1 and IVDD and conducted heterogeneity and pleiotropy tests. The results showed that the 3 genes had significant causal relationships with IVDD: CLEC11A (OR:1.072, 95% CI: 1.015 ~ 1.131, P = 0.011, Beta all=0.069), TREM1 (OR: 1.037, 95% CI: 1.003 ~ 1.072, P = 0.033, Beta all=0.036), and HMGN1 (OR: 1.085, 95% CI: 1.030 ~ 1.143, P = 0.002, Beta all=0.082) (Figure 9A and B), (Figure 10A–C). Heterogeneity tests and pleiotropy tests indicated no heterogeneity or pleiotropy. The mediation analysis results showed that CLEC11A might promote the risk of IVDD by promoting the expression of serum metabolites X-12731 levels (Mediated effect: 0.00581 (0.000297, 0.0113), Mediated proportion: 8.33% (0.426%, 16.2%)) and X-18901 levels (Mediated effect: 0.00645 (0.00184, 0.0111), Mediated proportion: 9.26% (2.64%, 15.9%)). The mediating effects of the remaining 6 serum metabolites did not meet the requirements of mediation Mendelian analysis and were excluded (Table 2).

Table 2 Mediation Effects and Proportion Mediated by Serum Metabolites

Figure 7 MR analysis of metabolites as mediators. (A) MR analysis of 3 key genes and 8 positive metabolites. nsnp: total number of instrumental variables used for analysis. OR > 1: exposure increases the risk of the outcome. OR < 1: exposure reduces the risk of the outcome. P < 0.05 is considered statistically significant. (B) Scatter plots of the MR analysis between 3 key genes and 8 positive metabolites. Bold numbers in the figures represent statistically significant results (p < 0.05).

Figure 8 MR analysis of metabolites as mediators. (A) Circle plot of the MR analysis between 8 metabolites and the risk of IVDD after adjustment. (B) Scatter plots of the MR analysis between 8 metabolites and the risk of IVDD after adjustment.

Figure 9 MR analysis of metabolites as mediators. (A) MR analysis of three key genes and the risk of IVDD. nsnp: total number of instrumental variables used for analysis. OR > 1: exposure increases the risk of the outcome. OR < 1: exposure reduces the risk of the outcome. beta all: total effect value of the mediation analysis. P < 0.05 is considered statistically significant. (B) Scatter plots of the MR analysis between three key genes and the risk of IVDD. Bold numbers in the figures represent statistically significant results (p < 0.05).

Figure 10 Mediation MR analysis. (A) Results of the mediation analysis with the inflammatory protein ARTN as a mediator of the causal relationship between CLEC11A and the risk of IVDD. (B) Results of the mediation analysis with the unidentified serum metabolite X-12731 as a mediator of the causal relationship between CLEC11A and the risk of IVDD. (C) Results of the mediation analysis with the unidentified serum metabolite X-18901 as a mediator of the causal relationship between CLEC11A and the risk of IVDD. Bold numbers in the figures represent statistically significant results (p < 0.05).

CLEC11A Mediates the Expression of ARTN, Promoting the Progression of IVDD

To elucidate the molecular mechanisms by which CLEC11A facilitates the progression of IVDD, we manipulated the expression levels of the CLEC11A gene in NPCs via overexpression and knockdown techniques. Western blot assays revealed that overexpression of CLEC11A resulted in a substantial increase in the protein levels of ARTN, IL1-β, and MMP13. Similarly, siRNA-mediated knockdown of CLEC11A significantly decreased the protein levels of these factors. Quantitative analyses of the protein levels confirmed the statistical significance of these results (P < 0.05) (Figure 11A). qPCR assays further supported these findings; overexpression of CLEC11A significantly upregulated the mRNA levels of ARTN, IL1-β, and MMP13 compared to the control group. Conversely, knockdown of CLEC11A using siRNA led to a notable reduction in the mRNA levels of these markers (P < 0.05) (Figure 11B). To further explore the role of ARTN in the progression of IVDD, we conducted additional genetic manipulation experiments on the ARTN gene in NPCs. Western blot analysis showed that ARTN overexpression significantly enhanced the protein levels of IL1-β and MMP13. Conversely, knockdown of ARTN resulted in a significant reduction in these protein levels, with quantitative analysis confirming the statistical significance of these changes (P < 0.05) (Figure 11C). qPCR assays further supported that the overexpression of ARTN markedly increased the mRNA levels of IL1-β and MMP13, while its knockdown significantly reduced their expression (P < 0.05) (Figure 11D).

Figure 11 Modulation of ARTN, IL1-β, and MMP13 by CLEC11A and ARTN in Nucleus Pulposus Cells (A) Western blot analyses show significant changes in protein levels of ARTN, IL1-β, and MMP13 following CLEC11A overexpression and knockdown. (B) qPCR results indicate that CLEC11A manipulation significantly alters mRNA levels of ARTN, IL1-β, and MMP13 compared to control. (C) Protein levels of IL1-β and MMP13 significantly fluctuate with ARTN genetic modifications, as confirmed by Western blot. (D) Corresponding qPCR data demonstrate significant ARTN-mediated changes in mRNA levels of IL1-β and MMP13. Mean ± standard deviation (SD). *p < 0.05, **p < 0.01, ***p < 0.001, ****p < 0.0001.

Discussion

In this study, we aimed to identify potential risk genes that may induce IVDD and explore the inflammatory proteins and metabolites mediating the causal effects of these genes on IVDD. To achieve this goal, we conducted a series of analyses, including transcriptome differential expression analysis, single-cell sequencing analysis, traditional and mediational Mendelian randomization, and in vitro validation. Our results showed that CLEC11A might increase the risk of IVDD by upregulating the expression of the inflammatory protein ARTN and two unidentified serum metabolites, X-12731 and X-18901.

By integrating MR analysis and transcriptomics analysis, we identified a set of genes significantly associated with IVDD. We found that increased expression of genes such as AEN, CLEC11A, HMGN1, TAF7, and TREM1 can significantly increase the risk of IVDD, while increased expression of LRRC25 has a protective effect in inhibiting IVDD. These findings reveal some potential molecular events and signaling pathways that may exist during the occurrence and development of IVDD.

Decreased number and vitality of intervertebral disc cells are considered important factors in initiating IVDD. AEN and TAF7 may play crucial roles in the early stages of IVDD by promoting apoptosis of intervertebral disc cells and regulating cell proliferation and differentiation. Studies have shown that the tumor suppressor gene p53 produces transcription factors that regulate cell metabolism and survival, playing a key role in intervertebral disc cell apoptosis, metabolism, and maintenance of microenvironmental homeostasis. AEN is a p53-positively regulated apoptosis-promoting protein that participates in the final stage of cell apoptosis by degrading chromatin DNA.²⁸ TAF7 is a key component of the RNA polymerase II transcription initiation complex TFIID, playing an important role in transcription initiation and elongation. Research has shown that TAF7, as an essential factor for transcription, is involved in cell cycle regulation;²⁹ simultaneously, TAF7 also acts as a transcription elongation factor, mediating the transition from initiation to elongation in multiple rounds of transcription.³⁰

Persistent activation of inflammatory responses is considered a major cause of accelerated IVDD. CLEC11A, HMGN1, and TREM1 play key roles in the progression stage of IVDD, mainly by mediating inflammatory responses in intervertebral disc tissue. Studies have shown that C-type lectins play important roles in innate immune responses, including mediating pathogen neutralization, leukocyte transport, and inflammatory reactions, by binding to polysaccharide residues on the surface of pathogens. However, dysregulation of C-type lectin signaling pathways can lead to inflammatory diseases such as arthritis.³¹ CLEC11A encodes a C-type lectin-like protein, suggesting that it may promote IVDD by regulating local inflammatory responses in the intervertebral disc. HMGN1, as an endogenous danger signal molecule, can induce inflammatory responses through TLR4-dependent pathways in its extracellular form. Research has found that HMGN1 is highly expressed in diabetic kidney injury, and its silencing can alleviate tubular cell damage and renal inflammation;³² furthermore, it can induce innate immune tolerance in peripheral monocytes through the TLR4 pathway.³³ TREM1, on the other hand, belongs to an inflammatory receptor expressed by myeloid cells and can amplify inflammatory responses mediated by neutrophils and monocytes by enhancing the release of pro-inflammatory cytokines and chemokines. Previous studies have shown that silencing TREM1 can significantly alleviate inflammatory responses and tissue damage in IVDD models.^34,35

In contrast, LRRC25 may play a protective role in the progression stage of IVDD by inhibiting inflammatory responses in intervertebral disc tissue. Research has found that LRRC25 can act as a negative regulator of the NF-κB inflammatory signaling pathway, inhibiting NF-κB activation and inflammatory responses by ubiquitinating and degrading the p65/RelA subunit.³⁶ In addition, LRRC25 can also suppress type I interferon antiviral inflammatory signaling by degrading the RIG-I receptor through the autophagy pathway.³⁷ LRRC25 deficiency can accelerate the progression of cardiac hypertrophy and dysfunction by promoting TGF-β1-induced fibrosis and inflammatory responses.³⁸ These findings suggest that LRRC25 may be an endogenous inflammatory suppressor that plays a general role in inhibiting chronic inflammation in tissues. Therefore, enhancing LRRC25 function may be a potential therapeutic strategy for slowing down and curbing the progression of IVDD.

Subsequently, we analyzed the molecular mechanisms of IVDD from the perspective of cellular heterogeneity using single-cell transcriptome sequencing technology. Our results revealed that CLEC11A, TREM1, and HMGN1 were significantly upregulated in the Chondrocytes and Monocytes of degenerated intervertebral disc tissues, suggesting that they may play key roles in the progression stage of IVDD by affecting chondrocyte function and inflammatory responses. These findings are consistent with our previous discoveries obtained through MR and transcriptomic analyses.

Our analysis first divided the cells into 8 clusters and identified the main cell types, including Chondrocytes, Monocytes, Endothelial_cells, Erythroblasts, and B_cells. Among them, Chondrocytes were further divided into 4 subclusters (Chondrocytes 1–4). Previous studies have shown that Chondrocytes play a key role in maintaining the structure and function of the intervertebral disc.³⁹ Chondrocyte apoptosis and matrix degradation are considered important pathological changes in IVDD.^40,41 In this study, we found that CLEC11A was highly expressed in Chondrocytes 1 and Chondrocytes 2 and significantly upregulated in severely degenerated samples. CLEC11A encodes a C-type lectin-like receptor involved in regulating cell proliferation and differentiation. Therefore, CLEC11A may participate in the pathological process of IVDD by affecting the function of Chondrocytes. TREM1 is another differentially expressed gene in Chondrocytes, mainly highly expressed in Chondrocytes 3 and significantly upregulated in the severely degenerated group. TREM1 is an immune-related receptor that plays an important role in inflammatory responses. Previous studies suggest that inflammation plays a crucial role in IVDD,³⁵ and abnormal expression of TREM1 may accelerate the progression of IVDD by promoting inflammatory responses. In addition, we also found that HMGN1 was significantly upregulated in Monocytes. HMGN1 is a non-histone chromatin-binding protein that can affect gene transcription. Monocyte infiltration is one of the pathological manifestations of IVDD,⁴² and the high expression of HMGN1 in Monocytes suggests that it may be involved in monocyte infiltration and inflammatory responses.

Our pseudotemporal analysis further revealed the dynamic change trajectory of chondrocytes during the IVDD process. The results showed that Chondrocytes 2 might be in the initial stage of degeneration, while Chondrocytes 1, 3, and 4 corresponded to the intermediate and terminal stages, respectively. This suggests that chondrocytes may undergo a continuous change in functional state, transforming from a normal steady state to a degenerated phenotype. The significant intersection of mildly and severely degenerated samples in the Chondrocytes 2 cluster further verified that this cluster might be the initial cell population of IVDD. CLEC11A was highly expressed in the Chondrocytes 2 cluster and significantly upregulated in severely degenerated samples, further indicating that CLEC11A might be a key gene triggering the early stage of IVDD. As the degeneration process progresses, the persistent activation of inflammatory responses may further accelerate chondrocyte damage and matrix degradation, leading to progressive loss of intervertebral disc structure and function. This finding is consistent with previous studies, which suggest that the decrease in the number and vitality of intervertebral disc cells is an important factor in initiating the degeneration process.⁴³

Subsequently, we used mediation MR analysis to further explore the mediating effects of inflammatory factors and metabolites on the relationship between TREM1/CLEC11A/HMGN1 genes and IVDD from a genetic perspective. We found that CLEC11A may increase the risk of IVDD by upregulating the levels of the inflammatory factor ARTN and serum metabolites X-12731 and X-18901. Inflammatory response is one of the core pathological processes that trigger and promote IVDD. Numerous studies have shown that pro-inflammatory cytokines such as TNF-α and IL-1β are significantly elevated in degenerated intervertebral disc tissues, inducing the expression of matrix metalloproteinases, degrading collagen, proteoglycans, and other extracellular matrix components, disrupting the structure and function of intervertebral discs, and accelerating the degeneration process.⁴⁴ Therefore, inhibiting inflammatory responses is considered a key strategy for slowing down and reversing IVDD. ARTN, as an important member of the GDNF family, plays a crucial regulatory role in various inflammatory diseases. A study found that ARTN is highly expressed in degenerated intervertebral disc tissues, especially in severely degenerated tissues, suggesting that ARTN may be an important promoting factor in the progression of IVDD.⁴⁵ This study confirmed the causal relationship between CLEC11A and an increased risk of IVDD by upregulating ARTN expression from a genetic epidemiological perspective. Consistent with these findings, in vitro cell experiments demonstrated that overexpression of CLEC11A enhances ARTN protein expression, thereby promoting inflammatory responses in IVDD. These results further support the key role of CLEC11A in the pathogenesis of IVDD. Metabolic disorders are another important pathological basis for IVDD. The nutrient supply of intervertebral discs mainly relies on diffusion, and metabolic waste is also removed through diffusion; therefore, the homeostasis of the metabolic microenvironment is crucial for the survival and function of intervertebral disc cells. Studies have found that metabolic products such as lactic acid accumulate in large amounts in degenerated intervertebral discs, causing a decrease in matrix pH, accelerating matrix degradation, and inducing multiple pathological reactions such as inflammation, apoptosis, and autophagy, ultimately leading to the loss of intervertebral disc structure and function.⁴⁶ Therefore, restoring metabolic balance and clearing metabolic waste are novel therapeutic strategies for slowing down the progression of IVDD. For example, urolithin A, as a gut microbiota metabolite, has various pharmacological activities such as anti-inflammatory and antioxidant effects, and can improve the metabolic microenvironment of intervertebral disc cells, reduce stress-induced extracellular matrix degradation and nucleus pulposus cell apoptosis, thereby slowing down IVDD.⁴⁷ The results of this study showed that CLEC11A can mediate the risk of IVDD by upregulating serum metabolites X-12731 and X-18901, which not only emphasizes the important role of metabolic factors in the pathogenesis of IVDD but also provides new potential targets for metabolic regulation and regenerative repair of IVDD. Although the chemical nature and biological functions of X-12731 and X-18901 are currently unclear, subsequent studies can use metabolomics, cheminformatics, and other means to deeply analyze their structural characteristics, construct their action patterns in the CLEC11A signaling network, explore their effects on the metabolism, inflammation, matrix synthesis, and other phenotypes of intervertebral disc cells. Furthermore, targeting CLEC11A may present a novel therapeutic strategy to attenuate inflammatory processes and decelerate the progression of IVDD, thereby enhancing patient outcomes and alleviating the associated economic burden.

This study employs a multi-omics approach to elucidate the potential molecular mechanisms underlying IVDD at the genomic, cellular, and metabolic levels; however, it is not without limitations. Firstly, selection bias and the limited strength of genetic instrumental variables may impact the analytical capacity of our findings. Secondly, we did not apply multiple testing corrections, such as Bonferroni adjustment, to the results of our MR analysis. This decision was made because our primary aim was to identify potential biomarkers or therapeutic targets for IVDD, and a strict Bonferroni criterion could potentially exclude meaningful indicators. Furthermore, while the mediation MR analysis suggests that the inflammatory protein ARTN and two unidentified serum metabolites, X-12731 and X-18901, mediate the causal relationship between CLEC11A and IVDD, our understanding of their biological functions remains limited at this stage.

Conclusion

In conclusion, this study utilized a comprehensive whole-genome multi-omics integration strategy to systematically elucidate the molecular pathogenesis of IVDD from genetic, transcriptomic, and experimental perspectives. We identified pivotal genes such as CLEC11A, TREM1, and HMGN1 that likely play crucial roles in IVDD by influencing chondrocyte function and inflammatory responses. Notably, CLEC11A was pinpointed as a central pathogenic gene in IVDD. Subsequent mediational Mendelian randomization analysis further revealed that CLEC11A could significantly elevate IVDD risk through the upregulation of inflammatory cytokine ARTN and the serum metabolites X-12731 and X-18901, underlining the significant roles of inflammation and metabolic dysregulation in the etiology of IVDD. Confirmatory in vitro cellular experiments demonstrated that CLEC11A exacerbates IVDD progression by enhancing ARTN expression and amplifying inflammatory responses. These insights not only provide a fresh perspective on the complex etiological mechanisms of IVDD but also lay an essential foundation for precision medicine and the optimization of intervention strategies for this disorder.