Comprehensive analysis of senescence-related genes and immune infiltration in intervertebral disc degeneration: a meta-data approach utilizing bulk and single-cell RNA sequencing data

1 Introduction

Lower back pain is one of the most common health problems, affecting up to 80% of people throughout their lives, causing severe disability worldwide, and resulting in significant medical costs (Qaseem et al., 2017). Intervertebral disc degeneration (IDD) is the leading cause of lower back pain (Cheng et al., 2018). Current treatments include medication, physical therapy, or surgical interventions that can partially relieve symptoms in the lower back and legs (Xin et al., 2022). They do not focus on replacing nucleus pulposus loss and restoring disc structure. This can lead to unsatisfactory results such as degeneration or recurrence of adjacent motion segments. Early diagnosis and timely treatment of IDD can effectively delay disease progression and significantly reduce the incidence of disability. Therefore, screening for diagnostic genes associated with IDD and elucidating their underlying pathogenesis could significantly prevent and treat IDD and may provide new avenues for the clinical treatment of IDD.

Senescence is a natural process, and there is a substantial link between senescence and IDD, which has been identified as a key factor in IDD (Cheng et al., 2022; Wu et al., 2022). Senescence leads to structural, metabolic, and biomechanical changes in the intervertebral disc that ultimately affect its function (Le Maitre et al., 2007). The senescent cells have several key features, such as sustained growth arrest, expression of anti-proliferative molecules (e.g., p16, INK4a) and activation of damage sensing signaling pathways (e.g., p38 MAPK and NF-kB) (He and Sharpless, 2017; Patil et al., 2019; Zhang et al., 2021). This in turn leads to the release of cytokines, chemokines, and other secreted phenotypic proteins associated with senescence, ultimately leading to inflammation and tissue degradation (Muñoz-Espín and Serrano, 2014). The senescence process is accompanied by changes in the immune microenvironment (Lanna et al., 2022). Studies have shown that immune dysregulation, including abnormal macrophage polarization, abnormal T-cell differentiation, abnormal expression of B cells, and improper T-cell differentiation are significantly associated with IDD (Song et al., 2022). Although both senescence and immune dysregulation are present in the IDD process, it is unclear if or how the two phenomena are linked. Therefore, identifying the molecular mechanisms of IDD senescence and further investigating the potential link between nucleus pulposus cell senescence and the immune microenvironment will provide new therapeutic concepts and alternative options for reversing the disease to provide better clinical treatment guidance.

Although previous study has provided preliminary evidence for the regulatory role of senescence-related genes in the development of IDD (Xu et al., 2023), they primarily focused on the biological functions of these genes in peripheral blood. The diagnostic value of these genes in nucleus pulposus tissue and their relationship with immune cells remain unclear. Thorough and comprehensive research on senescence-related genes in the nucleus pulposus tissue can enhance our understanding of IDD. Unlike previous publications that relied on experimental data generation, our research primarily focuses on meta-data analysis. We have extensively explored senescence-related biomarkers by publicly available datasets, including bulk and single-cell RNA sequencing data, through machine learning. In addition, we collected nucleus pulposus tissue samples for experimental verification of the analysis results to improve their accuracy and reliability. This approach allows us to assess a larger cohort and explore potential associations across various datasets, providing a broader perspective on IDD research.

2 Materials and methods

2.1 Bulk RNA data download

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.

2.2 Single cell data pre-processing and analysis

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.

2.3 Differential expression analysis

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).

2.4 Senescence-related dataset acquisition

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.

2.5 Model construction

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:

$\text{Riskscore}={\displaystyle \sum _{i=1}^n}{coef}_i*{\exp }_i$

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

2.6 Gene Ontology and Kyoto Encyclopedia of Genes and Genomes

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.

2.7 Gene set enrichment analysis and gene set variation analysis

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.

2.8 Nomogram construction

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.

2.9 Immune cell infiltration analysis

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.

2.10 MiRNA/TF-gene network construction

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).

2.11 Patient samples

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.

2.12 qPCR

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.

TABLE 1

TABLE 1. Primer sequences used for qPCR.

2.13 HE and Safranin O staining

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).

2.14 Fluorescence in situ hybridization (FISH)

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.

TABLE 2

TABLE 2. Probe sequences used for FISH.

2.15 Immunohistochemistry

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).

2.16 Statistical analysis

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.

3 Results

3.1 Flow chart of the study

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.

FIGURE 1

FIGURE 1. Flow chart of the study.

3.2 Datasets integration and batch effect removal

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.

FIGURE 2

FIGURE 2. Datasets integration and batch effect removal. (A) Box plots of gene expression levels for the combined two datasets (GSE34095 and GSE70362). (B) Principal component analysis (PCA) plots of the combined two datasets (GSE34095 and GSE70362). (C) Box plots of gene expression levels for the two datasets (GSE34095 and GSE70362) after removal of batch effects. (D) PCA plots of two datasets (GSE34095 and GSE70362) after removing batch effects.

3.3 Identification of senescence-related DEGs

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).

FIGURE 3

FIGURE 3. Identification of senescence-related differentially expressed genes (DEGs). (A) Volcano diagram of DEGs, red represents upregulation and blue represents downregulation. (B) Heat map of DEGs, red represents upregulation and blue represents downregulation. (C) Venn diagram of 26 senescence-related DEGs.

3.4 Lasso regression to establish senescence-related risk scores

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).

FIGURE 4

FIGURE 4. Lasso regression analysis. (A) Lasso regression analysis based on 26 senescence-related differentially expressed genes (DEGs). (B) Lasso regression analysis obtained 12 key senescence-related DEGs. (C) Box plot showing risk scores for the intervertebral disc degeneration (IDD) and control groups. Red represents the IDD group and blue represents the control group. (D) Chromosome map showing the chromosomal localization of these 12 key senescence-related DEGs.

3.5 TFs and miRNAs prediction of key senescence-related DEGs

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).

FIGURE 5

FIGURE 5. Transcription factors (TFs) and miRNAs prediction of key senescence-related differentially expressed genes (DEGs). (A) Box plots showing the expression levels of CABIN1, CDKN2B, CDKN2C, CXCL8, H1-5, MAP4K4, MAPK11, MDM2, MINK1, TFDP2, TINF2, and TNIK in the intervertebral disc degeneration (IDD) and control groups. Red represents the IDD group, blue represents the control group. *p < 0.05, **p < 0.01, ***p < 0.001. (B) Interaction network diagram showing 25 associated miRNAs and 45 TFs for 12 key senescence-related DEGs.

3.6 Validation of key senescence-related DEGs

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).

FIGURE 6

FIGURE 6. Validating key senescence-related DEGs. (A) Representative images were obtained using magnetic resonance imaging, and NP tissues obtained from patients with intervertebral disc degeneration were stained with HE, Safranin O, and immunohistochemical staining for P16 and P21. (B) qPCR validate the expression levels of CABIN1, CDKN2B, CDKN2C, CXCL8, H1-5, MAP4K4, MAPK11, MDM2, MINK1, TFDP2, TINF2, and TNIK in IDD and control nucleus pulposus tissue. Red represents the IDD group, blue represents the control group. *p < 0.05, **p < 0.01. Representative images of CDKN2C, MAP4K4, TFDP2, TNIK, and CXCL8 expression in NP tissue detected by FISH (C) and Immunohistochemistry (D).

3.7 ROC curve

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).

FIGURE 7

FIGURE 7. Receiver operating characteristic (ROC) curve. (A) ROC curve of CABIN1. (B) ROC curves of CDKN2B. (C) ROC curves of CDKN2C. (D) ROC curve of CXCL8. (E) ROC curves of H1-5. (F) ROC curve of MAP4K4. (G) ROC curve of MAPK11. (H) ROC curves of MDM2. (I) ROC curve of MINK1. (J) ROC curve of TFDP2. (K) ROC curve of TINF2. (L) ROC curves of TNIK.

3.8 Nomogram

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).

FIGURE 8

FIGURE 8. Clinical prediction model. (A) Nomogram plots based on CXCL8, MAP4K4, MINK1, and TNIK expression to identify intervertebral disc degeneration (IDD) and controls. (B) Receiver operating characteristic (ROC) curves of clinical prediction models. (C) Calibration plot of the clinical prediction model.

3.9 GO and KEGG enrichment analysis

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.

3.10 GSEA and GSVA enrichment analysis

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.

FIGURE 9

FIGURE 9. Gene set enrichment analysis (GSEA) and gene set variation analysis (GSVA). (A–C) KEGG SYSTEMIC LUPUS ERYTHEMATOSUS (A), KEGG TERPENOID BACKBONE BIOSYNTHESIS (B), KEGG RIBOSOME (C) were significantly enriched in the low-risk group. (D–F) KEGG CELL CYCLE (D), KEGG DNA REPLICATION (E), KEGG MISMATCH REPAIR (F) were significantly enriched in the high-risk group. (G) Heat map showing the 20 most significantly different Kyoto Encylopedia of Genes and Genomes (KEGG) pathways in the high- and low-risk groups. Red and blue represent pathways with high and low expression in that group, respectively. (H) Heat map showing the 20 most significantly different Gene Ontology (GO)-related functions in the high- and low-risk groups. Red represents pathways with high expression in that group, blue represents pathways with low expression in that group.

3.11 Immune infiltration analysis

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).

FIGURE 10

FIGURE 10. Risk score and immune cell infiltration. (A) Heat map showing the abundance of 28 immune cell infiltrates in intervertebral disc degeneration (IDD). (B) Box plot showing the difference in immune infiltration levels between high- and low-risk groups, with blue and red representing low- and high-risk groups, respectively. (C) Correlation analysis of eosinophil and risk scores. (D) Correlation analysis of effector memory CD8 T cells and risk scores. (E) Correlation analysis of central memory CD4 T cells and risk scores. (F) Correlation analysis between myeloid-derived suppressor cells (MDSCs) and risk score. (G) Correlation analysis between natural killer cells and risk score. (H) Correlation analysis of monocytes and risk score. (I) Correlation analysis between Type 1 T helper cells and risk score. (J) Correlation analysis between plasmacytoid dendritic cells and risk score. (K) Correlation analysis of natural killer T cells and risk score.

3.12 GSE165722 single-cell analysis

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).

FIGURE 11

FIGURE 11. GSE165722 single-cell sequencing analysis. (A) Reduced-dimensional cluster analysis divides cells into 16 clusters and presents them with t-distributed stochastic neighbor embedding (t-SNE) plots. (B) Cellular annotation heat map showing the expression levels of markers for 16 clusters. (C) Annotation of 16 clusters containing 7 cell types: B cells, nucleated erythrocytes, macrophages, neutrophils, myeloid cells, myeloid progenitor cells, and T or NK cells. (D) Histogram of cell proportions, where the horizontal coordinates represent cell proportions and the vertical coordinates represent sample names.

3.13 GSE165722-based senescence scoring

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).

FIGURE 12

FIGURE 12. Senescence score based GSE165722. (A) AddModuleScore calculated senescence score. Violin diagram of different cell expression. (B) Senescence score distributed stochastic neighbor embedding (t-SNE). Blue represents high scores, while white represents low scores. (C) Box plots showing the senescence scores of different samples in GSE165722. (D) Box plots showing senescence score for different grades of intervertebral disc degeneration (IDD), with red representing grade II/II and blue representing grade IV/V.

4 Discussion

IDD is the most common cause of lower back pain (Cheng et al., 2018). The pathophysiological mechanisms of IDD are complex and not fully understood. IDD is closely associated with senescence. Therefore, it is crucial to explore the key senescence-related genes and their pathogenic mechanisms to provide an experimental and theoretical basis for the repair of IDD.

In this study, 12 key senescence-related DEGs were identified. A risk score was established based on this, which could distinguish well between the control and IDD groups. In addition, a nomogram for predicting IDD prevalence was constructed by selecting genes with AUC >0.75 (CXCL8, MAP4K4, MINK1, and TNIK). The ROC curve and the calibration plot suggest that the nomogram has good predictive power for IDD. Calcineurin upregulation was observed in hypoxic nucleus pulposus cells (Huang et al., 2020). CABIN1, a natural calcineurin antagonist, is speculated to play a significant role in IDD. CDKN2B and CDKN2C regulate the cell cycle and are closely linked to cellular senescence (Gagrica et al., 2012; Park et al., 2023). These findings suggest that CDKN2B and CDKN2C could be potential targets for IDD. However, further experiments are needed to confirm the above hypothesis. CXCL8 plays a role in neurogenic pain caused by disc herniation, with its expression correlating to the severity of degenerative tissue changes in the nucleus pulposus (Phillips et al., 2013; Pedersen et al., 2015). MDM2 is highly expressed in degenerated intervertebral discs (Wang et al., 2021). Our study identified significant upregulation of CXCL8 and MDM2 expression in IDD, consistent with previous findings. Prior research demonstrates significant overexpression of MAPK11 in IDD (Yang et al., 2016). Conversely, our study suggests downregulation of MAPK11. Further clinical samples and cellular experiments are required to validate MAPK11 expression in IDD. No prior research has examined the role of H1-5 in IDD. Our study discovered low expression levels of H1-5, highlighting its potential as a promising molecular target for IDD. miR-547-3p regulates inflammatory cytokines by targeting MAP4K4, alleviating neuropathic pain from lumbar disc herniation (Yao et al., 2022). MINK1 and TNIK are associated with senescence, suggesting their potential role in IDD (Reizis, 2019; Yu et al., 2020). Our findings demonstrate downregulation of TFDP2 and TINF2 in IDD. However, their roles in IDD have not been reported and require further investigation.

To further investigate the biological function and pathway differences between patients at different risks, we performed functional and pathway enrichment analysis. The results of this study showed that the two groups of patients differed mainly in cellular response to chemical stress, extrinsic apoptotic signaling pathway, activation of protein kinase activity, and transcription coactivator activity. Oxidative stress has been confirmed as a significant contributor to IDD, with excessive reactive oxygen species causing cellular oxidation and chemical stress, leading to apoptosis of functional disc cells and hastening degeneration (Lin et al., 2020). Protein kinases are critical regulators of cellular function, governing multiple intracellular signaling pathways to maintain normal physiological cell activity (Kupka et al., 2020; Du et al., 2022). DNA-binding transcription factors interact with co-activator complexes to enhance gene expression (Dietrich et al., 2017). Mechanism-sensitive transcriptional coactivators MRTF-A and YAP/TAZ play a role in regulating the phenotype of nucleus pulposus cells through cell shape modulation (Fearing et al., 2019). IL-10 secretion was prominently enriched in the low-risk group, while the high-risk group exhibited significant enrichment in the mTOR pathway and positive regulation of ERAD pathways. Previous research has indicated that IL-10 delays IDD by inhibiting the p38 MAPK pathway, suggesting its protective role in IDD, consistent with our study findings (Ge et al., 2020). ERAD is the primary mechanism for clearing misfolded proteins from the endoplasmic reticulum and has potential as a target for IDD treatment (Wen et al., 2021). Recent studies demonstrate the involvement of mTOR signaling in regulating intervertebral disc cell functions like oxidative stress, inflammation, cellular senescence, and apoptosis (Chen et al., 2023). Further investigation is warranted to explore the potential of targeting mTOR signaling and promoting ERAD pathways as molecular interventions for IDD.

Although there have been some reports on IDD and the immune microenvironment (Saberi et al., 2021; Wang et al., 2021), the relationship between these immune cells and genes remains unknown. High eosinophil levels were found in the low-risk group, with a negative correlation between eosinophil infiltration and the risk score. Effector memory CD8 T, central memory CD4 T, MDSC, natural killer cells, monocyte, type 1 T helper, plasmacytoid dendritic, and natural killer T cells were significantly elevated in the high-risk group of patients, and their infiltration abundance was all positively correlated with the risk score. We discovered that eosinophil abundance is protective in IDD. However, additional research on eosinophils in IDD is required to confirm their role. Plasmacytoid dendritic cells (pDCs) are implicated in IFN-producing diseases (Li et al., 2019). However, the role of pDCs in IDD development has not been reported. Our study found significantly elevated pDC levels in the high-risk group, highlighting the need for further investigation on the pDC-IDD interaction. Studies have shown increased infiltration of CD8⁺ and CD4⁺ T cells in IDD, leading to exacerbation of the inflammatory response (Li et al., 2020; Lan et al., 2022). MDSCs were significantly enriched in severely degenerated nucleus pulposus tissue (Tu et al., 2022). In vitro assays demonstrated that NK cells and macrophages have immune function in the early stages of IDD (Stich et al., 2015). There is evidence of monocyte infiltration occurring during the early degenerative stages (Li et al., 2023). Research has shown that in IDD, activated T cells predominantly differentiate into Th1, Th2, and Th17 subsets, with Th1 regulating macrophage involvement in IDD (Francisco et al., 2022). The above studies suggest that patients with IDD of different degenerative periods and degrees have different patterns of immune cell infiltration, which is supported by our findings. These results further deepen the understanding that immune cell infiltration has on different functions in patients with different risks of IDD.

Senescent cells form a senescence associated secretory phenotype (SASP), which induces inflammation, recruits immune cells (such as macrophages, natural killer cells, and neutrophils), and spreads senescence to other cells, thus exerting long-range effects on surrounding cells and tissues (Kale et al., 2020). Through single cell sequencing analysis, we found that the senescence score was significantly high in neutrophils. Neutrophils are typically the first immune cells to arrive at the site of inflammation and respond to bacterial infections. IL-8, as a major SASP cytokine, attracts neutrophils and stimulates their release of antimicrobial granules (Prata et al., 2018). Lagnado et al. (2021) showed that neutrophils can trigger senescence in neighboring cells by delivering reactive oxygen species. Neutrophils show a strong correlation with cellular senescence, indicating potential benefits in countering neutrophil-induced senescence for aging and age-related diseases. This aligns with our findings. Further research is needed to explore the interaction between senescence-related genes and neutrophil infiltration in IDD.

In contrast to previous studies, our research encompasses a comprehensive exploration of mRNA, miRNA, and TFs, enabling a more in-depth understanding of the molecular mechanisms underlying IDD. This approach allows for an extensive examination of gene features and potential regulatory networks related to IDD, providing a unique contribution in comparison to studies that solely concentrate on mRNA. Furthermore, we employed PCR, in situ hybridization, and immunohistochemistry to experimentally validate diagnostic markers at both the transcriptional and protein levels. This rigorous approach enhances result reliability and provides valuable confirmation for our findings. Finally, single-cell sequencing analysis found that neutrophils had a very high senescence score, indicating that senescence-related genes may play a role in neutrophils in IDD, further deepening the mechanism research. This methodological difference allows for a more comprehensive examination of gene features and potential regulatory networks associated with IDD, offering new insights into the molecular mechanisms underlying the disease.

Through experimental verification and analysis, we observed discrepancies in the expression trends of certain genes compared to the preliminary analysis. These inconsistencies may be attributed to technical limitations, sample heterogeneity, disease complexity, and small sample size. Technical limitations arise from inherent errors and variability associated with qPCR, which can lead to different expression trends. Sample heterogeneity, including variations in age, gender, disease severity, and treatment status, can influence gene expression and contribute to divergent trends. The complexity of intervertebral disc degeneration (IDD) involves multiple regulatory mechanisms, resulting in varying gene expression among patients. Additionally, the limited sample size limits the accuracy and reliability of our findings, leading to unstable and inconsistent results.

There are some limitations in this study. First, although we have conducted validation using nucleus pulposus tissue, further biological experiments are needed to confirm the clinical value of this study due to the small sample size. In addition, the sample size of this study was limited due to the lack of high-quality datasets, and the analysis results may be biased. Therefore, in future studies, we will continue to expand our sample size and conduct basic experiments to validate our results for better clinical application.

In summary, this study has developed a diagnostic model that significantly associates immune cell infiltration with IDD based on key senescence-related DEGs. Early screening and effective prevention of high-risk populations from the perspective of senescence-related genes will have a profound impact on the management of IDD.

Data availability statement

The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding author.

Ethics statement

The studies involving humans were approved by the Ethics Committee of the Xi’an Daxing Hospital, affiliated with Yan’an University. The studies were conducted in accordance with the local legislation and institutional requirements. The participants provided their written informed consent to participate in this study.

Author contributions

Y-JD: Conceptualization, Writing–original draft. X-GW: Data curation, Writing–original draft. ZL: Formal Analysis, Writing–review and editing. BW: Investigation, Writing–original draft. JL: Methodology, Writing–review and editing. JM: Software, Writing–original draft. XX: Validation, Writing–original draft. XT: Formal Analysis, Writing–review and editing. Q-CL: Software, Writing–original draft. J-YL: Data curation, Software, Writing–original draft. YZ: Data curation, Writing–original draft. BY: Funding acquisition, Writing–review and editing.

Funding

The author(s) declare financial support was received for the research, authorship, and/or publication of this article. This work was supported by the National Science Foundation of China (No. 82260443).

Acknowledgments

We would like to thank Editage (www.editage.cn) for English language editing.

Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fmolb.2023.1296782/full#supplementary-material

References

Chen, H. W., Zhou, J. W., Zhang, G. Z., Luo, Z. B., Li, L., and Kang, X. W. (2023). Emerging role and therapeutic implication of mTOR signalling in intervertebral disc degeneration. Cell Prolif. 56, e13338. doi:10.1111/cpr.13338

PubMed Abstract | CrossRef Full Text | Google Scholar

Cheng, F., Yang, H., Cheng, Y., Liu, Y., Hai, Y., and Zhang, Y. (2022). The role of oxidative stress in intervertebral disc cellular senescence. Front. Endocrinol. (Lausanne) 13, 1038171. doi:10.3389/fendo.2022.1038171

PubMed Abstract | CrossRef Full Text | Google Scholar

Cheng, X., Zhang, L., Zhang, K., Zhang, G., Hu, Y., Sun, X., et al. (2018). Circular RNA VMA21 protects against intervertebral disc degeneration through targeting miR-200c and X linked inhibitor-of-apoptosis protein. Ann. Rheum. Dis. 77, 770–779. doi:10.1136/annrheumdis-2017-212056

PubMed Abstract | CrossRef Full Text | Google Scholar

David, C. C., and Jacobs, D. J. (2014). Principal component analysis: a method for determining the essential dynamics of proteins. Methods Mol. Biol. 1084, 193–226. doi:10.1007/978-1-62703-658-0_11

PubMed Abstract | CrossRef Full Text | Google Scholar

Davis, S., and Meltzer, P. S. (2007). GEOquery: a bridge between the gene expression omnibus (GEO) and BioConductor. Bioinforma. Oxf. Engl. 23, 1846–1847. doi:10.1093/bioinformatics/btm254

PubMed Abstract | CrossRef Full Text | Google Scholar

Dietrich, N., Schneider, D. L., and Kornfeld, K. (2017). A pathway for low zinc homeostasis that is conserved in animals and acts in parallel to the pathway for high zinc homeostasis. Nucleic Acids Res. 45, 11658–11672. doi:10.1093/nar/gkx762

PubMed Abstract | CrossRef Full Text | Google Scholar

Du, J., Xu, M., Kong, F., Zhu, P., Mao, Y., Liu, Y., et al. (2022). CB2R attenuates intervertebral disc degeneration by delaying nucleus pulposus cell senescence through AMPK/GSK3β pathway. Aging Dis. 13, 552–567. doi:10.14336/AD.2021.1025

PubMed Abstract | CrossRef Full Text | Google Scholar

Fearing, B. V., Jing, L., Barcellona, M. N., Witte, S. E., Buchowski, J. M., Zebala, L. P., et al. (2019). Mechanosensitive transcriptional coactivators MRTF-A and YAP/TAZ regulate nucleus pulposus cell phenotype through cell shape. FASEB J. 33, 14022–14035. doi:10.1096/fj.201802725RRR

PubMed Abstract | CrossRef Full Text | Google Scholar

Francisco, V., Pino, J., González-Gay, M. Á., Lago, F., Karppinen, J., Tervonen, O., et al. (2022). A new immunometabolic perspective of intervertebral disc degeneration. Nat. Rev. Rheumatol. 18, 47–60. doi:10.1038/s41584-021-00713-z

PubMed Abstract | CrossRef Full Text | Google Scholar

Gagrica, S., Brookes, S., Anderton, E., Rowe, J., and Peters, G. (2012). Contrasting behavior of the p18INK4c and p16INK4a tumor suppressors in both replicative and oncogene-induced senescence. Cancer Res. 72, 165–175. doi:10.1158/0008-5472.CAN-11-2552

PubMed Abstract | CrossRef Full Text | Google Scholar

Ge, J., Yan, Q., Wang, Y., Cheng, X., Song, D., Wu, C., et al. (2020). IL-10 delays the degeneration of intervertebral discs by suppressing the p38 MAPK signaling pathway. Free Radic. Biol. Med. 147, 262–270. doi:10.1016/j.freeradbiomed.2019.12.040

PubMed Abstract | CrossRef Full Text | Google Scholar

Hänzelmann, S., Castelo, R., and Guinney, J. (2013). GSVA: gene set variation analysis for microarray and RNA-seq data. BMC Bioinforma. 14, 7. doi:10.1186/1471-2105-14-7

CrossRef Full Text | Google Scholar

Hao, Y., Hao, S., Andersen-Nissen, E., Mauck, W. M., Zheng, S., Butler, A., et al. (2021). Integrated analysis of multimodal single-cell data. Cell 184, 3573–3587.e29. doi:10.1016/j.cell.2021.04.048

PubMed Abstract | CrossRef Full Text | Google Scholar

He, J., Xue, R., Li, S., Lv, J., Zhang, Y., Fan, L., et al. (2015). Identification of the potential molecular targets for human intervertebral disc degeneration based on bioinformatic methods. Int. J. Mol. Med. 36, 1593–1600. doi:10.3892/ijmm.2015.2389

PubMed Abstract | CrossRef Full Text | Google Scholar

He, S., and Sharpless, N. E. (2017). Senescence in health and disease. Cell 169, 1000–1011. doi:10.1016/j.cell.2017.05.015

PubMed Abstract | CrossRef Full Text | Google Scholar

Huang, B., He, Y., Li, S., Wei, X., Liu, J., Shan, Z., et al. (2020). The RCAN1.4-calcineurin/NFAT signaling pathway is essential for hypoxic adaption of intervertebral discs. Exp. Mol. Med. 52, 865–875. doi:10.1038/s12276-020-0441-x

PubMed Abstract | CrossRef Full Text | Google Scholar

Kale, A., Sharma, A., Stolzing, A., Desprez, P.-Y., and Campisi, J. (2020). Role of immune cells in. the removal of deleterious senescent cells. Immun. Ageing 17, 16. doi:10.1186/s12979-020-00187-9

PubMed Abstract | CrossRef Full Text | Google Scholar

Kazezian, Z., Gawri, R., Haglund, L., Ouellet, J., Mwale, F., Tarrant, F., et al. (2015). Gene expression profiling identifies interferon signalling molecules and IGFBP3 in human degenerative annulus fibrosus. Sci. Rep. 5, 15662. doi:10.1038/srep15662

PubMed Abstract | CrossRef Full Text | Google Scholar

Kolde, R., and Kolde, M. R. (2015). Package ‘pheatmap’.

Google Scholar

Kupka, J., Kohler, A., El Bagdadi, K., Bostelmann, R., Brenneis, M., Fleege, C., et al. (2020). Adrenoceptor expression during intervertebral disc degeneration. Int. J. Mol. Sci. 21, 2085. doi:10.3390/ijms21062085

PubMed Abstract | CrossRef Full Text | Google Scholar

Lagnado, A., Leslie, J., Ruchaud-Sparagano, M.-H., Victorelli, S., Hirsova, P., Ogrodnik, M., et al. (2021). Neutrophils induce paracrine telomere dysfunction and senescence in ROS-dependent manner. EMBO J. 40, e106048. doi:10.15252/embj.2020106048

PubMed Abstract | CrossRef Full Text | Google Scholar

Lan, T., Hu, Z., Guo, W., Yan, B., and Zhang, Y. (2022). Development of a novel inflammatory-associated gene signature and immune infiltration patterns in intervertebral disc degeneration. Oxid. Med. Cell Longev. 2022, 2481071. doi:10.1155/2022/2481071

PubMed Abstract | CrossRef Full Text | Google Scholar

Lanna, A., Vaz, B., D’Ambra, C., Valvo, S., Vuotto, C., Chiurchiù, V., et al. (2022). An intercellular transfer of telomeres rescues T cells from senescence and promotes long-term immunological memory. Nat. Cell Biol. 24, 1461–1474. doi:10.1038/s41556-022-00991-z

PubMed Abstract | CrossRef Full Text | Google Scholar

Leek, J. T., Johnson, W. E., Parker, H. S., Jaffe, A. E., and Storey, J. D. (2012). The sva package for removing batch effects and other unwanted variation in high-throughput experiments. Bioinformatics 28, 882–883. doi:10.1093/bioinformatics/bts034

PubMed Abstract | CrossRef Full Text | Google Scholar

Le Maitre, C. L., Freemont, A. J., and Hoyland, J. A. (2007). Accelerated cellular senescence in degenerate intervertebral discs: a possible role in the pathogenesis of intervertebral disc degeneration. Arthritis Res. Ther. 9, R45. doi:10.1186/ar2198

PubMed Abstract | CrossRef Full Text | Google Scholar

Li, N., Gao, Q., Zhou, W., Lv, X., Yang, X., and Liu, X. (2020). MicroRNA-129-5p affects immune privilege and apoptosis of nucleus pulposus cells via regulating FADD in intervertebral disc degeneration. Cell Cycle 19, 933–948. doi:10.1080/15384101.2020.1732515

PubMed Abstract | CrossRef Full Text | Google Scholar

Li, Q., Nirala, N. K., Chen, H. J., Nie, Y., Wang, W., Zhang, B., et al. (2019). The Misshapen subfamily of Ste20 kinases regulate proliferation in the aging mammalian intestinal epithelium. J. Cell. physiology 234, 21925–21936. doi:10.1002/jcp.28756

PubMed Abstract | CrossRef Full Text | Google Scholar

Li, W., Zhao, Y., Wang, Y., He, Z., Zhang, L., Yuan, B., et al. (2023). Deciphering the sequential changes of monocytes/macrophages in the progression of IDD with longitudinal approach using single-cell transcriptome. Front. Immunol. 14, 1090637. doi:10.3389/fimmu.2023.1090637

PubMed Abstract | CrossRef Full Text | Google Scholar

Liberzon, A., Birger, C., Thorvaldsdóttir, H., Ghandi, M., Mesirov, J. P., and Tamayo, P. (2015). The Molecular Signatures Database (MSigDB) hallmark gene set collection. Cell Syst. 1, 417–425. doi:10.1016/j.cels.2015.12.004

PubMed Abstract | CrossRef Full Text | Google Scholar

Lin, H., Wang, Y., Jing, K., Wu, T., Niu, Y., and Wei, J. (2020). Nuclear factor erythroid-2 related factor 2 inhibits human disc nucleus pulpous cells apoptosis induced by excessive hydrogen peroxide. Rev. Assoc. Med. Bras. 66, 986–991. doi:10.1590/1806-9282.66.7.986

PubMed Abstract | CrossRef Full Text | Google Scholar

Muñoz-Espín, D., and Serrano, M. (2014). Cellular senescence: from physiology to pathology. Nat. Rev. Mol. Cell Biol. 15, 482–496. doi:10.1038/nrm3823

PubMed Abstract | CrossRef Full Text | Google Scholar

Park, S. S., Lee, Y.-K., Park, S. H., Lim, S. B., Choi, Y. W., Shin, J. S., et al. (2023). p15INK4B is an alternative marker of senescent tumor cells in colorectal cancer. Heliyon 9, e13170. doi:10.1016/j.heliyon.2023.e13170

PubMed Abstract | CrossRef Full Text | Google Scholar

Patil, P., Dong, Q., Wang, D., Chang, J., Wiley, C., Demaria, M., et al. (2019). Systemic clearance of p16INK4a -positive senescent cells mitigates age-associated intervertebral disc degeneration. Aging Cell 18, e12927. doi:10.1111/acel.12927

PubMed Abstract | CrossRef Full Text | Google Scholar

Pedersen, L. M., Schistad, E., Jacobsen, L. M., Røe, C., and Gjerstad, J. (2015). Serum levels of the pro-inflammatory interleukins 6 (IL-6) and -8 (IL-8) in patients with lumbar radicular pain due to disc herniation: a 12-month prospective study. Brain Behav. Immun. 46, 132–136. doi:10.1016/j.bbi.2015.01.008

PubMed Abstract | CrossRef Full Text | Google Scholar

Pfirrmann, C. W. A., Metzdorf, A., Zanetti, M., Hodler, J., and Boos, N. (2001). Magnetic. Resonance. Classification of lumbar intervertebral disc degeneration. Spine 26, 1873–1878. doi:10.1097/00007632-200109010-00011

PubMed Abstract | CrossRef Full Text | Google Scholar

Phillips, K. L. E., Chiverton, N., Michael, A. L. R., Cole, A. A., Breakwell, L. M., Haddock, G., et al. (2013). The cytokine and chemokine expression profile of nucleus pulposus cells: implications for degeneration and regeneration of the intervertebral disc. Arthritis Res. Ther. 15, R213. doi:10.1186/ar4408

PubMed Abstract | CrossRef Full Text | Google Scholar

Prata, L. G. P. L., Ovsyannikova, I. G., Tchkonia, T., and Kirkland, J. L. (2018). Senescent cell clearance by the immune system: emerging therapeutic opportunities. Seminars Immunol. 40, 101275. doi:10.1016/j.smim.2019.04.003

PubMed Abstract | CrossRef Full Text | Google Scholar

Qaseem, A., Wilt, T. J., McLean, R. M., Forciea, M. A., and Denberg, T. D. (2017). Noninvasive treatments for acute, subacute, and chronic low back pain: a clinical practice guideline from the American college of physicians. Ann. Intern Med. 166, 514–530. doi:10.7326/M16-2367

PubMed Abstract | CrossRef Full Text | Google Scholar

Reizis, B. (2019). Plasmacytoid dendritic cells: development, regulation, and function. Immunity 50, 37–50. doi:10.1016/j.immuni.2018.12.027

PubMed Abstract | CrossRef Full Text | Google Scholar

Ritchie, M. E., Phipson, B., Wu, D., Hu, Y., Law, C. W., Shi, W., et al. (2015). Limma powers differential expression analyses for RNA-sequencing and microarray studies. Nucleic Acids Res. 43, e47. doi:10.1093/nar/gkv007

PubMed Abstract | CrossRef Full Text | Google Scholar

Ru, Y., Kechris, K. J., Tabakoff, B., Hoffman, P., Radcliffe, R. A., Bowler, R., et al. (2014). The multiMiR R package and database: integration of microRNA-target interactions along with their disease and drug associations. Nucleic Acids Res. 42, e133. doi:10.1093/nar/gku631

PubMed Abstract | CrossRef Full Text | Google Scholar

Saberi, M., Zhang, X., and Mobasheri, A. (2021). Targeting mitochondrial dysfunction with small molecules in intervertebral disc aging and degeneration. Geroscience 43, 517–537. doi:10.1007/s11357-021-00341-1

PubMed Abstract | CrossRef Full Text | Google Scholar

Shannon, P., Markiel, A., Ozier, O., Baliga, N. S., Wang, J. T., Ramage, D., et al. (2003). Cytoscape: a software environment for integrated models of biomolecular interaction networks. Genome Res. 13, 2498–2504. doi:10.1101/gr.1239303

PubMed Abstract | CrossRef Full Text | Google Scholar

Simon, N., Friedman, J., Hastie, T., and Tibshirani, R. (2011). Regularization paths for cox’s proportional hazards model via coordinate descent. J. Stat. Softw. 39, 1–13. doi:10.18637/jss.v039.i05

CrossRef Full Text | Google Scholar

Song, C., Cai, W., Liu, F., Cheng, K., Guo, D., and Liu, Z. (2022). An in-depth analysis of the immunomodulatory mechanisms of intervertebral disc degeneration. JOR Spine 5, e1233. doi:10.1002/jsp2.1233

PubMed Abstract | CrossRef Full Text | Google Scholar

Stich, S., Stolk, M., Girod, P. P., Thomé, C., Sittinger, M., Ringe, J., et al. (2015). Regenerative and immunogenic characteristics of cultured nucleus pulposus cells from human cervical intervertebral discs. PLoS One 10, e0126954. doi:10.1371/journal.pone.0126954

PubMed Abstract | CrossRef Full Text | Google Scholar

Subramanian, A., Tamayo, P., Mootha, V. K., Mukherjee, S., Ebert, B. L., Gillette, M. A., et al. (2005). Gene set enrichment analysis: a knowledge-based approach for interpreting genome-wide expression profiles. Proc. Natl. Acad. Sci. U. S. A. 102, 15545–15550. doi:10.1073/pnas.0506580102

PubMed Abstract | CrossRef Full Text | Google Scholar

Tu, J., Li, W., Yang, S., Yang, P., Yan, Q., Wang, S., et al. (2022). Single-cell transcriptome profiling reveals multicellular ecosystem of nucleus pulposus during degeneration progression. Adv. Sci. (Weinh) 9, e2103631. doi:10.1002/advs.202103631

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, L., He, T., Liu, J., Tai, J., Wang, B., Zhang, L., et al. (2021). Revealing the immune infiltration landscape and identifying diagnostic biomarkers for lumbar disc herniation. Front. Immunol. 12, 666355. doi:10.3389/fimmu.2021.666355

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, Z., Zhang, S., Zhao, Y., Qu, Z., Zhuang, X., Song, Q., et al. (2021). MicroRNA-140-3p alleviates intervertebral disc degeneration via KLF5/N-cadherin/MDM2/Slug axis. RNA Biol. 18, 2247–2260. doi:10.1080/15476286.2021.1898176

PubMed Abstract | CrossRef Full Text | Google Scholar

Wen, T., Xue, P., Ying, J., Cheng, S., Liu, Y., and Ruan, D. (2021). The role of unfolded protein response in human intervertebral disc degeneration: perk and IRE1-α as two potential therapeutic targets. Oxid. Med. Cell Longev. 2021, 6492879. doi:10.1155/2021/6492879

PubMed Abstract | CrossRef Full Text | Google Scholar

Wickham, H. (2009). ggplot2: elegant graphics for data analysis. New York, NY: Springer New York. doi:10.1007/978-0-387-98141-3

CrossRef Full Text | Google Scholar

Wu, J., Chen, Y., Liao, Z., Liu, H., Zhang, S., Zhong, D., et al. (2022). Self-amplifying loop of NF-κB and periostin initiated by PIEZO1 accelerates mechano-induced senescence of nucleus pulposus cells and intervertebral disc degeneration. Mol. Ther. 30, 3241–3256. doi:10.1016/j.ymthe.2022.05.021

PubMed Abstract | CrossRef Full Text | Google Scholar

Xin, J., Wang, Y., Zheng, Z., Wang, S., Na, S., and Zhang, S. (2022). Treatment of intervertebral disc degeneration. Orthop. Surg. 14, 1271–1280. doi:10.1111/os.13254

PubMed Abstract | CrossRef Full Text | Google Scholar

Xu, S., Fu, H., Weng, S., Gu, X., and Li, J. (2023). Derivation and comprehensive analysis of ageing-related genes in intervertebral disc degeneration for prediction and immunology. Mech. Ageing Dev. 211, 111794. doi:10.1016/j.mad.2023.111794

PubMed Abstract | CrossRef Full Text | Google Scholar

Yang, C., Cao, P., Gao, Y., Wu, M., Lin, Y., Tian, Y., et al. (2016). Differential expression of p38 MAPK α, β, γ, δ isoforms in nucleus pulposus modulates macrophage polarization in intervertebral disc degeneration. Sci. Rep. 6, 22182. doi:10.1038/srep22182

PubMed Abstract | CrossRef Full Text | Google Scholar

Yao, C., Ren, J., Huang, R., Tang, C., Cheng, Y., Lv, Z., et al. (2022). Transcriptome profiling of microRNAs reveals potential mechanisms of manual therapy alleviating neuropathic pain through microRNA-547-3p-mediated Map4k4/NF-κb signaling pathway. J. Neuroinflammation 19, 211. doi:10.1186/s12974-022-02568-x

PubMed Abstract | CrossRef Full Text | Google Scholar

Yu, D., Hu, J., Sheng, Z., Fu, G., Wang, Y., Chen, Y., et al. (2020). Dual roles of misshapen/NIK-related kinase (MINK1) in osteoarthritis subtypes through the activation of TGFβ signaling. Osteoarthr. Cartil. 28, 112–121. doi:10.1016/j.joca.2019.09.009

PubMed Abstract | CrossRef Full Text | Google Scholar

Yu, G., Wang, L.-G., Han, Y., and He, Q.-Y. (2012). clusterProfiler: an R package for comparing biological themes among gene clusters. OMICS 16, 284–287. doi:10.1089/omi.2011.0118

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhang, G.-Z., Liu, M.-Q., Chen, H.-W., Wu, Z.-L., Gao, Y.-C., Ma, Z.-J., et al. (2021). NF-κB signalling pathways in nucleus pulposus cell function and intervertebral disc degeneration. Cell Prolif. 54, e13057. doi:10.1111/cpr.13057

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhou, G., Soufan, O., Ewald, J., Hancock, R. E. W., Basu, N., and Xia, J. (2019). NetworkAnalyst 3.0: a visual analytics platform for comprehensive gene expression profiling and meta-analysis. Nucleic Acids Res. 47, W234–W241. doi:10.1093/nar/gkz240

PubMed Abstract | CrossRef Full Text | Google Scholar