We downloaded two bulk RNA datasets (GSE34095 and GSE70362) from the GEO database (https:// www.ncbi.nlm.nih.gov/geo) via the GEOquery package (Davis and Meltzer, 2007). The GSE34095 dataset is from Homo sapiens, and the data platform is GPL96, which contains six samples, including three control and three disease samples (He et al., 2015). The GSE70362 dataset is from H. sapiens, and the data platform is GPL17810, containing 48 samples, including 16 control and 32 disease samples (Kazezian et al., 2015). The two datasets were standardized by limma package (Ritchie et al., 2015), and the entire samples of both datasets were de-batched using the sva package (Leek et al., 2012) after integrating the two datasets. The integrated data set consisted of 54 samples, 19 control and 35 IDD samples, and were included in this study.

We downloaded single cell dataset GSE165722 from the GEO database, which contains eight samples, namely, S1 and S2 (Grade II), S3 and S4 (Grade III), S5 and S6 (Grade IV), and S7 and S8 (Grade V). Using the Seurat package (Hao et al., 2021), the expression matrix of the GSE165722 dataset was created as a Seurat object. We filtered the cells with >20% mitochondrial gene content as well as cells with features <200 or >4,000. We normalized the dataset’s sequencing depth using “SCTransform” to remove mitochondrial and cell cycle effects, identifying 3,000 highly variable genes. Principal Component Analysis (PCA) (David and Jacobs, 2014) was applied to identify significant principal components and visualize p-value distribution using the Elbowplot function. 15 principal components (PCs) were selected for t-SNE analysis to reduce dimensionality. We constructed K-nearest neighbors based on distances in PCA space using default parameters from the “FindNeighbors” function and 15 PC dimensions. Calling the “FindClusters” function divided cells into 17 clusters with a resolution of 0.5. Additionally, the “RunTSNE” function enabled dimensionality reduction for dataset visualization and exploration in cell annotation. Finally, AddModuleScore calculated scores for each cellular senescence.

Differential analysis of genes in different groups was performed using the limma package. p < 0.05 was set as the threshold for DEGs. Where log₂FC > 0 and log₂FC < 0 mean they are upregulated or downregulated in the disease, respectively. The results of the difference analysis were presented by plotting heat maps with the pheatmap package (Kolde and Kolde, 2015) and volcano maps with ggplot2 (Wickham, 2009).

The senescence-related dataset REACTOME CELLULAR SENESCENCE v2023.1 was obtained from the MSigDB database (Liberzon et al., 2015) containing 198 genes and intersected with the DEGs obtained by differential analysis for subsequent analytical studies.

Least absolute shrinkage and selection operator (LASSO) regression is a machine learning algorithm commonly used to build models today, using regularization to address the occurrence of overfitting during curve fitting and improve the accuracy of the model. We used the glmnet package (Simon et al., 2011) with integrated GSE34095 and GSE70362 for model construction with parameters set.seed 1) and family = “binomial” to further screen for variables associated with IDD based on key senescence-related DEGs. Finally, the risk score equation was calculated by optimizing the expression values of the genes and the associated regression coefficients: Riskscore = ∑ i = 1 n c o e f i * ⁡ exp i

The genes obtained by the LASSO algorithm were defined as key senescence-related DEGs.

The integrated GSE34095 and GSE70362 disease samples were divided into a high-risk score group and a low-risk score group. Differential analysis of genes in different groups was performed using limma. p < 0.05 was set as the threshold for DEGs. Gene Ontology (GO) analysis is a common approach for large-scale functional enrichment studies, including biological process (BP), molecular function (MF) and cellular component (CC). The Kyoto Encyclopedia of Genes and Genomes (KEGG) is a widely used database that stores information about genomes, biological pathways, diseases, and drugs. GO analysis and KEGG enrichment analysis of DEGs were performed using the clusterProfiler package (Yu et al., 2012), and results with p < 0.05 were considered statistically significant.

To investigate the differences in biological processes between different subgroups, based on the integrated GSE34095 and GSE70362 disease samples, we used gene set enrichment analysis (GSEA). GSEA is an algorithm used to evaluate the trend of distribution of genes of a predefined gene set in a table of genes ordered by their phenotypic correlation and thus determine their contribution to the phenotype (Subramanian et al., 2005). Gene set variation analysis (GSVA) is a non-parametric unsupervised analysis method that evaluates the enrichment of gene sets in the microarray nuclear transcriptome by converting the expression matrix of genes between samples into the expression matrix of gene sets between samples (Hänzelmann et al., 2013). This is used to assess whether different pathways are enriched across samples. We downloaded the c2.all.v7.2.symbols.gmt gene set from the MSigDB database for GSEA and GSVA analysis, and results with p < 0.05 was considered significantly enriched.

We screened genes with area under curve (AUC) > 0.7 by ROC curves, and based on the results of multifactorial analysis, multiple predictors were integrated and assigned according to certain proportions to visualize the interrelationship between each variable on outcome prediction in the form of a graph. We used multifactorial logistic regression to predict the incidence of IDD based on these senescence-related genes and plotted the Nomogram.

We performed ssGSEA analysis (Supplementary Table S1) on IDD samples based on immune cell markers by GSVA package, estimated the composition and abundance of 28 immune cells, and compared the differences in immune cells between high- and low-risk groups for IDD. The Spearman algorithm was used to assess the infiltration abundance of significant immune cells for correlation analysis with the risk score.

Networkanalyst is an online visual analysis platform for gene expression analysis and meta-analysis (Zhou et al., 2019). We used Networkanalyst to analyze key senescence-related DEGs associated transcription factors (TFs) based on the JASPAR database. Key senescence-related DEGs associated miRNAs were analyzed by the multiMiR package (Ru et al., 2014). Based on the above results, the correlation network was plotted using Cytoscape software (Shannon et al., 2003).

Human degenerative nucleus pulposus tissues were obtained from 6 patients (4 females and 2 males; age 60.33 ± 4.76 years; Grade IV) with degenerative disc disease undergoing surgery. The control samples were taken from 6 patients (3 females and 3 males; age 30.17 ± 7.61 years; Grade Ⅱ) undergoing surgery due to scoliosis or thoracolumbar fracture. The degenerative grade of NP tissues was determined using the MRI-based Pfirrmann grading system (Pfirrmann et al., 2001). Informed consent was obtained from patients preoperatively. This study protocol was approved by the Ethics Committee of the Xi’an Daxing Hospital, affiliated with Yan’an University.

Total RNA was extracted using the RNA extraction solution (cat. no. G3013, Servicebio, China), and the RNA concentration and purity were measured using a Nanodrop 2,000 spectrophotometer (Thermo Scientific, United States). The RNA samples were reverse transcribed into cDNA using a reverse transcription kit (cat. no. G3337, Servicebio, China), and the cDNA was used as a template for amplifying the target genes. The reaction was performed with 40 amplification cycles using the following protocol: denaturation at 95°C for 30 s, annealing at 60°C for 30 s, and extension at 72°C for 60 s. Samples were analyzed in triplicate, and the mRNA expression levels were calculated using the 2^−ΔΔCT method, with GAPDH serving as the internal reference. The sequences of the primers are listed in Table 1.

HE staining (cat. no. G1003, Servicebio, China) and Safranin O staining (cat. no. G1053, Servicebio, China) were performed according to manufacturer’s instructions, and the results were recorded under a microscope (Nikon, Japan).

RNA probes were designed and obtained from Servicebio (Wuhan, China), see Table 2 for probe details. Initially, paraffin slices underwent dewaxing and dehydration. Next, the slides were immersed in a boiling retrieval solution, treated with proteinase K for digestion, and then incubated with prehybridization solution for 1 h at 37°C. Following removal of the prehybridization solution, the slides were exposed to probe hybridization solution in a humidity chamber overnight at 40°C. Subsequently, the slides were washed, signal probe hybridization solution was added, and nuclear counterstaining with DAPI was performed according to the manufacturer’s instructions. Finally, the slices were placed under a fluorescence microscope (Nikon, Japan) for observation and image capture.

After deparaffinization and rehydration of paraffin-embedded tissue sections, they were blocked with 3% BSA blocking solution (cat. no. GC305010, Servicebio, China) at room temperature for 30 min. The sections were then incubated overnight at 4°C with primary antibodies: MAP4K4 (Rabbit, cat. no. 55247-1-AP, Proteintech, China), CXCL8 (Rabbit, cat. no. 17038-1-AP, Proteintech, China), TFDP2 (Rabbit, cat. no. 11500-1-AP, Proteintech, China), TNIK (Mouse, cat. no. 67948-1-Ig, Proteintech, China), CDKN2C (Rabbit, cat. no. AF0620, Affinity Biosciences, China). After washing the sections, they were incubated with secondary antibodies (cat. no. GB23303, Servicebio, China) for 50 min. Finally, the sections were rinsed with water, counterstained with hematoxylin, and mounted with a coverslip. The results were examined under an optical microscope (Nikon, Japan).

All data calculations and statistical analyses were performed using R software (https://www.r-project.org/, version 4.0.2). For the comparison of two groups of continuous variables were analyzed by Mann-Whitney U test (Wilcoxon rank sum test). The cut-off for statistical significance was set at p < 0.05.

We merged the GSE34095 and GSE70362 datasets then removed the batch effect. The integrated dataset was used to identify DEGs using the limma package and intersection was taken with senescence-related genes to identify senescence-related DEGs associated with IDD. Key senescence-related DEGs were identified by LASSO regression and risk scores were established based on this. Next, IDD patients were divided into high- and low-risk groups to compare the differences in pathways, function, and immune cell infiltration between the two groups. Finally, in the single-cell dataset GSE165722, key senescence-related DEGs were scored to evaluate their expression in IDD single-cell samples. The flow chart is shown in Figure 1. The GEO datasets are shown in Supplementary Table S2.

We performed a preliminary merge of the GSE34095 and GSE70362 datasets, showing the expression levels of 54 samples with box plots (Figure 2A) and the expression distribution of 54 samples with PCA plots (Figure 2B). We removed batch effects for the combined dataset and showed the expression levels of 54 samples using box plots (Figure 2C) and the expression distribution of 54 samples using PCA plots (Figure 2D). The results showed good correction with no batch effect in the sample.

A total of 2,122 DEGs were obtained, including 975 upregulated and 1,147 downregulated genes, and all DEGs were displayed in a volcano map (Figure 3A). The 30 genes with the largest log₂FC and 30 with the smallest log₂FC were plotted in the heat map (Figure 3B). A total of 26 senescence-related DEGs were obtained from the intersection of senescence-related genes and DEGs (Figure 3C).

To further identify key senescence-related DEGs for IDD, we performed LASSO regression analysis based on 26 senescence-related DEGs (Figure 4A) and obtained 12 key senescence-related DEGs, including, CABIN1, CDKN2B, CDKN2C, CXCL8, H1-5, MAP4K4, MAPK11, MDM2, MINK1, TFDP2, TINF2, and TNIK (Figure 4B). Based on these key genes we developed a diagnostic model. Risk score = −0.747× CABIN1 + 0.887× CDKN2B+ 1.06× CDKN2C + 1.21× CXCL8–0.0238× H1-5 + MAP4K4×2.19–0.356× MAPK11 + 0.935× MDM2–0.148× MINK1–0.0520× TFDP2–0.567× TINF2–1.36× TNIK. This risk score could distinguish well between the control and disease groups (Figure 4C). Subsequently, used chromosome mapping, we presented the chromosomal localization of these 12 genes (Figure 4D).

The expression levels of the 12 key senescence-related DEGs in the IDD and control groups were demonstrated by box plots. Five of these genes were significantly highly expressed in IDD: CDKN2B (p < 0.05), CDKN2C (p < 0.01), CXCL8 (p < 0.01), MAP4K4 (p < 0.001), and MDM2 (p < 0.05). Six genes were significantly low expressed in IDD: H1-5 (p < 0.01), MAPK11 (p < 0.05), MINK1 (p < 0.01), TFDP2 (p < 0.01), TINF2 (p < 0.05), and TNIK (p < 0.05) (Figure 5A). A total of 25 relevant miRNAs and 45 TFs were obtained through the JASPAR database and multiMiR, and the network was mapped (Figure 5B).

To screen high-quality NP tissues for subsequent experiments, we collected control (grade II) and IDD (grade IV) NP tissues using the MRI-based Pfirrmann grading system. Grade II NP tissues exhibited high water content and appeared as high signal intensity in T2-weighted MRI, presenting a more gelatinous texture under macroscopic observation. In contrast, grade IV NP tissues showed low signal intensity in T2-weighted MRI and varying degrees of fibrotic changes when observed macroscopically (Figure 6A). HE staining revealed distinct characteristics of the tissues associated with increasing IDD grade. Grade II NP cells displayed a uniformly dark-stained extracellular matrix (ECM) with a few small vacuolated cells. Conversely, grade IV cells exhibited a uniformly light-stained ECM with heavily aggregated and large vacuolated cells (Figure 6A). Safranin O staining indicated reduced levels of proteoglycan and collagen fiber in the ECM of grade IV tissues (Figure 6A). Immunohistochemical staining was employed to detect the expression of P16 and P21. The results revealed a higher expression of P16 and P21 in the IDD group, suggesting increased cellular senescence (Figure 6A). In the qPCR results, the expression levels of CDKN2C (p < 0.01), MAP4K4 (p < 0.01), TFDP2 (p < 0.05), and TNIK (p < 0.01) in the IDD group were significantly higher compared to the control group. However, CXCL8 showed a significant decrease (p < 0.01) in the IDD group (Figure 6B). The FISH and Immunohistochemistry results were generally consistent with the qPCR (Figures 6C, D).

We analyzed the ROC curves of key senescence-related DEGs for identifying IDD and controls, and labeled the corresponding AUC. AUC of CABIN1 = 0.639 (Figure 7A), AUC of CDKN2B = 0.689 (Figure 7B), AUC of CDKN2C = 0.746 (Figure 7C), AUC of CXCL8 = 0.753 (Figure 7D), AUC of H1-5 = 0.740 (Figure 7E), AUC of MAP4K4 = 0.880 (Figure 7F), AUC of AUC of MAPK11 = 0.698 (Figure 7G), AUC of MDM2 = 0.701 (Figure 7H), AUC of MINK1 = 0.771 (Figure 7I), AUC of TFDP2 = 0.731 (Figure 7J), AUC of TINF2 = 0.680 (Figure 7K) and AUC of TNIK = 0.783 (Figure 7L).

The incidence of IDD was predicted by constructing models using genes with AUC >0.75 (CXCL8, MAP4K4, MINK1, and TNIK), performing multifactorial logistic regression and plotting nomogram (Figure 8A), and the model was evaluated. We plotted the ROC curve (AUC = 0.905) (Figure 8B) and the cxalibration plot of the prediction model, suggesting that the prediction model has good predictive power for IDD (Figure 8C).

We identified 2,368 DEGs in the high- and low-risk groups for IDD and performed GO and KEGG enrichment analyses. These DEGs mainly affect biological processes such as cellular response to chemical stress, extrinsic apoptotic signaling pathway, and activation of protein kinase activity (Supplementary Figure S1A), molecular functions such as protein heterodimerization activity, transcription coactivator activity and small GTPase binding (Supplementary Figure S1B) and cellular components such as proteasome complex, endopeptidase complex, and proteasome accessory complex (Supplementary Figure S1C). Pathways such as pathways of neurodegeneration−multiple diseases, amyotrophic lateral sclerosis, and Alzheimer’s disease (Supplementary Figure S1D). Detailed information about the significantly enriched GO and KEGG categories is provided in Supplementary Tables S3, S4.

We performed a GSEA analysis (Supplementary Table S5) to further analyze the differences between the high- and low-risk groups for IDD. Among them, KEGG SYSTEMIC LUPUS ERYTHEMATOSUS (P.adj = 2.4e-06) (Figure 9A), KEGG TERPENOID BACKBONE BIOSYNTHESIS (P.adj = 0.011) (Figure 9B) and KEGG RIBOSOME (Figure 9C) (P.adj = 0.028) were significantly enriched in the low-risk group. KEGG METABOLISM OF XENOBIOTICS BY CYTOCHROME P450 (P.adj = 1.1e-03) (Figure 9D), KEGG PROTEASOME (P.adj = 2.4e-04) (Figure 9E), KEGG STEROID HORMONE BIOSYNTHESIS (P.adj = 1.6e-04) (Figure 9F) were significantly enriched in the high-risk group. We also performed GSVA analysis and drew a heat map showing the 20 most relevant KEGG pathways (Figure 9G). The pathways KEGG LIMONENE AND PINENE DEGRADATION and KEGG TERPENOID BACKBONE BIOSYNTHESIS were significantly enriched in the low-risk group and the pathways KEGG MTOR SIGNALING PATHWAY and KEGG UBIQUITIN MEDIATED PROTEOLYSIS were significantly enriched in the high-risk group. Figure 9H shows the 20 most relevant GO functions, GO INTERLEUKIN 10 SECRETION was significantly enriched in the low-risk group, GO POSITIVE REGULATION OF ERAD PATHWAY and GO PROTEIN LOCALIZATION TO CHROMATIN were significantly enriched in the high-risk group.

We used ssGSEA to calculate the levels of 28 immune cell infiltrates in patients in the high- and low-risk groups for IDD, plotted a heat map for display (Figure 10A), and plotted a box plot comparing the differences in immune cell infiltration levels between patients in the high-risk and low-risk groups. Nine immune cells were significantly different, with eosinophils being significantly elevated in patients in the low-risk group (p < 0.05); Effector memory CD8 T, central memory CD4 T, myeloid-derived suppressor cells (MDSCs), natural killer, monocyte, type 1 T helper cell, plasmacytoid dendritic cell, and natural killer T cells were significantly elevated in patients in the high-risk group (p < 0.05) (Figure 10B). Subsequently, we calculated the correlation between the abundance of immune cell infiltration and the risk score. Eosinophils were negatively correlated with the risk score (cor = −0.34, p = 0.048) (Figure 10C); Effector memory CD8 T cells were positively correlated with risk score (cor = 0.51, p = 0.0021) (Figure 10D); Central memory CD4 T cells were positively correlated with risk score (cor = 0.52, p = 0.0016) (Figure 10E), MDSCs were positively correlated with risk score (cor = 0.30, p = 0.079) (Figure 10F); Natural killer cells were positively correlated with risk score (cor = 0.34, p = 0.048) (Figure 10G); Monocytes were positively correlated with risk score (cor = 0.38, p = 0.042) (Figure 10H), Type 1 T helper cells were positively correlated with risk score (cor = 0.43, p = 0.010) (Figure 10I); Plasmacytoid dendritic cells were positively correlated with risk score (cor = 0.30, p = 0.084) (Figure 10J); Natural killer T cells were positively correlated with risk score (cor = 0.34, p = 0.046) (Figure 10K).

We performed single-cell data analysis on eight IDD samples from GSE165722 and divided the cell population into 16 clusters using tSNE analysis (Figure 11A), and each cluster was annotated by markers of the cells and plotted in bubble plots for display (Figure 11B). The cell cluster was divided into B cells, nucleated erythrocytes, macrophages, neutrophils, myeloid cells, myeloid progenitor cells and T or NK cells according to the results of cell annotation (Figure 11C), and a proportional histogram was drawn to show the proportion of cells in these eight IDD samples (Figure 11D).

Senescence scores were calculated for individual cells (B cells, nucleated red blood cells, macrophages, neutrophils, myeloid cells, myeloid progenitor cells, and T or NK cells) based on key senescence-related DEGs (CABIN1, CDKN2B, CDKN2C, CXCL8, H1-5, MAP4K4, MAPK11, MDM2, MINK1, TFDP2, TINF2, and TNIK) using the AddModuleScore function (Figure 12A). Using the calculated senescence score to plot tSNE, we found that the senescence score was significantly highly expressed in neutrophils (Figure 12B). We show the senescence scores of these eight IDD samples using box plots (Figure 12C) and that IDD patients with IV/V classification have lower senescence scores relative to IDD patients with II/III classification. This result suggests that these genes may play a role in neutrophils in IDD (Figure 12D).