Abstract Background Low back pain represents a major global health issue, with intervertebral disc degeneration (IVDD) being one of its primary causes. Disc degeneration involves complex processes such as inflammation, matrix degradation, and cell death, yet the underlying mechanisms remain poorly understood. Single-cell RNA sequencing offers a powerful approach to elucidate cellular heterogeneity and dynamic changes in IVDD, providing valuable insights for early diagnosis and targeted therapeutic strategies. Methods The Harmony algorithm was used to integrate four independent single-cell sequencing datasets. Subtype identification, differential expression analysis, enrichment analysis, and cell proportion analysis were conducted to explore functional alterations in various nucleus pulposus cell (NPC) subpopulations. Molecular docking was employed to evaluate the stability of aspirin targeting GPX4. In vitro and in vivo experiments were performed to assess the therapeutic effects of aspirin on IVDD. Results Eight distinct NPC subtypes were identified based on cellular heterogeneity and their associated marker genes. The CDKN1A⁺aNPC subtype increased progressively with disease severity, while the matrix-supporting ABI3BP⁺mNPC and SOD3⁺mNPC subtypes significantly decreased in advanced degeneration. Concurrently, there was an increase in ECM remodeling-related LTBP1⁺mNPCs. Within the CDKN1A⁺aNPC, GPX4 was notably downregulated, suggesting the activation of ferroptosis. Molecular docking results revealed a high affinity of aspirin for GPX4. Additionally, aspirin inhibited ferroptosis and ameliorated disc structural damage. Conclusion The increased proportion of CDKN1A⁺aNPC cells serves as an early warning feature for the progression of IVDD. Aspirin stabilizes the targeting of GPX4, thereby inhibiting ferroptosis and exerting therapeutic effects on IVDD. Keywords: low back pain, intervertebral disc degeneration, single-cell analysis, cell death, glutathione peroxidase 4 Introduction Low back pain (LBP) adversely impacts individuals across all age groups and socioeconomic strata, with a lifetime prevalence ranging from 11% to 84% worldwide.[36]^1–3 LBP significantly compromises quality of life, is the leading cause of disability, and results in substantial losses in societal productivity while directly escalating healthcare costs.[37]^1 Although the etiology of LBP is multifactorial and complex, intervertebral disc degeneration (IVDD) has been identified as its primary driver.[38]^4^,[39]^5 IVDD represents an age- and injury-related pathological process within the intervertebral disc (IVD), driven by a series of intricate molecular mechanisms that evolve over time and ultimately lead to severe clinical manifestations.[40]^6^,[41]^7 The bidirectional crosstalk between the IVD and the adjacent bone marrow critically shapes the microenvironment and determines the progression of pathology.[42]^8^,[43]^9 Loss of homeostatic balance within the disc microenvironment triggers IVDD, characterized by a catabolic, hypoxic milieu and cellular senescence, which in turn drives immunometabolic alterations.[44]^10^,[45]^11 Current treatments for IVDD include conservative management and surgical interventions tailored to alleviate patient symptoms. However, these approaches primarily provide symptomatic relief and fail to reverse IVDD or restore the mechanical functionality of the spine. The development of intervertebral disc degeneration is influenced by several factors, including genetic predisposition, aging, mechanical injury, and nutritional deficiencies. Dysfunction of nucleus pulposus cells (NPCs), which reside in the gel-like central region of the intervertebral disc, is recognized as a key initiating factor in IVDD.[46]^12^,[47]^13 Pathological changes associated with NPC dysfunction primarily include cellular senescence and apoptosis, progressive extracellular matrix (ECM) degradation, fibrosis of the annulus fibrosus (AF), and inflammatory responses. These processes are marked by elevated levels of pro-inflammatory cytokines, such as tumor necrosis factor (TNF)-α, interleukin (IL)-1α/β, IL-6, and IL-17, which drive matrix degradation, chemokine production, and alterations in cell phenotype.[48]^7 The resultant imbalance between catabolic and anabolic activities leads to disc degeneration, herniation, and nerve root pain. Recent studies suggest that conditions observed in type 2 diabetes mellitus (T2DM), such as nutrient deprivation and hyperglycemia within the intervertebral disc, impose significant cellular stress on disc cells. This stress activates death receptors and triggers apoptosis via endoplasmic reticulum and mitochondrial pathways.[49]^14–16 Elevated glucose levels further induce reactive oxygen species (ROS) production and mitochondrial damage in AF,[50]^17 NPCs, and notochordal cells. Oxidative stress, a key mediator of ferroptosis, has been implicated in various diseases. Ferroptosis is a form of oxidative damage characterized by iron-dependent lipid peroxidation, driven by intracellular iron or lipoxygenase-catalyzed oxidation of unsaturated fatty acids in cell membranes. Hallmarks of ferroptosis include glutathione (GSH) depletion and inactivation of GPX4, the central regulator of the glutathione antioxidant system.[51]^18^,[52]^19 However, the mechanistic interplay between oxidative stress and ferroptosis in IVDD remains poorly understood. Elucidating the roles of oxidative stress and regulated cell death in the catabolic and injury processes of IVDD may unveil novel therapeutic targets for the treatment of symptomatic disc diseases. This study aims to elucidate the critical mechanisms underlying IVDD and explore potential new therapeutic targets. By integrating single-cell sequencing data from four independent IVDD cohorts and analyzing the heterogeneity of NPCs, we have conducted an in-depth investigation into the dynamic changes of NPCs during the progression of IVDD. The study focuses on the potential role of ferroptosis in IVDD, proposing that targeting ferroptosis through the regulation of GPX4 could represent a novel therapeutic strategy. Our findings highlight the significance of cellular heterogeneity in IVDD, particularly in processes such as matrix remodeling and cell death. Furthermore, we suggest that ferroptosis may play a pivotal role in the progression of IVDD, offering new insights for future therapeutic approaches. Methods Pre-Processing of Single-Cell Data Single-cell transcriptomic data were obtained from previous studies.[53]^20–23 All samples underwent stringent quality control following standard protocols implemented in R (version 4.4.0). The raw count matrix was imported using the “Read10X” function from the Seurat package and converted into a sparse matrix format. Individual datasets were integrated into a single combined object using the “merge” function, and unique cell labels were generated using the “RenameCells” function. Potential doublets were identified and removed using the Scrublet algorithm. Cells with fewer than 100 detected genes or genes expressed in fewer than three cells were excluded during the quality filtering process. Gene expression normalization across cells was performed using the “LogNormalize” method with a scaling factor of 10,000. To identify the top 2000 most variable genes, we used the “FindVariableFeatures” function. To account for unwanted technical variation, including UMI and mitochondrial gene content, these factors were regressed out using the “ScaleData” function. Principal component analysis was conducted on the highly variable features to reduce dimensionality, and the first 30 principal components were selected for subsequent analyses. Batch effects were corrected using the Harmony algorithm. For further dimensionality reduction, UMAP was employed. Clustering was performed based on the similarity of edge weights between cells using the Louvain method, with different resolution parameters (ranging from 0.1 to 2) tested to determine the optimal clustering configuration. A resolution of 0.6 produced the most robust and distinct clusters. Cell clusters were annotated by identifying differentially expressed genes within each cluster using the “FindAllMarkers” function, applying a nonparametric Wilcoxon rank-sum test with Bonferroni correction. Annotation of cell types was based on known surface markers, published gene lists, and the cell taxonomy database ([54]https://ngdc.cncb.ac.cn/celltaxonomy/), ensuring an accurate and comprehensive classification. Finally, only nucleus pulposus cells were extracted for further analysis. Enrichment Analysis Enrichment analysis of differentially expressed genes was performed using the “GSEABase” “ClusterProfiler” packages”, and “org.Hs.eg.db” database. The analysis utilized data from the GO and KEGG databases. Pathways with a P-value < 0.05 were considered significantly enriched. Analysis of Transcription Factor Activity The transcription factor activity was assessed using the SCENIC approach. First, GENIE3 was applied to construct a regulatory network by identifying co-expression patterns between transcription factors and their potential target genes. This network was then enhanced by integrating motif relationships and ranking the regulatory potential of motifs binding to genes, which led to the formation of a “regulon.” A regulon is defined as a set of target genes that are regulated by transcription factors through direct binding to upstream motifs. Finally, the activity of each regulon was evaluated across all cells using AUCell, which assesses the enrichment of regulon activity in individual cells. Establishment of a Rat Tail Intervertebral Disc Degeneration Model by Needle Puncture and Intervention Administration Twelve male Sprague-Dawley rats, eight weeks old and weighing approximately 300 grams, were used to establish a caudal intervertebral disc degeneration model. Inclusion criteria required the absence of congenital spinal deformities, prior trauma, infection, metabolic disease, immunodeficiency, and abnormal physiological parameters. Exclusion criteria included: (1) fractures, deformities, or congenital defects of the caudal vertebrae or intervertebral discs; skin damage, infection, or scarring in the surgical area; (2) recent exposure (within two weeks) to surgery, drug administration, or radiation, or prior participation in experiments; (3) restricted mobility, abnormal gait, pain behavior, significant weight loss (>10%), or abnormal food and water intake; and (4) procedural accidents during modeling, such as nerve or vascular injury from needle insertion. Following a two-week acclimatization period, the rats were randomly assigned to three groups: control, needle puncture + aspirin, and needle puncture only. Under 4% isoflurane anesthesia (maintained at 2% via a face mask), rats were placed in a prone position and immobilized. The intervertebral disc spaces at C6/7, C7/8, and C8/9 were identified and marked. The skin was sterilized with povidone-iodine, and a 21-gauge hypodermic needle was inserted dorsally into the C7/8 and C8/9 discs to a depth of 5 mm, guided by a depth-limiting cap. The needle was rotated 360° for 30 seconds before withdrawal to induce degeneration. For intervention, a 31-gauge microsyringe was used to inject 2 µL of aspirin solution into the C7/8 disc and 2 µL of phosphate-buffered saline (PBS) into the C8/9 disc along the puncture track. A similar procedure was performed at the C6/7 disc, where 2 µL of PBS was injected as a control. Following the procedure, injection sites were examined for bleeding and re-sterilized. Isoflurane was discontinued, and animals were monitored until full recovery before being returned to standard housing with free access to food and water. Radiographic and MRI Evaluation of Rat Tail Vertebrae Four weeks after the induction of intervertebral disc degeneration and corresponding interventions, radiographic (X-ray) and magnetic resonance imaging (MRI) examinations of the rat tail vertebrae were performed. The rats were anesthetized with 4% isoflurane for induction and maintained with 2% isoflurane inhalation. They were positioned in the prone position and secured with adhesive tape to ensure the tail remained straight during imaging. Sequential X-ray and MRI scans were conducted. Preparation of Rat Intervertebral Disc Tissue Sections Following radiographic and MRI evaluations, rats were euthanized by CO₂ inhalation, and tail vertebrae were harvested for histological analysis. The tail was transected at the upper vertebral body of the C6/7 intervertebral space. Skin was removed, and the tail was rinsed thoroughly with phosphate-buffered saline (PBS) before fixation in 4% paraformaldehyde for 48 hours. After fixation, residual paraformaldehyde was removed by repeated PBS washes, and the tail was decalcified in freshly prepared 10% EDTA-2Na solution, with the decalcification solution replaced every 2–3 days for a total of 8–12 weeks. Decalcification was considered complete when a 21G needle could penetrate the vertebral body without resistance. The intervertebral discs at C6/7, C7/8, and C8/9, along with adjacent vertebral bodies, were isolated and trimmed into tissue blocks. Each block was oriented with the transverse process trimmed to create a flat surface and placed into embedding cassettes with appropriate labeling. Tissue Processing and Embedding: Decalcified samples were dehydrated and embedded in paraffin using the following sequence: 70% ethanol (48 h) → 80% ethanol (2 h) → 90% ethanol (2 h) → 95% ethanol (2 h) → absolute ethanol I (2 h) → absolute ethanol II (2 h) → absolute ethanol III (1.5 h) → xylene I (1 h) → xylene II (1 h) → xylene III (0.5 h) → molten paraffin I (3 h) → molten paraffin II (3 h). Tissue blocks were oriented with the trimmed surface at the mold base, embedded in paraffin, and cooled on a refrigeration unit until fully solidified. Sectioning: Paraffin blocks were sectioned at a thickness of 5 µm using an ultramicrotome. Sections were floated on 40°C water, stretched, and mounted onto slides. Slides were dried at 60°C for 2 hours to complete preparation. Hematoxylin and Eosin (H&E) Staining Paraffin-embedded tissue sections (5 µm) were deparaffinized in xylene, rehydrated through a graded ethanol series, and rinsed in distilled water. Sections were stained with hematoxylin for 3–5 minutes, differentiated in 1% acid alcohol, blued in tap water, and counterstained with eosin for 1–3 minutes. After dehydration through graded ethanol and clearing in xylene, coverslips were mounted with a resin-based medium. Stained sections were examined under a light microscope. Safranin O–Fast Green Staining Paraffin-embedded tissue sections (5 µm) were deparaffinized in xylene and rehydrated through a graded ethanol series. Sections were stained with Fast Green solution (0.1%) for 3 minutes and rinsed briefly in 1% acetic acid to remove excess stain. Safranin O solution (0.1%) was then applied for 5 minutes. After rinsing in distilled water, sections were dehydrated through graded ethanol, cleared in xylene, and mounted with a resin-based medium. Stained sections were observed under a light microscope. Cell Culture Human intervertebral disc nucleus pulposus cells were purchased from Procell (#CP-H097, Wuhan, China) and passed to the third generation for experimental purposes. Cells were cultured in DMEM/F12 medium (#C11330500BT, Gibco) supplemented with 10% FBS (#10099141C, Gibco) and 100 μg/mL penicillin/streptomycin (#C0222, Beyotime). Cells were maintained in a humidified incubator at 37°C with 5% CO[2]. The cell model of Intervertebral disc degeneration (IVDD) was established using LPS (#L-2880, Sigma-Aldrich) according to the scheme proposed by Fan et al.[55]^24 CCK-8 Cells were plated onto a 96-well plate at a density of 1×10^4 cells per well, with a medium volume of 100 μL. After 48 hours of treatment with LPS, 10μL CCK-8 solution was added to each well and incubated at 37°C for 1 hour. Absorbance was measured at 450 nm using a VersaMax microplate reader (TECAN M1000 Pro). Flow Cytometry Assay For cell death assay, cells were stained with 1 μM Propidium Iodide (#ST511, Beyotime) and analyzed by flow cytometry after incubation. For ROS assay, cells were incubated with 1 μM DCFH-DA (#S0033S, Beyotime) at 37°C for 30 min, followed by flow cytometry analysis. Immunofluorescent Staining Cells were anchored to the coverslip and washed thrice in PBS before being fixed in 4% paraformaldehyde for 15 min, than rinsed in PBS, stained with BODIPY 493/503 (#C2053S, Beyotime) for 30 min at 37°C, and imaged using confocal microscopy (OLYMPUS FV3000). Western Blot Proteins were extracted using RIPA lysis buffer (#P0013B, Beyotime). Total protein level quantified using the BCA method (#P0009, Beyotime), the denatured protein was separated by SDS-PAGE, transferred to PVDF membranes, blocked with 5% fat-free milk, and probed with primary antibodies overnight at 4°C. Membranes were then incubated with HRP-conjugated secondary antibodies and visualized using chemiluminescent reagents. The following primary antibodies were used in this study: GPX4 Rabbit mAb (#59735), GAPDH Rabbit mAb (#2118), Anti-rabbit IgG, HRP-linked Antibody (#7074), all purchased from CST. Statistical Analysis Statistical analysis was performed using GraphPad Prism software. A t-test was applied to analyze group differences for normally distributed variables, with the Mann–Whitney U-test used for non-normally distributed variables. A p-value of less than 0.05 was considered statistically significant. Results Multi-Cohort Single-Cell Analysis Reveals Functional Remodeling of Nucleus Pulposus Cells in Intervertebral Disc Degeneration We integrated single-cell sequencing data from NPCs samples of intervertebral discs with varying degrees of degeneration across four independent cohorts.[56]^20^,[57]^21^,[58]^23 A total of 24 samples passed quality control, comprising 91,131 NPCs. Based on the Pfirrmann grading system, samples were categorized into mild degenerative disc (MDD, grades I–II, n=9) and severe degenerative disc (SDD, grades III–V, n=15) groups. Differential gene expression analysis between SDD and MDD NPCs identified 521 differentially expressed genes (DEGs) (|log[2]FC| > 0.25, adj. p-value < 0.05), including 284 upregulated and 237 downregulated genes. GO enrichment analysis revealed significant enrichment of these DEGs in pathways related to apoptosis, vascular development, skeletal system development, and collagen fibril organization ([59]Figure 1A). Consistent with the activation of angiogenesis during disc degeneration, our results further confirmed the vascularization processes in SDD. Furthermore, disruptions in collagen and fiber metabolism indicated ongoing ECM remodeling. KEGG pathway analysis showed significant enrichment of DEGs in pathways such as protein processing in the endoplasmic reticulum, cytoskeleton organization, ferroptosis, focal adhesion, and ECM-receptor interactions ([60]Figure 1B). Focal adhesion, a primary interface for cell-ECM interactions, involves the interplay between cytoskeletal structures and ECM receptors. Our findings suggest that NPCs in the SDD group exhibit matrix-fibrotic characteristics. Additionally, the protein processing pathway in the endoplasmic reticulum may regulate ferroptosis-related protein degradation, accelerating ferroptosis.[61]^25 Gene set enrichment analysis further demonstrated the activation of matrix-fibrotic features in NPCs of SDD, including processes related to cell adhesion, collagen fibril organization, and cell junction assembly ([62]Figure 1C). In contrast, GSEA revealed significant suppression of cellular oxidant detoxification in SDD, potentially impairing ROS-induced stress responses and reducing p53-mediated ferroptosis inhibition,[63]^26 thereby promoting ferroptosis. Our study reveals substantial remodeling of NPC characteristics during disc degeneration, involving key biological processes such as matrix reconstruction, cell adhesion, and ferroptosis. Figure 1. [64]Figure 1 [65]Open in a new tab (A) GO enrichment analysis of differentially expressed genes. (B) KEGG pathway enrichment analysis of differentially expressed genes. (C) GSEA enrichment analysis. Single-Cell Atlas of NPCs via Cellular Heterogeneity Analysis To further explore the heterogeneity of NPCs, we performed cellular subclustering, identifying three major types: matrix-NPCs (mNPC), apoptosis-NPCs (aNPC), and inflammation-NPCs (iNPC) ([66]Figure 2A). Based on specific cell marker genes, these NPCs were further categorized into eight subtypes: COL9A3+mNPC, SOD3+mNPC, COL14A+mNPC, MMP14+mNPC, CDKN1A+aNPC, ABI3BP+mNPC, LTBP1+mNPC, and SLC7A2+iNPC ([67]Figure 2A). UMAP analysis revealed the distribution patterns of these subtypes under different degenerative conditions ([68]Figure 2B), with differential gene expression identified for each subtype ([69]Figure 2C). We observed that ABI3BP+mNPCs exhibited upregulation of COL2A1 expression, accompanied by downregulation of COL1A2 and COL3A1. COL2A1, encoding type II collagen, is a key structural component of cartilage, with degradation mediated by matrix metalloproteinases (MMPs), contributing to cartilage degeneration. In contrast, COL1A2 and COL3A1 regulate collagen degradation via interaction with MMPs (eg, MMP-1, MMP-2). This suggests that ABI3BP+mNPCs may help maintain ECM balance and homeostasis in the early stages, thus protecting the nucleus pulposus structure and disc ecology. We also found aberrant expression of multiple CXCL family proteins in mNPCs such as COL14A+mNPC and LTBP1+mNPC.These proteins play crucial roles in ECM remodeling during IVDD, primarily affecting tissue repair and regeneration through modulation of cell migration, inflammation, angiogenesis, and fibrosis. This indicates that COL14A+mNPC and LTBP1+mNPC may regulate ECM remodeling during IVDD progression. Furthermore, we observed significant downregulation of GPX4 expression in late-stage apoptosis-type NPCs (CDKN1A+aNPC). GPX4, a key antioxidant enzyme, regulates cellular redox homeostasis to prevent ferroptosis.[70]^27^,[71]^28 The reduced expression of GPX4 suggests the activation of ferroptosis in CDKN1A+aNPCs.In summary, our study highlights the heterogeneous characteristics of NPCs in IVDD samples, including ECM remodeling and the activation of ferroptosis. These findings provide new insights into the complex matrix degeneration and apoptotic responses in nucleus pulposus tissue, offering a novel perspective on the pathogenesis of IVDD. Figure 2. [72]Figure 2 [73]Open in a new tab (A) Gene expression bubble plot for cell annotation. (B) Dimensionality reduction landscape of different cell subgroups at various disease stages. (C) Differentially expressed genes in each cell subgroup. Dynamic Shifts in Cell Subtypes Reflect Intervertebral Disc Degeneration Progression We analyzed the proportional changes of different NPC subtypes across various stages of IVDD ([74]Figure 3A). We found that the ferroptosis-associated subtype CDKN1A+aNPC increased progressively with disease progression, while matrix-supporting cells like ABI3BP+mNPC and SOD3+mNPC significantly decreased in the late degeneration stages. Conversely, the proportion of ECM remodeling-related cells, such as LTBP1+mNPC, increased, suggesting that dynamic changes in NPC proportions are key features of IVDD. Cell preference analysis (Ro/e)[75]^29 revealed that CDKN1A+aNPC and LTBP1+mNPC predominated in stages IV–V of degeneration, while ABI3BP+mNPC and SOD3+mNPC were concentrated in stages I–III ([76]Figure 3B), consistent with the observed proportional trends. Functional enrichment analysis showed that ABI3BP+mNPC and SOD3+mNPC were enriched in growth factor and skeletal development pathways, highlighting their role in maintaining disc homeostasis. CDKN1A+aNPC was enriched in apoptosis and programmed cell death pathways, potentially driving degenerative progression by disrupting tissue structure. Matrix-related subtypes (eg, COL14A+mNPC and MMP14+mNPC) were primarily involved in ECM synthesis and degradation, emphasizing the importance of matrix remodeling ([77]Figure 3C). Further analysis of transcription factor activity revealed that KLF12 and ELK1 were downregulated in early degeneration-associated cells (eg, ABI3BP+mNPC and SLC7A2+iNPC), suggesting suppressed matrix function and stress responses. In contrast, NR2F2 was downregulated in late degeneration-associated cells, indicating enhanced inflammatory and apoptotic responses. Additionally, the angiogenesis factor RBPJ was activated in LTBP1+mNPC but inhibited in SOD3+mNPC, revealing dynamic changes in angiogenesis during degeneration ([78]Figure 3D). Figure 3. [79]Figure 3 [80]Open in a new tab (A) Analysis of cell proportions in different disease stages for each cell subgroup. (B) Ro/e analysis of cell proportion preferences. (C)