Abstract Background and objectives Intervertebral disc degeneration (IVDD) is a leading cause of chronic low back pain, characterized by extracellular matrix (ECM) degradation, excessive inflammation activation, and increased cell apoptosis. LGR6, a receptor known for its role in tissue regeneration, has recently been implicated in modulating macrophage efferocytosis, a process critical for clearing apoptotic cells and maintaining tissue homeostasis. This study aimed to investigate the role of LGR6 in regulating IVDD progression and to focus on its impact on macrophage efferocytosis, ECM regulation, and apoptosis in nucleus pulposus cells (NPCs). Methods A comprehensive bioinformatic analysis was performed using datasets [36]GSE56081 and [37]GSE70362 to identify differentially expressed genes (DEGs) and gene modules associated with IVDD. Principal component analysis (PCA), volcano plots, and hierarchical clustering were utilized to assess gene expression patterns. Weighted Gene Co-Expression Network Analysis (WGCNA) was employed to identify gene modules correlated with IVDD, and integrative analysis pinpointed key genes and pathways. In vitro, LGR6 expression in macrophages was manipulated through shRNA interference and overexpression assay. The effects of LGR6 on macrophage efferocytosis, ECM synthesis, and apoptosis were assessed. An in vivo IVDD model was established in mice via disc puncture to evaluate the impact of LGR6 modulation on disc degeneration. Results Bioinformatic analysis revealed distinct gene expression profiles between control and IVDD samples, with key gene modules identified by WGCNA showing strong correlations with IVDD. Integrative analysis highlighted critical pathways, including ECM-receptor interaction and efferocytosis, that are potentially regulated by several key genes including SERPINA1, THBS4, ELMO1, LGR6, and ITGB8. Of those genes, LGR6 appeared to be the gene most closely related to IVDD severity. In addition, the mRNA level and protein level of LGR6 in macrophages co-cultured with IL-1β-treated NPCs were raised significantly, compared to the control group. In vitro, LGR6 overexpression enhanced macrophage efferocytosis. Meanwhile, under co-culturing with IL-1β-treated NPCs, LGR6 overexpression in macrophages led to increased expression of ECM components such as COL2A1 and decreased expression of matrix-degrading enzymes like MMP13, indicating a protective effect against matrix degradation. Additionally, LGR6 overexpression inhibited IL-1β-induced apoptosis in NPCs by upregulating anti-apoptotic proteins (BCL2) and downregulating pro-apoptotic markers (cleaved caspase 3 and BAX). Conversely, LGR6 knockdown impaired macrophage efferocytosis and exacerbated NPCs apoptosis. In the mouse IVDD model, promoting efferocytosis resulted in improved ECM integrity and reduced apoptosis; and suppressing efferocytosis caused opposite effect, further supporting the protective role of LGR6-related efferocytosis in IVDD. Conclusions LGR6 significantly contributes to the protective effects on IVDD by modulating macrophage efferocytosis, enhancing ECM synthesis, and reducing apoptosis in NPCs. These findings highlight that LGR6 could be a promising therapeutic target for IVDD, with its dual role in regulating immune responses and preserving tissue integrity. Future studies are necessary to evaluate the clinical potential of LGR6-based therapies in treating degenerative disc diseases. Supplementary Information The online version contains supplementary material available at 10.1186/s12967-025-06427-0. Keywords: Intervertebral disc degeneration, LGR6, Macrophages, Efferocytosis, Apoptosis, Nucleus pulposus Introduction Intervertebral disc degeneration (IVDD) is a prevalent dysfunction affecting millions of people worldwide, imposing enormous socioeconomic burden worldwide [[38]1]. Previous studies indicated that up to 80% of people will experience low back pain during their lives, with IVDD being the major contributing factor [[39]1, [40]2]. Characterized by pathological changes such as disc dehydration, loss of disc height, and structural disorganization [[41]3], considerable treatment challenges existed in treating IVDD. The complexity of IVDD involves genetic, biochemical, and biomechanical factors that contribute to its development and progression, making it a multifactorial disease that requires comprehensive research approaches [[42]4]. Despite extensive research, the precise mechanisms underlying its onset and progression remain elusive, necessitating further investigation at the molecular biology level. Inflammatory cell infiltration is a typical feature of IVDD, with macrophages playing a pivotal role during the degenerative process, such as waste scavenging and inflammation regulation [[43]5]. Efferocytosis, a process by which macrophages clear apoptotic cells, is essential for maintaining the tissue homeostasis and limiting chronic inflammation [[44]6]. Via efferocytosis, macrophages recognize and engulf apoptotic cells through specific receptors, such as MerTK and phosphatidylserine receptors, to remove the potentially harmful cellular debris [[45]7]. This process not only prevents secondary cellular apoptosis and the release of pro-inflammatory cellular contents but also promotes the resolution of inflammation by stimulating anti-inflammatory pathways [[46]8]. In the context of degenerative diseases, efficient efferocytosis is crucial for mitigating chronic inflammation and facilitating tissue repair [[47]9]. Impaired efferocytosis has been implicated in the progression of various chronic inflammatory conditions, including multiple sclerosis, Alzheimer’s disease, glioblastoma, and Rett syndrome [[48]9]. However, the role of efferocytosis in IVDD has not yet be well investigated [[49]10]. Investigating the role of efferocytosis and its critical regulators in IVDD could facilitate the understanding of the disease and the identification of new therapeutic targets. By exploring the pathways and genes involved in macrophage-mediated efferocytosis and their impact on disc degeneration, this research aims to identify novel therapeutic targets that could be leveraged to develop more effective treatments for IVDD. Understanding these mechanisms is crucial for devising strategies to mitigate the impact of IVDD on patients’ quality of life and reduce its substantial socioeconomic burden. Methods and materials Gene expression analysis The study utilized datasets [50]GSE56081 and [51]GSE70362 to analyze gene expression characteristics related to IVDD [[52]11, [53]12]. The gene expression levels across samples in these datasets were presented using boxplots. Principal Component Analysis (PCA) was employed to visualize the distribution of samples based on their gene expression profiles, revealing distinct patterns between the two datasets. After correcting for batch effects, the differences in expression levels between the datasets were reduced, as shown in the post-correction boxplots and PCA plots. A volcano plot illustrated the differential expression of genes between the datasets, marking significantly upregulated and downregulated genes. Hierarchical clustering of gene expression profiles for control and IVDD groups was visualized using a heatmap, highlighting potential biomarkers or therapeutic targets associated with IVDD. Weighted gene co-expression network analysis Weighted Gene Co-Expression Network Analysis (WGCNA) was performed to identify gene modules related to IVDD. Hierarchical clustering was used to organize samples based on their similarity, measured using Euclidean distance and average linkage criteria. Outliers, identified as samples branching off at higher heights in the dendrogram, were noted for their potential to affect downstream analyses. Dynamic module identification was conducted across different cohorts to detect groups of co-expressed and potentially co-regulated genes. The soft threshold power was determined to achieve a scale-free topology, with an optimal power of 6 providing a good fit. Module preservation was validated by evaluating the correlation between module membership and gene significance for IVDD, with specific focus on modules showing significant correlations. The relationships between module membership and gene significance were assessed to determine the relevance of identified gene modules to IVDD. Integrative analysis An integrative analysis was conducted to identify key genes and pathways associated with IVDD. A Venn diagram was used to display the overlap genes among differentially expressed genes (DEGs), efferocytosis-related genes (ECRG), and genes identified by WGCNA. LASSO regression modeling was employed to identify important predictors for IVDD, with optimization of the regularization parameter determined through cross-validation. Model performance was assessed using the Receiver Operating Characteristic (ROC) curve, which provided the area under the curve (AUC) as a measure of discrimination ability. Boxplots compared the risk scores of genes between control and IVDD groups, while a bar plot ranked genes by their importance in a random forest model. Feature selection was guided by cross-validation to determine the optimal number of genes for accurate IVDD classification. Pathway enrichment analysis was performed based on the selected genes, providing insights into the biological mechanisms underlying IVDD. Development and validation of a predictive model A predictive model for IVDD was developed and validated using multiple key genes as predictors. Nomogram analysis was conducted to estimate the risk of IVDD, with each predictor assigned a score based on its expression level. The total score, obtained by summing individual scores, corresponded to a predicted probability of IVDD risk. Calibration curves were used to assess the agreement between predicted probabilities and observed outcomes, with the mean absolute error indicating the model’s accuracy. Decision curve analysis (DCA) was used to evaluate the clinical utility of the predictive model by comparing the net benefit at different risk thresholds. The model’s discrimination ability was assessed using the ROC curve, with the AUC indicating the model’s performance. Additionally, correlation coefficients heatmaps were generated to show the relationships between key genes, and individual ROC curves were used to evaluate the discriminatory power of each gene. This comprehensive approach ensured the development of a robust and clinically useful predictive model for IVDD. Gene set enrichment analysis Gene Set Enrichment Analysis (GSEA) was performed to identify key pathways involved in IVDD. The analysis utilized the gene expression data to determine the enrichment of predefined gene sets within the dataset. Pathways were considered significantly enriched based on a p-value threshold, with adjustments made for multiple testing to control the false discovery rate. Specifically, the apoptosis pathway, inflammatory response pathway, reactive oxygen species (ROS) pathway, and TNFα signaling via NF-κB pathway were analyzed. Enrichment scores were calculated to quantify the degree of overrepresentation of these pathways in the IVDD samples compared to controls. Immune cell infiltration analysis Immune cell infiltration analysis was conducted to examine the involvement of immune cells in IVDD. The proportions of various immune cell types in the samples were estimated using computational algorithms designed for deconvoluting gene expression data into immune cell type fractions. A bar plot was generated to display the estimated proportions of different immune cell types across samples, providing a visual comparison between control and IVDD groups. A heatmap was used to illustrate the relative abundance of these immune cell types, highlighting differences in immune cell infiltration between the two groups. Additionally, scatter plots were created to assess the correlations between the expression levels of key genes (LGR6, THBS4, SERPINA1, ELMO1, and ITGB8) and the proportions of specific immune cell types. Nucleus pulposus isolation and culture Isolation and culture of human primary nucleus pulposus cells (NPCs) were conducted in accordance with Institutional Review Board approval and ethical guidelines. This study was approved by the Ethical Committee of Changzheng Hospital (Approval No. 2022SL217). Lumbar intervertebral disc samples were procured from patients, with no systemic diseases or infections, aged 18–60 years undergoing spinal surgery (patient 1, male, 37-year-old, L3/4 and L4/5, Pfirrmann grade II; patient 2, female, 47-year-old, L4/5, Pfirrmann grade II; patient 3, female, 41-year-old, L3/4, Pfirrmann grade II; patient 4, male, 52-year-old, L3/4, Pfirrmann grade II), with informed consent obtained from all participants. The samples were collected from January 2023 to June 2024 at Shanghai Changzheng Hospital. Exclusion criteria included systemic diseases or infections. For cell isolation, disc samples were placed in antibiotic-containing PBS post-surgery. The annulus fibrosus and cartilaginous endplates were excised to get the nucleus pulposus, which was then diced into 1 mm³ pieces and subjected to enzymatic digestion with 0.2% collagenase type II for 18–24 h at 37 °C and 5% CO2. Post-digestion, the tissue suspension was filtered through a 70-µm strainer, centrifuged at 300 g for 5 min, and the pellet resuspended in complete culture medium composed of Dulbecco’s Modified Eagle Medium/F-12 (DMEM/F-12) supplemented with 10% fetal bovine serum (FBS) and 1% Penicillin-Streptomycin. Cells were seeded in collagen type I-coated dishes and cultured under humidified conditions with 5% CO2 at 37 °C. Media changes were performed every 3–4 days to maintain cell purity. Upon reaching 80–90% confluence, cells were detached with trypsin-EDTA, counted, and subcultured at a density of 5,000 cells/cm² for subsequent experiments. Cell viability and proliferation were evaluated using trypan blue exclusion and MTT assays, respectively. Immunofluorescence staining for COL2A1 and ACAN confirmed cell identity. IVDD model establishment The disc puncture model was used to induce IVDD in mice [[54]13]. C57 BL/6 mice (8 weeks-old, male, 18–25 g) were purchased from Shanghai Research Center for Model Organisms. Mice were anesthetized using inhaled anesthetics like isoflurane, which were then placed in a prone position with the tail extended to expose the caudal vertebral column, specifically targeting the caudal intervertebral discs (C5-C6 or C6-C7) by palpation. The tail area was sterilized using 75% ethanol to minimize infection risk. A needle was then carefully inserted through the skin into the center of the intervertebral disc, penetrating about 3 mm to avoid complete disruption. The needle was rotated 180° and kept for 30s within the disc to cause degenerative alternations, then slowly withdrawn to prevent additional damage or leakage. The mice were monitored until full recovery from anesthesia. Mice were randomly assigned to groups in random order and then treated with pioglitazone (Piog, used for efferocytosis enhancement) by oral administration (10 mg/kg/day, for 8 weeks) [[55]14–[56]16], Annexin V (Ann, used for efferocytosis suppression) through micro-injection into the NP tissues (0.4 mg/kg/week, for 8 weeks) [[57]17, [58]18], or normal saline according to the groups (sham + NC, sham + Piog, IVDD + NC, IVDD + Piog; sham + NC, sham + Ann, IVDD + NC, IVDD + Ann; n = 6 per group). IVDD was assessed over 8 weeks using post-mortem histological analysis. The mice were euthanized following ethical guidelines, and the caudal intervertebral discs were harvested for further analysis. Each experiment was repeated independently three times to ensure statistical reliability. Co-culture of NPCs and THP-1 THP-1 cells, a human monocytic cell line, were cultured in RPMI 1640 medium supplemented with 10% FBS, 1% penicillin-streptomycin, and 1% L-glutamine at 37 °C in a humidified atmosphere containing 5% CO2. To induce differentiation into macrophages, THP-1 cells were treated with 100 ng/mL of phorbol 12-myristate 13-acetate (PMA) for 24 h [[59]19]. After differentiation, the medium was replaced with fresh RPMI 1640 medium without PMA, and the cells were incubated for an additional 24 h to allow for recovery and stabilization of the macrophage phenotype. Differentiated THP-1 macrophages were harvested using cold PBS and gently scraped from the culture dishes to avoid activation or damage. NPCs were seeded into 24-well plates at a density of 1 × 10^5 cells per well and allowed to adhere overnight. Each experimental condition was performed in triplicate, and co-culture experiments were conducted using three independent biological replicates to validate the findings. The following day, macrophages were added to the NPC cultures at a 1:1 ratio (1 × 10^5 macrophages per well) and co-cultured under standard conditions (37 °C, 5% CO2). NPCs were subjected to different experimental treatments prior to co-culture, to evaluate the effect of these treatments on macrophage-NPC interactions (n = 3 per group). To investigate the role of LGR6 in macrophages, LGR6 expression was manipulated in THP-1 macrophages using shRNA interference or overexpression techniques. The shRNA targeting LGR6 (LGR6 shRNA: CCUGGAACUGUCUCACAAUTT, AUUGUGAGACAGUUCCAGGTT) was synthesized and obtained from Genomeditech (Shanghai, China). For shRNA transfection, a final concentration of 30 nM shRNA was used in the culture medium. For LGR6 overexpression, macrophages were transfected with the LGR6 overexpression plasmid, also obtained from Genomeditech (Shanghai, China). And 5 µg of LGR6 overexpression plasmid was used per well in a 6-well plate. Macrophages were transfected using Lipofectamine 2000 (Invitrogen, USA) according to the manufacturer’s protocol. Assessment of efferocytosis using live-cell analysis The level of macrophage-mediated phagocytosis targeting apoptotic NPCs was assessed using a live-cell analysis (LCA) system. The co-cultures were monitored in real-time using the LCA system equipped with phase-contrast and fluorescence microscopy capabilities. Phagocytosis was quantified by measuring the internalization of fluorescently labeled apoptotic NPCs by macrophages. Macrophages and apoptotic NPCs were pre-labeled with DID (Cat#C1039, Beyotime, China) and YO-PRO-1 (Cat#C2022, Beyotime, China) before co-culture, respectively, and the dynamics of efferocytosis were recorded over a period of 24 h (n = 3 per group). The rate of efferocytosis was quantified as the percentage of macrophages containing internalized fluorescent NPCs. Reverse transcription quantitative polymerase chain reaction Reverse Transcription Quantitative Polymerase Chain Reaction (RT-qPCR) was performed to quantify gene expression in groups (n = 3 per group). Total RNA was extracted using a standardized extraction kit, with subsequent assessment of RNA purity and concentration using a spectrophotometer. For cDNA synthesis, 1 µg of total RNA was reverse-transcribed using a commercially available reverse transcription kit following the manufacturer’s protocol. Quantitative PCR was carried out on a real-time PCR system using SYBR Green dye for detection. Each reaction mixture contained 10 ng of cDNA, 0.2 µM of each primer, and SYBR Green PCR Master Mix in a total volume of 20 µL. Thermal cycling conditions included an initial denaturation at 95 °C for 10 min, followed by 40 cycles of 95 °C for 15 s and 60 °C for 1 min. Gene expression levels were calculated using the 2^^−ΔΔCt method, with β-actin or GAPDH serving as endogenous controls for normalization. Each sample was analyzed in triplicate, and no-template controls were included to ensure the absence of contamination. The specificity of amplification was confirmed by melt curve analysis. Statistical comparisons of gene expression levels between IVDD and control groups were performed using appropriate statistical tests to determine significance. The primers used in this study are presented in the Table [60]1. Table 1. The primers used in this study Genes Forward Reverse SERPINA1 GTGGTTTCTGAGCCAGCAG CCCTGTCCTCGTCCGTATT THBS4 GTTGCAGAACCTGGCATTCAG CCCTGGACCTGTCTTAGACTTCA LGR6 ACGGCTTACCTGGACCTCA GCTTGTCCTGGGATGTGTGA ELMO1 AGCATCTGGACAGTCAACACGG AGGCAGCGCCACAGAAATTTAC ITGB8 TGT GTG TGG AAG GTG TGA GTG TGT GAG GGT GAA GGC ATT G MMP13 ACTTTGTTGCCAATTCCAGG TTTGAGAACACGGGGAAGAC COL2A1 TGAGGGCGCGGTAGAGACCC TGCACACAGCTGCCAGCCTC BCL2 ACTTCGCCGAGATGTCC ATGACCCCACCGAACTC CASP3 GAGCACTGGAATGTCATCTCGCTCTG AGACCGAGATGTCATTCCAGTGCTT BAX CCCGAGAGGTCTTTTTCCGAG CCAGCCCATGATGGTTCTGAT MERTK AGA TCG TGT CTG ATC CCA TCT GAG GTT GAA GGC TGT GTT TCT AXL CAGCGCAGCCTGCATGT TTGGCGTTATGGGCTTCG TYRO3 GAGGATGGGGGTGAAACC ACTGTGAAAAATGGCACACCT CX3CR1 TGACTGGCAGATCCAGAGGTT GTAGAATATGGACAGGAACAC GAPDH AGAAGGCTGGGGCTCATTTG AGGGGCCATCCACAGTCTTC ACTB AAGGTGACAGCAGTCGGTT TGTGTGGACTTGGGAGAGG [61]Open in a new tab Western blotting Western blotting was conducted to analyze protein expression in groups (n = 3 per group). Cells were lysed in RIPA buffer supplemented with protease and phosphatase inhibitors, while tissue samples were homogenized in the same buffer. Lysates were centrifuged at 12,000 g for 20 min at 4 °C to remove cellular debris, and supernatants were collected. Protein concentrations were determined using a bicinchoninic acid (BCA) protein assay kit. Equal amounts of protein (30 µg for cells and 50 µg for tissues) were separated by SDS-PAGE on 10% polyacrylamide gels and transferred onto polyvinylidene fluoride (PVDF) membranes. Membranes were blocked with 5% non-fat dry milk in Tris-buffered saline with 0.1% Tween-20 (TBST) for 1 h at room temperature. Primary antibodies specific to target proteins and the loading control (β-actin) were incubated with the membranes overnight at 4 °C at appropriate dilutions. Following washes with TBST, membranes were incubated with horseradish peroxidase (HRP)-conjugated secondary antibodies for 1 h at room temperature. Protein bands were visualized using an enhanced chemiluminescence (ECL) detection system and quantified by densitometry using ImageJ software. Relative protein expression levels were normalized to β-actin. Statistical analysis was performed to determine significant differences in protein expression between IVDD and control groups, with appropriate statistical tests used to assess these differences. Primary antibodies used in this study included LGR6 (Cat#17658-1-AP, Proteintech, China), THBS4 (Cat#ab263898, Abcam, UK), SERPINA1 (Cat#16382-1-AP, Proteintech, China), ELMO1 (Cat#ab174298, Abcam, UK), ITGB8 (Cat#ab317020, Abcam, UK), GAPDH (Cat#60004-1-Ig, Proteintech, China), COL2A1 (Cat#AF0135, Affinity, China), MMP13 (Cat#18165-1-AP, Proteintech, China), Bcl2 (Cat#AF6139, Affinity, China), Cleaved-Caspase 3 (Cat#AF7022, Affinity, China), Bax (Cat#AF0120, Affinity, China), MERTK (Cat#DF7344, Affinity, China), AXL (Cat#AF7793, Affinity, China), TYRO3 (Cat#DF10363, Affinity, China), CX3CR1 (Cat#DF7096, Affinity, China). Cell counting Kit-8 assay The cell viability was assessed using the Cell Counting Kit-8 (CCK8, Dojindo, Japan) according to the manufacturer’s instructions. Briefly, cells were seeded into 96-well plates at a density of 5 × 10⁴ cells per well in 100 µL of complete culture medium and incubated overnight at 37 °C in a humidified atmosphere containing 5% CO₂. After the experimental treatments, 10 µL of CCK8 reagent was added to each well, followed by incubation for 2 h under the same conditions. The absorbance was measured at 450 nm using a microplate reader (Bio-Rad). Wells without cells served as blanks for background absorbance subtraction. The results were normalized to the untreated control group, and cell viability was calculated as a percentage relative to the control. Immunofluorescence analysis for cells and tissue samples Immunofluorescence analysis was performed to investigate the localization and expression of target proteins in both cells (n = 3 per group) and IVDD tissue samples (n = 6 per group). Cells were cultured on coverslips and fixed with 4% paraformaldehyde for 15 min, followed by permeabilization with 0.1% Triton X-100 for 10 min. Tissue samples were fixed in 10% formalin, embedded in paraffin, and sectioned into 4-µm thick slices. Sections were deparaffinized in xylene, rehydrated through graded ethanol, and subjected to antigen retrieval by heating in citrate buffer (pH 6.0) for 20 min. Both cells and tissue sections were blocked with 5% bovine serum albumin (BSA) for 30 min to prevent non-specific binding. Primary antibodies specific to the target proteins were applied and incubated overnight at 4 °C, which included LGR6 (Cat# 17658-1-AP, Proteintech, China), Bcl2 (Cat#AF6139, Affinity, China), COL2A1 (Cat#AF0135, Affinity, China), and MMP13 (Cat#18165-1-AP, Proteintech, China). After washing with phosphate-buffered saline (PBS), samples were incubated with fluorophore-conjugated secondary antibodies for 1 h at room temperature in the dark. Nuclei were counterstained with DAPI for 5 min. Coverslips and tissue sections were mounted with an anti-fade mounting medium and examined under a fluorescence microscope. Images were captured and analyzed using image analysis software to quantify the fluorescence intensity. The relative expression levels and subcellular localization of the target proteins were compared between IVDD and control groups. Statistical analysis was performed to determine significant differences in protein expression and localization, with appropriate statistical tests employed to ensure rigor and reproducibility. Statistical analysis All statistical analyses were performed using GraphPad Prism software (version 10). For comparisons between two groups, a two-tailed unpaired Student’s t-test was employed to assess statistical significance. For multiple group comparisons, one-way ANOVA followed by post-hoc Tukey’s test was utilized to determine significance across different experimental conditions. Each experiment was repeated at least three times independently to confirm the consistency of the findings. The level of statistical significance was set at P < 0.05 for all analyses. Significance levels in the figures were denoted as follows: P < 0.05 (*), P < 0.01 (**), P < 0.001 (***), and P < 0.0001 (****). Non-significant results were denoted as “ns.” Results Gene expression analysis of [62]GSE56081 and [63]GSE70362 [64]GSE56081 and [65]GSE70362 were used to analyze the gene expression characteristics related to IVDD. The expression levels of genes across samples in datasets ([66]GSE56081 and [67]GSE70362) were presented in the boxplot. The boxplot distribution indicated that the expression levels were generally higher in [68]GSE56081 compared to [69]GSE70362 (Fig. [70]1A). PCA revealed the distribution of samples based on the gene expression profiles in the primary dataset. The [71]GSE56081 samples (blue) were clearly separated from the [72]GSE70362 samples (orange), suggesting distinct gene expression patterns between the two groups (Fig. [73]1B). The boxplot presented gene expression levels in dataset [74]GSE56081 (blue) and [75]GSE70362 (orange) after batch effect correction and the results showed that post-correction, the difference in expression levels between the two datasets was reduced (Fig. [76]1C). PCA plot of samples after batch effect correction showed that after correction, there was more overlap between the samples from the two datasets, suggesting that the batch effects have been reduced, though some separation still exists (Fig. [77]1D). Volcano plot illustrated the differential expression of genes between the two datasets ([78]GSE56081 and [79]GSE70362). Genes that were significantly upregulated are marked in orange, and downregulated genes were marked in blue. The volcano plot revealed a considerable number of differentially expressed genes between the two datasets, indicating potential biological significance (Fig. [80]1E). The heatmap showed hierarchical clustering of gene expression profiles for control (blue) and IVDD groups (orange). The expression levels of multiple genes across samples were presented in the heatmap, with a gradient from blue (low expression) to orange (high expression). The clustering pattern revealed distinct gene expression signatures between the control and IVDD groups, suggesting potential biomarkers or therapeutic targets associated with IVDD (Fig. [81]1F). Fig. 1. [82]Fig. 1 [83]Open in a new tab Analysis of gene expression in IVDD datasets [84]GSE56081 and [85]GSE70362. (A) Boxplot of gene expression levels across samples in datasets [86]GSE56081 and [87]GSE70362. (B) PCA plot of gene expression profiles in the primary dataset. (C) Boxplot of gene expression levels in datasets [88]GSE56081 and [89]GSE70362 after batch effect correction. (D) PCA plot of samples after batch effect correction. Increased overlap between the datasets suggests reduced batch effects, with some separation persisting. (E) Volcano plot of differential gene expression between [90]GSE56081 and [91]GSE70362. (F) Heatmap of hierarchical clustering of gene expression profiles for control and IVDD groups. The clustering pattern indicates distinct gene expression signatures between the control and IVDD groups WGCNA for identifying gene modules related to IVDD The clustering was performed using hierarchical methods, which organize samples into a tree-like structure based on their similarity. The similarity was usually measured using a distance metric (such as Euclidean distance) and a linkage criterion (such as average linkage). In the dendrogram, each leaf (terminal node) represents a sample, and branches represent the clustering process. The height at which two samples (or clusters) were joined indicates their similarity, with lower heights indicating higher similarity. Outliers were identified as samples that branch off at a higher height, indicating they were less similar to the majority of the samples. These outliers could be due to technical variations, data entry errors, or biological differences. Detecting outliers was crucial because they can skew the results of downstream analyses, such as module detection in WGCNA or differential expression analysis. By identifying and potentially removing or adjusting for outliers, the integrity and reliability of the analysis were improved (Fig. [92]2A). Dynamic module identification was performed in the different cohorts (Fig. [93]2B). Modules are thought to represent groups of genes that are co-expressed and potentially co-regulated, indicating they may participate in similar biological pathways or processes. The preservation of these modules can further validate their biological relevance. Highly preserved modules are more likely to represent fundamental biological processes. The soft threshold power in the cohort was 5 (Fig. [94]2C). The plot demonstrated the scale-free topology fit index as a function of the soft-thresholding power (Fig. [95]2C). The network achieved a scale-free topology at a power of 6, where the R^2 value was close to 0.9, indicating a good fit. The plot of mean connectivity showed mean connectivity as a function of the soft-thresholding power. As the power increases, mean connectivity decreases, confirming the scale-free nature of the network. For the cohort, 8 co-expression modules were clustered, with the brown module (Cor = -0.22 vs. 0.22, P = 0.09) and green module (Cor = -0.34 vs. 0.34, P = 0.009) having the stronger positive correlations with control score and IVDD score (Fig. [96]2D). The relationships between module membership and gene significance for IVDD within the green module and brown module were presented. The correlation coefficient of 0.11 (P = 0.00046) indicated a mild positive relationship between module membership and gene significance (Fig. [97]2E). This meant that genes with higher connectivity within the green module tended to have higher significance for IVDD. The correlation coefficient of -0.033 (P = 0.24) indicated a possible negative relationship between module membership and gene significance. However, in this case, there was no significant relationship between these metrics, implying that the genes in the brown module were not particularly relevant to IVDD (Fig. [98]2F). Fig. 2. [99]Fig. 2 [100]Open in a new tab Analysis of sample clustering, module identification, and their relationships with IVDD scores. (A) Hierarchical clustering dendrogram of samples based on gene expression similarity. Each leaf represents an individual sample, and the height at which samples are joined indicates their similarity. Outliers are identified as samples branching off at higher heights, suggesting lower similarity to the majority of samples. (B) Dynamic module identification across different cohorts. Modules represent clusters of genes that are co-expressed and potentially co-regulated, suggesting involvement in similar biological processes. (C) Selection of soft threshold power for network construction. (D) Clustering of 8 co-expression modules in the cohort. (E) Relationship between module membership and gene significance for IVDD within the green module. (F) Relationship between module membership and gene significance for IVDD within the brown module Integrative analysis identifying key genes and pathways associated with IVDD The Venn diagram displayed the overlap between different sets of genes: DEGs, ECRG, and genes from the green module identified by WGCNA green. The intersections highlighted genes common to multiple datasets, emphasizing their potential importance in IVDD. Notably, 11 genes were common across all three sets, suggesting they are key players in the disease process (Fig. [101]3A). The LASSO regression model showed the trajectories of the LASSO (Least Absolute Shrinkage and Selection Operator) regression coefficients for various genes as the regularization parameter (lambda) changes. Genes with non-zero coefficients were selected as important predictors for IVDD (Fig. [102]3B). The results of cross-validation for the LASSO model showed that the optimal lambda, marked by the vertical dashed line, minimizes the cross-validation error, balancing model complexity and predictive accuracy. Lambda value was used to select the most relevant genes for predicting IVDD, contributing to a parsimonious and robust predictive model. The optimal lambda value was determined by the point where the binomial deviance is minimized. The minimum binomial deviance Log(lambda) was around “-3”, which corresponded to the value where the red dots and error bars were at their lowest point (Fig. [103]3C). The ROC curve assessed the performance of the prediction model. The AUC was 0.919, indicating excellent discrimination between IVDD and control samples (Fig. [104]3D). And the high sensitivity and specificity values suggested that the model is robust in distinguishing between the two groups (Fig. [105]3D). The boxplot compared the risk score of genes between the control group and the IVDD group. Higher risk score in the IVDD group suggested it may be involved in the disease mechanism. The significant difference underscored the gene’s potential role in IVDD (P < 0.001) (Fig. [106]3E). The bar plot ranked genes by their importance in the random forest model, using the mean decrease in Gini index as a metric. This analysis helped prioritize genes for further study based on their predictive power. Genes with higher values, such as LGR6, THBS4, SERPINA1, ELMO1, ITGB8, and IL6R, were critical in classifying samples as IVDD or control (Fig. [107]3F). The plots in Fig. [108]3G show the accuracy and error rates of the model with different numbers of features (genes) selected by cross-validation. The left plot indicated that accuracy peaks at 9 features, while the right plot shows the corresponding error rate decreases (Fig. [109]3G), suggesting an optimal set of 9 genes for the most accurate classification of IVDD. Then the Venn diagram was used to illustrate the overlap of genes selected by different feature selection methods: RF, svm_RFE, and LASSO. Five genes were common across all methods, highlighting their robustness and potential as key biomarkers for IVDD (Fig. [110]3H). The comparison of expression levels of multiple key genes, including LGR6, THBS4, SERPINA1, ELMO1, and ITGB8, between control and IVDD groups were carried out. The significant differences in expression levels suggested that when compared with those in the control group, these genes in IVDD group were differentially expressed, reinforcing their potential role in IVDD (Fig. [111]3I). In addition, the results of pathway enrichment analysis for the selected genes showed that many pathways were significantly enriched, including ECM-receptor interaction, regulation of actin skeleton, efferocytosis, and phagosome, indicating that these pathways may play critical roles in the pathogenesis of IVDD and could be targets for therapeutic intervention (Fig. [112]3J). Fig. 3. [113]Fig. 3 [114]Open in a new tab Analysis of gene overlap, predictive modeling, and pathway enrichment in IVDD. (A) Venn diagram illustrating the overlap between DEGs, ECRG, and genes from the WGCNA green module. Eleven genes are common to all three sets, indicating their potential significance in IVDD. (B) LASSO regression model showing the trajectories of regression coefficients for various genes as the regularization parameter (lambda) changes. (C) Cross-validation results for the LASSO model, with the optimal lambda value indicated by the vertical dashed line, minimizing cross-validation error and binomial deviance. (D) ROC curve and corresponding area under the curve for the prediction model, demonstrating excellent discrimination between IVDD and control samples. (E) Boxplot comparing gene risk scores between control and IVDD groups. (F) Bar plot ranking genes by their importance in the random forest model, with the mean decrease in Gini index as a metric. Genes such as LGR6, THBS4, SERPINA1, ELMO1, ITGB8, and IL6R are prioritized for further study. (G) Plots showing the accuracy and error rates of the model with varying numbers of features selected by cross-validation. (H) Venn diagram illustrating the overlap of genes selected by RF, svm_RFE, and LASSO methods. (I) Comparison of expression levels for key genes (LGR6, THBS4, SERPINA1, ELMO1, ITGB8) between control and IVDD groups, showing significant differences that reinforce their potential role in IVDD. (J) Pathway enrichment analysis results for the selected genes, showing significant enrichment in pathways Development and validation of a predictive model for IVDD Multiple predictors (LGR6, THBS4, SERPINA1, ELMO1, ITGB8) was used to estimate the risk of IVDD by nomogram analysis. Each predictor was assigned a score based on its expression level. The total score, obtained by summing individual scores, corresponded to a predicted probability of IVDD risk. Lower level of LGR6, THBS4, or SERPINA1 was relevant to higher IVDD risk, and higher level of ELMO1 or ITGB8 was closely associated with higher risk of IVDD (Fig. [115]4A). The calibration curve assessed the agreement between predicted probabilities and observed outcomes. The solid line represented the bias-corrected model, while the dashed line represented the apparent model. The closer these lines are to the ideal 45-degree line, the better the model’s predictive accuracy. The mean absolute error of 0.036 suggested good calibration (Fig. [116]4B). The decision curve analysis (DCA) evaluated the clinical utility of the predictive model by comparing the net benefit at different risk thresholds. The net benefit of using the model to predict IVDD was higher than those of single gene, especially at risk thresholds between 0.1 and 0.7 (Fig. [117]4C). This indicated that the model provides a meaningful benefit for clinical decision-making within this range. The ROC curve evaluated the model’s ability to discriminate between IVDD and non-IVDD cases. The AUC of 0.898 (95% CI: 0.812–0.985) indicated excellent discrimination of the model performance (Fig. [118]4D). The correlation coefficients heatmap showed the relationship of key genes (LGR6, THBS4, SERPINA1, ELMO1, ITGB8). Positive correlations were indicated by shades of red, and negative correlations by shades of blue. Significant positive correlations were observed between several pairs of genes, suggesting potential co-regulation or shared pathways (Fig. [119]4E). ROC curves for individual key genes revealed that the key genes show good discriminatory power (LGR6: 0.7928; THBS4: 0.7824; SERPINA1: 0.8275; ELMO1: 0.7568; ITGB8: 0.749), indicating their effective prediction for IVDD (Fig. [120]4F). Fig. 4. [121]Fig. 4 [122]Open in a new tab Development and validation of a predictive model for IVDD. (A) Nomogram for estimating the risk of IVDD using multiple predictors (LGR6, THBS4, SERPINA1, ELMO1, ITGB8). (B) Calibration curve assessing the agreement between predicted probabilities and observed outcomes. (C) DCA demonstrating the net benefit of using the predictive model at different risk thresholds. (D) ROC curve evaluating the model’s ability to discriminate between IVDD and non-IVDD cases. (E) Heatmap of correlation coefficients showing the relationships between key genes (LGR6, THBS4, SERPINA1, ELMO1, ITGB8). (F) ROC curves for individual key genes, with AUC values indicating their discriminatory power GSEA of key pathways in IVDD The GSEA result showed that the apoptosis pathway is significantly enriched with a p-value of 0.0203 and an adjusted p-value of 0.01691, suggesting that apoptosis is actively involved in IVDD (Fig. [123]5A). The significant enrichment of apoptosis-related genes indicated that programmed cell death plays a crucial role in IVDD, which could contribute to the loss of disc cells and structural integrity. The enrichment of inflammatory response pathway was statistically significant, highlighting the role of inflammation in IVDD pathogenesis (Fig. [124]5B). The strong enrichment of inflammatory response genes underscored the importance of inflammation in the pathogenesis of IVDD. Chronic inflammation could lead to pain and further tissue injury. The reactive oxygen species (ROS) pathway showed high enrichment score, indicating oxidative stress involvement in IVDD (Fig. [125]5C). The enrichment of ROS-related genes suggested that oxidative stress was a significant factor in IVDD. Oxidative damage to cells and extracellular matrix components could accelerate disc degeneration. The TNFα signaling via NF-κB pathway was also highly enriched, underscoring the critical role of TNFα signaling-related inflammation in IVDD (Fig. [126]5D). Fig. 5. [127]Fig. 5 [128]Open in a new tab GSEA of key pathways in IVDD. (A) Enrichment of the apoptosis pathway. (B) Enrichment of the inflammatory response pathway. (C) Enrichment of the ROS pathway. (D) Enrichment of the TNFα signaling via NF-κB pathway Immune cell infiltration analysis in IVDD The bar plot showed the estimated proportion of various immune cell types across different samples. Each bar represented a sample, and the colors within each bar indicate the proportion of different immune cell types. The distribution of immune cell types varied among samples, with certain cell types like macrophages and mast cells showing more prominent proportions (Fig. [129]6A). The heatmap illustrated the relative abundance of different immune cell types in control and IVDD samples. The differential abundance of cells in the control group and the IVDD group suggested that specific immune cell types such as macrophages may play roles in the pathophysiology of IVDD, potentially contributing to inflammation and tissue degeneration (Fig. [130]6B). The scatter plots showed the correlation between the expression of key genes (LGR6, THBS4, SERPINA1, ELMO1, ITGB8) and various immune cell types (Fig. [131]6C). LGR6 was negatively correlated with macrophages M0. THBS4 and SERPINA1 were positively correlated with CD8 T cells (SERPINA1 is also negatively correlated with memory B cells). ELMO1 was positively correlated with macrophages M0 but negatively with T cells follicular helper. And ITGB8 was negatively correlated with plasma cells. The correlations suggested that these key genes may be involved in modulating the immune response in IVDD. Fig. 6. [132]Fig. 6 [133]Open in a new tab Analysis of immune cell infiltration in intervertebral disc degeneration (IVDD). (A) Bar plot depicting the estimated proportion of various immune cell types across different samples. Each bar represents an individual sample, and the color segments within each bar indicate the proportion of different immune cell types. (B) Heatmap illustrating the relative abundance of different immune cell types in control and IVDD samples. (C) Scatter plots revealing the correlation between the expression of key genes (LGR6, THBS4, SERPINA1, ELMO1, ITGB8) and various immune cell types. These correlations imply that the key genes may be involved in regulating the immune response in IVDD LGR6 is significantly upregulated in THP-1 co-cultured with IL-1β-treated NPCs According to the above results, efferocytosis might play an important role in modulating IVDD. IL-1β is commonly used to treat the NPCs to induce inflammation, cell death, and degenerative changes [[134]20, [135]21]. Firstly, in this present study, IL-1β was used to induce NPCs apoptosis. The results of RT-qPCR demonstrated a dose-dependent increase in BAX and CASP3 mRNA levels, coupled with a decrease in BCL2 mRNA expression (Fig. [136]7A-C). Western blot analysis confirmed the corresponding changes at the protein level of these genes. There was a marked increase in the protein expression of both BAX and cleaved caspase 3, alongside a reduction in BCL2 protein level, indicating apoptosis activation (Fig. [137]7D-G). The TUNEL assay, which quantitatively measures apoptosis by detecting DNA fragmentation, revealed a significant increase in TUNEL-positive cells with increasing concentrations of IL-1β, confirming the induction of apoptosis in a dose-dependent manner (Fig. [138]7H-I). These findings above emphasized the pro-apoptotic effects of IL-1β. Fig. 7. [139]Fig. 7 [140]Open in a new tab IL-1β-induced apoptosis in NPCs and subsequent changes in efferocytosis-related gene expression in THP-1-derived macrophages. (A-C) IL-1β induces a dose-dependent shift towards a pro-apoptotic environment in NPCs. Data represent mean ± SD (n=3; one-way ANOVA) of relative mRNA levels for BCL2 (A), BAX (B), and CASP3 (C) following treatment with various concentrations of IL-1β. (D-G) Western blot analysis and quantification of protein levels for BCL2, BAX, and cleaved caspase 3 in NPCs treated with IL-1β(n=3; one-way ANOVA). (H-I) TUNEL assay reveals a significant increase in apoptotic cells following IL-1β treatment (n=3; one-way ANOVA). (J-N) mRNA levels of efferocytosis-related genes SERPINA1 (J), THBS4 (K), ELMO1 (L), LGR6 (M), and ITGB8 (N) in THP-1-derived macrophages co-cultured with IL-1β-treated NPCs (n=3; one-way ANOVA). (O-T) Western blot analysis and quantification of protein levels for SERPINA1 (O-P), THBS4 (O-Q), ELMO1 (O-R), LGR6 (O-S),and ITGB8 (O-T) in THP-1-derived macrophages co-cultured with IL-1β-treated NPCs (n=3; one-way ANOVA). Data are presented as mean ± SD. Statistical significance is denoted as follows: “*”, p < 0.05; “**”,p < 0.01; “***”, p < 0.001; “****”, p < 0.0001; “ns”, no significance. Then, macrophages were cocultured with IL-1β-treated NPCs. The mRNA and protein levels of genes potentially associated with macrophage efferocytosis, including SERPINA1, THBS4, ELMO1, LGR6, and ITGB8 were firstly detected. The relative mRNA levels of SERPINA1 showed no significant differences across various concentration (Fig. [141]7J). There was a significant decrease in the mRNA expression of THBS4 in the 50 ng/mL group compared to 0 ng/mL group, but no significant change was observed between 0 ng/mL and 100 ng/mL groups (Fig. [142]7K). The mRNA level of ELMO1 remained unchanged across all groups, suggesting that ELMO1 mRNA expression was not significantly influenced by IL-1β treatment (Fig. [143]7L). LGR6 showed a significant increase in mRNA level in 100 ng/mL group compared to both the 0 ng/mL group and 50 ng/mL group, indicating a strong association between LGR6 expression in macrophages and IVDD progression (Fig. [144]7M). ITGB8 mRNA levels did not exhibit significant changes (Fig. [145]7N). Furthermore, the protein levels of these genes were assessed through Western blot analysis. SERPINA1 protein level in the 100 ng/mL group was slightly decreased compared to that in the 0 ng/mL group (Fig. [146]7O-P). THBS4 protein levels also did not significantly differ, despite the mRNA decrease in the moderate condition, suggesting potential post-transcriptional regulation (Fig. [147]7O, Q). Although significant changes in ELMO1, LGR6, and ITGB8 protein levels were observed only in 100 ng/mL group, LGR6 showed 3.5-fold increase, higher than ELMO1 and ITGB8 (Fig. [148]7O, R-T). According to the above results, LGR6 showed the most significant increase in high concentration of IL-1β group, indicating its potential upregulation as IVDD progresses. Thus, we chose LGR6 gene for the following experiments. The effects of LGR6 overexpression and knockdown on macrophages efferocytosis According to the findings above, the role of LGR6 on efferocytosis of macrophages was explored. Macrophages were co-cultured with NPCs that were treated with or without IL-1β to assess the efferocytosis activity of macrophages. LGR6 in macrophages was overexpressed and the expression level of LGR6 was assessed. The results showed that both the mRNA level and protein level of LGR6 were significantly raised after overexpression with co-culturing with NPCs (Fig. [149]8A, F-G), demonstrating that the LGR6 gene was successfully overexpressed. Then, the mRNA levels of other efferocytosis-related genes, including MERTK, AXL, and TYRO3, and macrophage marker CX3CR1 were detected and the results showed that these indicators were significantly upregulated, which were especially raised with both LGR6 overexpression (IL-1β + OE-LGR6 group) compared to the IL-1β + OE-NC group (Fig. [150]8B-E), implying that LGR6 overexpression in an inflammatory environment enhances efferocytosis. The results of western blot showed that IL-1β-treated NPCs co-culturing enhanced the expression of MERTK, AXL, TYRO3, and CX3CR1 in macrophages, and LGR6 overexpression in macrophages further enhances protein levels under IL-1β stimulation, emphasizing the impact of LGR6 on efferocytosis (Fig. [151]8F, H-K). Then, the macrophage-mediated efferocytosis targeting apoptotic NPCs was assessed using a LCA system. The results showed that IL-1β treatment in NPCs, at different time points (1 h, 12 h, 18 h, 24 h), significantly enhanced the efferocytosis level, indicating that efferocytosis was activated by IL-1β-induced apoptosis (Fig. [152]8L-M). In addition, the IL-1β-induced raise of efferocytosis level was further enhanced by LGR6 overexpression at 12 h, 18 h, and 24 h, implying a synergistic interaction between LGR6 and inflammatory signals in promoting macrophage efferocytosis (Fig. [153]8L-M). Fig. 8. [154]Fig. 8 [155]Open in a new tab Effects of LGR6 on efferocytosis in the context of inflammation-induced NPCs apoptosis. (A-E) Relative mRNA expression levels of LGR6, MERTK, AXL, TYRO3, and CX3CR1 in macrophages (n = 3; one-way ANOVA). (F-K) Western blot and quantification of protein expression levels for MERTK, AXL, TYRO3, and CX3CR1. (L) The representative images of efferocytosis in macrophages co-cultured with NPCs in OE-NC or OE-LGR6 models over time (1 h, 12 h, 18 h, 24 h) (n = 3; one-way ANOVA). (M) Quantification of efferocytosis at different time points. (N-S) Western blot and protein quantification of LGR6, MERTK, AXL, TYRO3, and CX3CR1 (n = 3; one-way ANOVA). (T) The representative images of efferocytosis in macrophages co-cultured with NPCs in sh-NC or sh-LGR6 models over time (1 h, 12 h, 18 h, 24 h) (n = 3; one-way ANOVA). (U) Quantification of efferocytosis at different time points (n = 3; one-way ANOVA). Data are presented as mean ± SD. Statistical significance is denoted as follows: “*”, p < 0.05; “**”, p < 0.01; “***”, p < 0.001; “****”, p < 0.0001; “ns”, no significance To validate the role of LGR6 in regulating macrophage efferocytosis, LGR6 expression in macrophages was inhibited. Western blot results showed that LGR6 knockdown (sh-LGR6) significantly lowered the protein level of LGR6 with or without IL-1β-stimulation in NPCs (Fig. [156]8N-O). In addition, the results showed that LGR6 downregulation reduced the expression of MERTK, AXL, TYRO3, and CX3CR1 with the most pronounced effects observed under IL-1β stimulation (Fig. [157]8N, P-S), underscoring LGR6’s role in regulating efferocytosis during inflammatory microenvironment. Efferocytotic activity under LGR6 knockdown with and without IL-1β exposure over time (1 h, 12 h, 18 h, 24 h) was further assessed through LCA system. The results showed that under IL-1β treatment LGR6 knockdown significantly decreased efferocytotic levels compared to control at 12 h, 18 h, and 24 h (Fig. [158]8T-U), highlighting its crucial role in macrophage efferocytosis. The role of macrophage LGR6 in modulating ECM degradation and NPCs apoptosis To assess the effect of regulating macrophage LGR6 on NPCs, anabolism and catabolism of ECM were detected. The result showed that MMP13 mRNA level in NPCs did not significantly change with macrophage LGR6 overexpression alone but was significantly upregulated with IL-1β treatment in NPCs, which was lowered in the IL-1β + OE-LGR6 group, indicating that LGR6 overexpression in macrophages lowered the MMP13 mRNA level in NPCs under inflammatory situation (Fig. [159]9A). COL2A1 mRNA level was significantly decreased with IL-1β treatment, which was reversed by LGR6 overexpression, suggesting a protective role of macrophage LGR6 overexpression in regulating IVDD progress (Fig. [160]9B). The increased MMP13 protein level and reduced COL2A1 protein level induced by IL-1β stimulation were reversed by LGR6 overexpression (Fig. [161]9F, G-H), consistent with the mRNA results. And the immunofluorescence results showed that under IL-1β exposure, LGR6 overexpression raised the IF intensity of COL2A1, suggesting the ECM protection effect of macrophage LGR6 overexpression for NPCs (Fig. [162]9L-M). Fig. 9. [163]Fig. 9 [164]Open in a new tab The impact of LGR6 overexpression and knockdown on ECM degradation and NPCs apoptosis. (A-E) Relative mRNA expression levels of MMP13, COL2A1, BCL2, Cleaved Caspase 3, and BAX in overexpression models (n = 3; one-way ANOVA). (F-K) Western blot and quantification of protein expression levels for MMP13, COL2A1, BCL2, cleaved Caspase 3, and BAX in overexpression models (n = 3; one-way ANOVA). (L) Immunofluorescence results of BCL2 and COL2A1 in OE-NC and OE-LGR6 models, both with and without IL-1β treatment (n = 3; one-way ANOVA). (M-N) Relative fluorescence intensity of BCL2 and COL2A1 (n = 3; one-way ANOVA). (O-P) TUNEL assay and quantification (n = 3; one-way ANOVA). (Q-V) Western blot and quantification of protein levels for, MMP13, COL2A1, BCL2, Cleaved Caspase 3, and BAX (n = 3; one-way ANOVA). Data are presented as mean ± SD. Statistical significance is denoted as follows: “*”, p < 0.05; “**”, p < 0.01; “***”, p < 0.001; “****”, p < 0.0001; “ns”, no significance In addition, apoptosis was assessed in NPCs. Compared to the IL-1β + OE-NC group, BCL2 mRNA expression increased significantly in the IL-1β + OE-LGR6 group (Fig. [165]9C). The mRNA levels of CASP3 and BAX with significant upregulation were observed under IL-1β stimulation, and OE-LGR6 lowered IL-1β-induced upregulation of CASP3 and BAX (Fig. [166]9D-E). Under IL-1β stimulation, BCL2 protein levels were significantly elevated with LGR6 overexpression, indicating enhanced anti-apoptotic signaling (Fig. [167]9F, I). In contrast, IL-1β stimulation-induced increase in both cleaved caspase 3 and BAX proteins were significantly downregulated by LGR6 overexpression (Fig. [168]9F, J-K). What’s more, TUNEL assay showed a marked increase in TUNEL positive cells in the IL-1β-treated group, which was reversed by LGR6 overexpression (Fig. [169]9O-P). Besides, relative IF intensity measurement for BCL2 revealed increased anti-apoptotic signals with LGR6 overexpression (Fig. [170]9L, N). The above results indicated a potential anti-apoptotic effect by LGR6 overexpression in macrophages. Then, to confirm the role of LGR6 in regulating ECM and apoptosis in NPCs, LGR6 was downregulated in NPCs. Western blot results show that the knockdown of LGR6 results in increased MMP13 and decreased COL2A1 protein levels under IL-1β stimulation, suggesting protective effects against matrix degradation (Fig. [171]9Q, R-S). The knockdown of LGR6 under IL-1β condition resulted in raised level of pro-apoptotic protein levels (cleaved caspase 3 and BAX) and lowered anti-apoptotic protein level (BCL2) (Fig. [172]9Q, T-V). The above results emphasize the dual role of LGR6 in suppressing apoptosis and protecting ECM from degradation. Comprehensive analysis of the effects of modulating macrophage efferocytosis on IVDD To investigate the effect of modulating macrophage efferocytosis on IVDD, efferocytosis enhancers (pioglitazone, Piog) and suppressors (Annexin V, Ann) were used to treat cells and mice [[173]15–[174]18]. We performed the CCK8 assay to examine the cell viability of macrophage treated with Piog or Ann at a set concentration gradient. The results showed that the pioglitazone at the concentration of 1 µM significantly lowered the cell viability and the Annexin V at the concentration of 300 µM significantly lowered the cell viability (Supplementary Fig. [175]1A-B). Then the final concentrations of pioglitazone and Annexin V used to treat macrophages were 0.5 µM and 200µM, respectively. The mRNA levels of MERTK and AXL were significantly increased in the IL-1β group compared to the control group; in addition, the mRNA levels of MERTK and AXL in the IL-1β + Piog group were further significantly raised, compared to the IL-1β group, suggesting that pioglitazone could enhance efferocytosis (Fig. [176]10A-B). Protein levels of MERTK and AXL were significantly elevated in the IL-1β + Piog group compared to the IL-1β group (Fig. [177]10E-G), consistent with mRNA findings. The above findings demonstrated that efferocytosis was promoted by pioglitazone treatment. In addition, the effect of efferocytosis enhancement on IVDD was assessed. The mRNA levels of both COL2A1 and MMP13 in the Piog group were not significantly changed compared to those in the Ctrl group; however, mRNA level of COL2A1 was significantly elevated and MMP13 was significantly lowered in the IL-1β + Piog group compared to those in the IL-1β group, indicating that pioglitazone could promote matrix synthesis and suppress catabolic enzyme expression in degenerative discs (Fig. [178]10C-D). Similarly, an increase in COL2A1 protein level in the IL-1β + Piog group further supported the protective role of pioglitazone in enhancing matrix synthesis (Fig. [179]10E, H). MMP13 protein level was reduced in the IL-1β + Piog group, aligning with the mRNA results indicating reduced matrix degradation (Fig. [180]10E, I). HE and SOFG staining showed no damaged effect on discs in the sham + Piog group compared to the sham + NC group (Fig. [181]10J). In addition, the Piog + IVDD group showed improved structural integrity compared to IVDD + NC group (Fig. [182]10J). The histological score was significantly better in the Piog-treated IVDD group than that in the IVDD group, indicating reduced degeneration (Fig. [183]10K). And IVDD mice with Piog treatment significantly preserved the nucleus pulposus area, with a notable increase compared to the IVDD model without pioglitazone treatment (Fig. [184]10L). Immunofluorescence for COL2A1 showed significantly higher intensity in the Piog + IVDD group compared to that in the IVDD group, suggesting increased extracellular matrix production after pioglitazone treatment (Fig. [185]10J, M). MMP13 immunofluorescence intensity was significantly lower in the Piog-treated IVDD group compared to the IVDD group without pioglitazone treatment, supporting reduced catabolic activity (Fig. [186]10J, N). Fig. 10. [187]Fig. 10 [188]Open in a new tab Pioglitazone and Annexin V modulate macrophage efferocytosis in IVDD progression. (A-D) mRNA expression analysis showing the relative mRNA expression levels of MERTK, AXL, COL2A1, and MMP13 in control (Ctrl), pioglitazone (Piog), IL-1β, and IL-1β + Piog groups (n = 3; one-way ANOVA). (E-I) Representative Western blot images and quantification of MERTK, AXL, COL2A1, and MMP13 protein levels normalized to β-actin or GAPDH as loading controls in different treatment groups (n = 3; one-way ANOVA). (J-N) Representative images of HE, SOFG staining, and immunofluorescence staining for COL2A1 and MMP13 (n = 6). (K) Histological score of intervertebral discs (n = 6; one-way ANOVA). (L) Nucleus pulposus area size in intervertebral discs of various groups (n = 6; one-way ANOVA). (M-N) Relative IF intensity of COL2A1 and MMP13 in each group (n = 6; one-way ANOVA). (O-R) Quantitative real-time PCR analysis showing the relative mRNA expression levels of MERTK, AXL, COL2A1, and MMP13 in control (Ctrl), anakinra (Ann), IL-1β, and IL-1β + Ann groups (n = 3; one-way ANOVA). (S-U) Representative Western blot results and quantification of (T) MERTK and AXL protein levels, (U) COL2A1 and MMP13 protein levels normalized to β-actin or GAPDH (n = 3; one-way ANOVA). (V) Representative images of HE, SOFG staining, and immunofluorescence staining for COL2A1 and MMP13 (n = 6; one-way ANOVA). (W) Histological score of intervertebral discs (n = 6; one-way ANOVA). (X) Nucleus pulposus area size of intervertebral discs (n = 6; one-way ANOVA). (Y-Z) Relative IF intensity of COL2A1 and MMP13 of discs in each group (n = 6; one-way ANOVA). Data represent mean ± SD. *p < 0.05, **p < 0.01, ***p < 0.001, ****p < 0.0001, ns = not significant The mRNA levels of both MERTK and AXL in the IL-1β + Ann group were significantly lowered compared to that in the IL-1β group (Fig. [189]10O-P). In addition, the IL-1β + Ann group showed significantly lower protein levels of both MERTK and AXL than the IL-1β group (Fig. [190]10S-T), indicating that Annexin V suppressed the efferocytosis. Ann significantly suppressed COL2A1 mRNA expression and raised MMP13 mRNA level under IL-1β treatment (Fig. [191]10Q-R); and Ann lowered COL2A1 protein level and increased MMP13 protein expression under IL-1β treatment (Fig. [192]10S, U). The above results imply that Annexin V could promote ECM degradation through suppressing efferocytosis. HE and SOFG staining indicated that the damage of disc structure, induced by IVDD, was further promoted by Ann treatment (Fig. [193]10V). The histological score in the IVDD + Ann group was significantly higher than that in the IVDD group (Fig. [194]10W), and Ann treatment in IVDD model lowered nucleus pulposus area size (Fig. [195]10X). What’s more, the immunofluorescence for COL2A1 was lower in Ann-treated IVDD mice compared to IVDD mice, indicating lowered matrix production (Fig. [196]10Y); and MMP13 immunofluorescence intensity was raised in Ann-treated groups, supporting enhanced matrix degradation (Fig. [197]10Z). Discussion The primary aim of this study was to investigate the relationship between macrophage efferocytosis and IVDD and to identify the key regulators. Our major findings demonstrated that LGR6 plays a crucial role in regulating IVDD through the enhancement of macrophage efferocytosis. These findings not only deepen our understanding of the molecular mechanisms underlying IVDD but also suggest that targeting LGR6 could be a novel therapeutic approach for managing degenerative disc diseases. Efferocytosis, a process by which macrophages clear apoptotic cells, is a critical mechanism for maintaining tissue homeostasis [[198]22]. Recent research has highlighted the importance of efficient efferocytosis in various tissues, particularly in preventing chronic inflammation and subsequent tissue damage [[199]23]. Recent studies have shown that defects in efferocytosis are linked to the progression of several chronic inflammatory diseases, including periodontitis, intestinal inflammation, and osteoarthritis [[200]24–[201]26]. These conditions share common pathological features with IVDD, such as chronic inflammation and tissue degeneration, further underscoring the relevance of efferocytosis in degenerative processes. Our study reveals that enhancing efferocytosis through LGR6 overexpression significantly reduces the progression of IVDD. This finding is consistent with emerging research that suggests enhancing efferocytosis may be a viable therapeutic strategy for treating degenerative diseases. For instance, recent advances in nanotechnology have led to the development of nanoparticles that can enhance efferocytosis [[202]27], offering new therapeutic avenues for conditions characterized by impaired clearance of apoptotic cells. Additionally, recent work has shown that manipulating efferocytosis pathways, such as the MER receptor tyrosine kinase pathway, can reduce tissue damage and promote regeneration in models of tissue injury and degeneration [[203]28]. This aligns with previous research indicating that impaired efferocytosis leads to unresolved inflammation and contributes to the pathogenesis of various degenerative diseases [[204]29]. Moreover, studies on the role of efferocytosis in neurodegenerative diseases like Alzheimer’s disease have shown that enhancing this process can mitigate neuroinflammation and slow disease progression [[205]30], providing further evidence of the protective effects of efferocytosis across different tissues. In our present study, the findings shows that by promoting efferocytosis, LGR6 helps to clear apoptotic debris, thereby reducing the inflammatory burden and slowing down the degenerative process. This emerging understanding of efferocytosis as a critical regulatory mechanism in tissue homeostasis highlights its potential as a therapeutic target in a wide range of degenerative conditions, including IVDD. Our study provides novel insights into the molecular mechanisms by which LGR6 regulates efferocytosis. We observed that LGR6 overexpression upregulates the expression of key efferocytosis receptors, including MERTK, AXL, TYRO3, and CX3CR1, in macrophages. This upregulation enhances the ability of macrophages to engulf and clear apoptotic NPCs [[206]31], thereby reducing the pro-inflammatory signals that drive disc degeneration. Recent research has further elucidated the intricate signaling pathways involved in efferocytosis, highlighting the role of these receptors in not only facilitating the clearance of apoptotic cells but also in modulating immune responses. For instance, MERTK and AXL have been identified as crucial mediators of the anti-inflammatory phenotype in macrophages, promoting tissue repair and resolution of inflammation in various contexts, including subarachnoid hemorrhage and cancer [[207]32, [208]33]. Recent studies have also revealed that the signaling pathways downstream of these receptors, such as the PI3K-Akt and MAPK pathways, play essential roles in the regulation of macrophage polarization and the efferocytosis process [[209]28, [210]34]. These pathways help to balance pro-inflammatory and anti-inflammatory signals, ensuring that macrophages efficiently clear apoptotic cells without triggering excessive inflammation. In the context of IVDD, the upregulation of these pathways through LGR6 could provide a dual benefit: enhancing efferocytosis while simultaneously suppressing inflammatory responses that contribute to disc degeneration. Conversely, LGR6 knockdown resulted in decreased expression of these receptors, leading to impaired efferocytosis and accelerated degeneration. This is consistent with recent findings in other models of chronic inflammation and tissue degeneration, where loss of efferocytosis receptors like MERTK has been linked to worsened disease outcomes. For example, in models of tissue fibrosis, reduced expression of MERTK has been associated with increased inflammation, fibrosis, and tissue damage [[211]35], underscoring the importance of these receptors in maintaining tissue homeostasis. These findings suggest that LGR6 exerts its protective effects in IVDD by modulating the expression of efferocytosis-related genes, thereby enhancing the clearance of apoptotic cells and maintaining tissue integrity. The involvement of LGR6 in regulating these key pathways adds a new layer of complexity to our understanding of its role in tissue homeostasis and offers potential therapeutic avenues. Targeting LGR6 or its downstream signaling pathways could be a promising strategy to enhance efferocytosis and mitigate the progression of degenerative diseases such as IVDD. LGR6 has been implicated in various degenerative diseases due to its role in tissue regeneration and repair [[212]36]. Recent research has expanded our understanding of LGR6, particularly its involvement in stem cell maintenance and wound healing [[213]37, [214]38]. In the context of skeletal health, studies have shown that LGR6 promotes the repair of bone fractures and enhances the regeneration of osteochondral tissues [[215]36, [216]39], further solidifying its role in tissue repair mechanisms. In our study, LGR6 not only modulated efferocytosis but also played a direct role in regulating ECM synthesis and apoptosis in NPCs. The upregulation of LGR6 has been linked to the activation of Wnt/β-catenin signaling pathway [[217]40], which might lead to increased synthesis of ECM components and the inhibition of catabolic processes [[218]41]. We found that LGR6 overexpression led to increased expression of ECM components such as COL2A1 and reduced expression of matrix-degrading enzymes like MMP13, indicating a protective effect against matrix degradation. Additionally, LGR6 has been reported to interact with other signaling molecules like TGF-β [[219]42], which is known to play a pivotal role in ECM production and fibrosis. This interaction suggests that LGR6 may exert its protective effects through multiple pathways that converge on ECM regulation. Furthermore, LGR6 overexpression inhibited IL-1β-induced apoptosis in NPCs by upregulating anti-apoptotic proteins like BCL2 and downregulating pro-apoptotic markers such as cleaved caspase 3 and BAX. These findings are consistent with recent studies that have highlighted the regenerative role of LGR6 in other tissues, such as the skin, where LGR6 has been shown to protect against apoptosis and promote cell survival under stress conditions [[220]43, [221]44]. The anti-apoptotic effects of LGR6 are likely mediated through its influence on key signaling pathways, including the PI3K/Akt pathway [[222]45], which is known to regulate cell survival and apoptosis in response to inflammatory stimuli. However, our study provides the first evidence of LGR6’s role in modulating efferocytosis, adding a new dimension to its protective mechanisms in degenerative diseases. This novel finding opens up new research avenues, suggesting that LGR6 could serve as a dual-function therapeutic target that not only enhances tissue regeneration and ECM protection but also modulates immune responses to mitigate inflammation and degeneration. Given the protective role of LGR6 in IVDD through its modulation of macrophage efferocytosis, ECM homeostasis, and apoptosis inhibition, LGR6-based therapeutic strategies hold significant potential for IVDD treatment. One promising approach could be the development of small-molecule agonists or biologics that enhance LGR6 signaling to stimulate efferocytosis and maintain tissue integrity. This strategy could mitigate chronic inflammation and prevent progressive disc degeneration by promoting the efficient clearance of apoptotic cells. Another potential avenue is gene therapy-based strategies to enhance LGR6 expression within intervertebral disc tissues. Recent advances in viral vector-based gene delivery, such as adeno-associated viruses (AAVs), have shown promise in achieving localized and sustained gene expression in degenerative musculoskeletal tissues [[223]46, [224]47]. The use of AAV-LGR6 vectors might provide a targeted approach to upregulate LGR6 expression in degenerating discs, potentially restoring macrophage-mediated tissue repair mechanisms. Furthermore, pharmaceutical agents that indirectly modulate LGR6 activity by regulating upstream or downstream signaling pathways could offer another therapeutic strategy. A recent study by Cohen et al. revealed that the activation of the Wnt/β-catenin signaling pathway, which is known to be influenced by LGR6 [[225]48, [226]49], may enhance ECM synthesis and reduce catabolic enzyme activity [[227]50]. Small molecules targeting Wnt/β-catenin pathway activators could be explored as potential therapeutics in conjunction with LGR6-targeted approaches. Additionally, biomaterial-based approaches may provide a complementary strategy for LGR6-based therapies. Injectable hydrogels loaded with bioactive compounds that enhance LGR6 activity could be developed to provide sustained therapeutic effects within the degenerating disc microenvironment, which might help minimize systemic side effects while ensuring prolonged LGR6 activation for tissue repair. Despite the significant findings, our study has several limitations. First, while we demonstrated the role of LGR6 in modulating efferocytosis and IVDD in vitro and in animal models, the translation of these findings to human subjects requires further investigation. Second, the study primarily focused on the effects of LGR6 at the molecular and cellular levels, leaving the broader systemic implications less explored. Additionally, while we identified key efferocytosis-related receptors influenced by LGR6, the downstream signaling pathways involved remain to be fully elucidated. Future studies should aim to address these limitations by exploring the clinical relevance of LGR6 in human IVDD and investigating the detailed signaling mechanisms involved. Conclusions In conclusion, our study provides compelling evidence that LGR6 modulates IVDD through its regulation of macrophage efferocytosis. By enhancing the clearance of apoptotic cells and protecting against ECM degradation and apoptosis, LGR6 emerges as a potential therapeutic target for IVDD. These findings not only advance the understanding of the molecular mechanisms driving IVDD but also open new avenues for the development of targeted therapies aimed at mitigating degenerative disc diseases through the modulation of efferocytosis and related pathways. Further research is warranted to explore the translational potential of LGR6-based interventions in clinical settings. Electronic supplementary material Below is the link to the electronic supplementary material. [228]Supplementary Material 1^ (217KB, jpg) Acknowledgements