Nucleus pulposus cell senescence is abnormally activated in intervertebral disc degeneration

To investigate the activation state of cell senescence in IVDD, we firstly extracted NP cells from the public dataset PRJCA014236 and performed geneset scoring with the universal senescence geneset Senmayo. Interestingly, among the five cell types, NP cells did not exhibit the highest senescence score at baseline (Fig. 1A–B), which indicated that a high baseline NP cell senescence score is not a specific characteristic of intervertebral disc degeneration. We further compared the senescence score of NP cells in different grades of IVDD and found the Senmayo score was significantly elevated as the grade of NP cell degeneration increases (Fig. 1C–E). At the same time, we integrated the senescence score in different samples and separated the samples into a high score group and a low score group, with 9 cases in the former and 5 cases in the latter (Fig. 1F). Interestingly, higher scores were mainly distributed in samples with high levels of degeneration. The UMAP plot revealed the characteristics of cells in the intervertebral disc microenvironment in patients with high scores and low scores (Fig. 1G). The stacking plots exhibited that patients in the high score group have higher infiltration levels of most cell types, including Annulus fibrosus, Endothelial cells, and immune cells (Fig. 1G), except for nucleus pulposus cells. Considering NP cells are the initiating cells and functional cells during the degeneration of the intervertebral disc, we extracted NP cells for subsequent analyses.

Fig. 1

The senescence landscape of NP cells in PRJCA014236 dataset. A UMAP plot of Senmayo score in PRJCA014236 dataset. B Violin plot of Senmayo score in different cell types. C Violin plot of Senmayo score in control NP cells and moderately degenerated NP cells. D Violin plot of Senmayo score in control NP cells and severely degenerated NP cells. E Violin plot of Senmayo score in moderately and severely degenerated NP cells. F Senmayo score of different samples. G The UMAP plot and stacking plot revealing the characteristic and proportion of high (left) and low (right) Senmayo score group

Identification of hub genes and characteristics of high senescence score

To explore the underlying hub genes in high senescence NP cells, we conducted hdWGCNA with the Senmayo score. We firstly selected the optimal soft threshold to make excellent consistency between the network’s topology and the biological information. Subsequently, we constructed the hierarchical co-expression network visualized by the dendrogram through measuring gene expression similarity, calculating the topological overlap matrix, and performing hierarchical clustering (Fig. 2A). Next, we obtained representative genes from each module, and the correlation of each gene within modules was exhibited by the connectivity plot (Fig. 2B). Meanwhile, we calculated the correlation between modules and senescence scores and found that Module-5, 14, 16, and 18 were significantly correlated with high senescence scores (Fig. 2C). We then extracted the top 25 central genes in the above modules and applied function enrichment analyses in GO and KEGG databases. In GO biological process pathways, the classical pathways involving inflammation and matrix metabolism, such as cellular responses to tumor necrosis factor (TNF) and interleukin-1, hyaluronan metabolic process was highly enriched. Meanwhile, the collagen-containing extracellular matrix (ECM), endoplasmic reticulum lumen, and similar cellular components participating in the degeneration of NP cells were also significantly enriched [26, 27]. In the category of molecular function, the cytokine activity [28] and extracellular matrix structural constituent involved in the senescence of NP cells were also activated (Fig. 2D). Moreover, the hub genes also showed prominent associations with pathways involved in senescence and ECM metabolism of NP cells, such as TNF signaling pathway [5], transforming growth factor-β (TGF-β) signaling pathway [29] and ECM-receptor interaction pathway [30] (Fig. 2E), which verified the correlation between the hub genes and senescence of NP cells.

Fig. 2

HdWGCNA analysis of NP cells with high Senmayo score in PRJCA014236 dataset. A Co-expression network diagram of NP cell with high Senmayo score. B Modules and associated top hub genes displayed according to the hdWGCNA pipeline. C Correlation diagram between modules and Senmayo score. D Bar diagram of GO enrichment of hub genes, including biological process (BP), cellular component (CC) and molecular function (MF). E Alluvial diagram and bubble diagram of KEGG enrichment analysis of hub genes

Validation of the abnormal activation of NP cell's senescence in IVDD

To validate the robustness of the above outcomes, we collected another single cell dataset, GSE244889. After quality control, dimensionality reduction, and clustering, we also annotated cell types with specific cell markers (Fig. 3A–B). Consistently, there was no high baseline NP cell senescence score based on the Senmayo geneset in this dataset (Fig. 3C), and the activated senescence state was also observed in high-grade degeneration samples (Fig. 3D–F). Subsequently, we conducted hdWGCNA analysis in this dataset and found Module-2, 3, 4, 7, 8, 9, 10, and 12 were positively correlated with senescence scores (Fig. 3G–I). Function enrichment analysis revealed the extracted hub genes were similarly enriched in the pathway regarding specific characteristics of NP cells, such as cartilage development, extracellular matrix organization, and structural constituent (Fig. 4A). Meanwhile, the above hub genes were also associated with pathways involved in the activation of senescence in the nucleus pulposus, for example, secretory granule lumen [31] and ATP-dependent protein folding chaperone [32]. KEGG enrichment analysis gave similar results that TGF-β signaling pathway, Longevity regulating pathway, and ECM-receptor interaction were obviously enriched (Fig. 4B), which validated the effectiveness of the above hub genes. We then intersected the hub genes in two datasets and obtained 15 genes including NBL1, CRTAC1, ANK3, SOD2, BMP2, COL15A1, TNC, TMSB4X, PAPSS2, IL1R1, GLIS3, CD44, INHBA, GPRC5A, and IER3 (Fig. 4C and Supplemental Table 3).

Fig. 3

The senescence landscape of NP cells in GSE244889 dataset. A UMAP plot of cell types in GSE244889 dataset. B UMAP plot of Senmayo score in GSE244889 dataset. C Violin plot of Senmayo score in different cell types. D Violin plot of Senmayo score in control NP cells and moderately degenerated NP cells. E Violin plot of Senmayo score in control NP cells and severely degenerated NP cells. F Violin plot of Senmayo score in moderately and severely degenerated NP cells. G Co-expression network diagram of NP cell with high Senmayo score. H Modules and associated top hub genes displayed according to the hdWGCNA pipeline. I Correlation diagram of modules and Senmayo score

Fig. 4

Biological characteristics of hub genes in the GSE244889 dataset. A Bar diagram of GO enrichment analysis of hub genes. B Alluvial diagram and bubble diagram of KEGG enrichment analysis of hub genes. C Venn diagram of common hub genes (genes indicated by the arrow only for exhibition)

Construction of the novel signature called NP_Senescence

Considering the classic Senmayo geneset containing too many genes, which limited it applicability, we constructed a novel signature for nucleus pulposus called NP_Senescence containing above 15 genes. We validated this signature in above two single cell datasets, and found the significant increase of the score in the high Pfirrmann grade samples. Compared with the Senmayo score, the difference in NP_Senescence score between different grades was more obviously, not only in PRJCA014236 (Fig. 5A–C), but also in GSE244889 (Fig. 5D–F). For verifying the effectiveness of the score, we applied SenCID analysis, a machine learning program which could accurately and automatically identify six major senescence identities (SID1-6) and calculate corresponding identity scores (SenCID scores) spanning 30 cell types, including chondrocytes. In PRJCA014236 dataset, the SenCID scores across five identities exhibited an overall increase trend with degeneration levels, except SID2 (Fig. 5G). The vast majority of NP cells were identified as the SID3 category, which was reported to show enrichment for mitochondria and redox reactions. The remaining one thirds of NP cells were recommended as the SID5 category with the specific characteristic of high activity about BCL2 family genes, which participated in the senescence-related NP degeneration [33]. Moreover, consistent with above results, NP cells in GSE244889 dataset were also divided into SID3 pattern, which also showed an increasing trend as the grade of degradation increases (Fig. 5H). Interestingly, it is reported that among the 30 cell types, chondrocytes in osteoarthritis were also detected as SID3 cells, which, to a certain extent, explained the homologous between chondrocytes and NP cells. We further calculated the correlation between NP_Senescence and SenCID score of corresponding categories, and the results showed a significant correlation in both two datasets (Fig. 5I–K), which indicated the effectiveness of our novel signature NP_Senencence.

Fig. 5

Quantification system based on the novel signature NP_Senescence. A–C Violin plot of NP_Senescence in PRJCA014236 dataset. D–F Violin plot of NP_Senescence in GSE244889 dataset. G Boxplot of SenCID score in PRJCA014236 dataset. H Boxplot of SenCID score in GSE244889 dataset. I Correlation diagram of NP_Senescence and SID3 score in PRJCA014236 dataset. J Correlation diagram of NP_Senescence and SID5 score in PRJCA014236 dataset. K Correlation diagram of NP_Senescence and SID3 score in GSE244889 dataset

Identification of metabolism-specific phenotypes in senescent NP cells

Cell senescence is closely related to metabolic regulation, which can shape the specific senescent phenotype. The imbalance in the redox system, proteostasis, and organelle homeostasis, as well as changes in related metabolites, are involved in the regulation of cell senescence. Therefore, we adopted scFEA, a single-cell level flux estimation algorithm based on a multi-layer neural network, to calculate the underlying metabolites in the process of NP cell senescence for exploring the specific metabolic pattern. Based on the NP_Senescence score, we calculated the metabolic flux of the samples and groups. From an overall perspective, in the PRJCA014236 dataset, metabolic pathways including hyaluronic acid synthesis, transporters, spermine metabolism, and purine synthesis were significantly activated in the high NP_Senescence group (Fig. 6A). In the process of hyaluronic acid synthesis, the synthesis of two monomers, UDP-N-acetylglucosamine and UDP-glucuronic acid, required for hyaluronic acid was obviously activated (Fig. 6B). However, this is inconsistent with our common knowledge that extracellular matrix anabolism is downregulated in senescent cells, which indicated that negative feedback may increase monomer synthesis in senescent cells, but hyaluronic acid synthesis and transport may be blocked to some extent during senescence. Moreover, the transport of multiple metabolites, like tyrosine and arginine, and spermine metabolism were also upregulated (Fig. 6C–D). Cells are driven into the senescent state by a series of factors, including telomere shortening and DNA damage [34]. Increased synthesis of β-alanine and deoxycytidine in pyrimidine metabolism was observed in senescent cells (Fig. 6E), which suggested that abnormal DNA damage repair may be in an active state. At the same time, in purine synthesis, there appeared to be an increased demand for various purines in high NP_Senescence score cells, compared with the low score senescent cells. However, the synthesis of deoxyadenosine, whose derivatives have been reported to have anti-senescence effects [35], was significantly inhibited (Fig. 6F). Afterwards, correlation analysis between metabolite concentrations and NP_Senescence revealed that adenosine monophosphate (AMP) and glucose-1-phosphate were negatively correlated with NP_Senescence (Fig. 6G–H), while β-alanine and UDP-glucuronic acid were positively correlated (Fig. 6I–J). In the GSE244889 dataset, the main activated metabolic pathways included hyaluronic acid synthesis, transporters, glycolysis, TCA cycle, and purine synthesis, while pyrimidine synthesis was partially inhibited (Fig. 7A). Similar to the results of the PRJCA014236 dataset, the synthesis of UDP-N-acetylglucosamine and UDP-glucuronic acid was enhanced, as well as glucose-6-phosphate and glucose-1-phosphate, in the process of hyaluronic acid synthesis (Fig. 7B). The transport of multiple metabolites, glycolysis, TCA cycle, and purine synthesis were observed to be abnormally active, while the process of pyrimidine synthesis exhibits dual effects of both activation and inhibition (Fig. 7C–F). Meanwhile, the concentration of AMP and glucose-1-phosphate also showed a negative correlation with NP_Senescence, while the level of β-alanine and UDP-glucuronic acid elevated as the novel signature increased (Fig. 7G–J). On the whole, the senescent NP cells recognized by our novel score may be blocked during hyaluronic acid synthesis and transport and exhibit the characteristic of extraordinarily active DNA mismatch repair.

Fig. 6

Metabolic characteristic of NP cells based on scFEA algorithm in PRJCA014236 dataset. A Heatmap of specific metabolic landscape of NP cells grouped by NP_Senescence. B–F Volcano diagram of metabolic pathways hyaluronic acid synthesis, transporters, spermine metabolism, pyrimidine synthesis and purine synthesis. G–J Correlation diagrams of NP_Senescence with AMP, glucose-6-phosphate, β-aianine and UDP-glucuronic acid

Fig. 7

Metabolic characteristic of NP cells based on scFEA algorithm in GSE244889 dataset. A Heatmap of specific metabolic landscape of NP cells grouped by NP_Senescence. B–F Volcano diagram of metabolic pathways hyaluronic acid synthesis, transporters, glycolysis TCA cycle, pyrimidine synthesis and purine synthesis. G–J Correlation diagrams of NP_Senescence with AMP, glucose-6-phosphate, β-aianine and UDP-glucuronic acid

External validation on diagnostic efficacy of NP_Senescence

For validating the robustness of the novel signature, we collected the other datasets about IVDD, such as GSE70362, GSE230809, and GSE34095 datasets, to validate the diagnostic efficacy. We found four markers in the signature, namely IL1R1, SOD2, CD44, and CRTAC1, also had a higher predictive accuracy of diagnosis on the external validation set GSE70362, respectively the AUC values of 0.68, 0.773, 0.641, and 0.633 (Fig. 8A–D). On the other hand, in the external single-cell dataset, GSE230809, we also extracted NP cells after conventional dimensionality reduction, clustering, and annotation (Fig. 8E). The specific signature was still observed excellent predictive efficacy, with elevated senescence scores in the high-grade degeneration group (Fig. 8F). Meanwhile, the expression of the markers found above, like IL1R1 and SOD2, was also significantly increased in the severely degenerated group (Fig. 8G–H), which verified their potential to be diagnostic markers. We further validated the effectiveness of our senescence evaluating system in another dataset GSE34095, which only included three pairs of degenerated intervertebral disc and non-degenerated intervertebral disc samples. Due to the limited patient number, the results indicated the trend of increased expression in markers IL1R1 and CRTAC1 regardless of insignificant statistical difference (Fig. 8I–J). Together, these findings indicate that NP_Senescence is a potential and stable predictor of the diagnosis in IVDD.

Fig. 8

Validation of the NP_Senescence in external validation set. A–D ROC curves of the marker genes IL1R1, SOD2, CD44 and CRTAC1 in GSE70362 datasets. E UMAP plot of cell types in GSE230809 dataset. F Violin plot of NP_Senescence in GSE230809 dataset. G–H Violin plots of IL1R1 and SOD2 expressions in GSE230809 dataset. I–J Box plots of IL1R1 and SOD2 expressions in GSE34095 dataset