NP samples from four patients in lumbar discectomy surgery were collected and minced with a razor blade into 1 mm³ fragments on a dish under aseptic conditions. NP fragments were suspended in 10 ml of 1×TrypLE Express and incubated in 37°C water bath with shaking for 30 min. Centrifuge the tissue at 500g for 5 minutes at room temperature and wash twice by adding 25 ml PBS. Discard the supernatant and add 10 mL of digestion buffer consisting of 2 mg/mL Collagenase type II (Sigma, c6885), 0.1mg/ml DNase I (Roche, 4716728001) in RPMI medium (Coring, 10-040-CV), and incubated in 37°C water bath with shaking for 4-6 hours. Terminate the enzymatic digestion when the tissue pieces have entirely disintegrated.

The suspension was passed through a 100 μm filter (Falcon, 3523260) and centrifuged (400g, 10 min, 4°C). Pelleted cells were resuspended in red blood cell lysis buffer (Solarbio, R1010) for 5 min at room temperature, followed by an additional filter step through a 40 μm filter (Falcon,3523240). Cells were collected by centrifugation (400g, 10 min, 4°C) and resuspended in PBS (BI, 02-024-1ACS) containing 0.04% BSA (Sigma, B2064).

Cells were manually counted by Trypan blue (Thermo, T10282) and AO-PI (LUNA, D23001) after each centrifugation (400g, 10 min, 4°C) and resuspended. Single cells were processed using Chromium Controller (10× Genomics) according to the manufacturer’s protocol.

By using Chromium Next GEM Single Cell 3′ Kit v3.1 (10x Genomics, 1000268) and Chromium Next GEM Chip G Single Cell Kit (10xGenomics, 1000120), we performed single cell 3′ gene expression profiling. The cell suspension was loaded onto the Chromium single cell controller (10x Genomics) to generate single-cell gel beads in the emulsion according to the manufacturer’s protocol. Captured cells were lysed and the released RNA were barcoded through reverse transcription in individual GEMs. Cell-barcoded 3′ gene expression libraries were sequenced on an Illumina NovaSeq6000 system by Shanghai Biochip Co., Ltd., Shanghai, China.

Raw sequencing data was processed by Cell Ranger (10xGenomics, version cellranger-6.0.0) pipeline and aligned to the human reference genome (GRCh38) to generate Gene-Barcode matrices of gene expression. The Gene-Barcode matrices containing the barcoded cells and gene expression counts were imported into the Seurat R toolkit (R version 4.0.3 (2020–10–10)). Cells with gene number (<300 & > 6000) or high mitochondrial transcript ratio (<25%), and genes expressed in less than 3 cells were all excluded. After removing unwanted cells from the dataset, all samples were combined with function “merge”. Next, we employed a global-scaling normalization method “LogNormalize” to normalize the feature expression measurements (UMI counts) for each cell by the total expression and the canonical correlation analysis (CCA) method to remove the batch effect.

Highly variable genes (top 4000) were extracted to perform the principal component analysis (PCA) and top 30 of significant principal components were used for cluster analysis. Clusters were visualized using the Uniform Manifold Approximation and Projection (UMAP) and t-Distributed Stochastic Neighbor Embedding (t-SNE). Marker genes for each cluster, cell type and subgroup were identified by contrasting gene expression of cells from certain cluster, cell type or subgroup to that of others using the Seurat FindMarkers function. The types of cell clusters were annotated by specific marker genes of published cell types.

GO enrichment, KEGG enrichment of markers from cluster or cell type were performed using clusterProfiler software with Benjamini-Hochberg multiple testing adjustment, using marker genes from Seurat software with wilcox test and log-scale foldchange (logfc.threshold >= 0.25). The results were visualized using R package.

GSEA was performed by using clusterProfiler software, which used predefined gene sets from KEGG Pathway Database. The genes in selected gene sets from distinct cell types of four samples were used. Gene expression data was calculated by mean umi count of genes in one cell type and the rest cell types, respectively.

Cell types of blood cell and stromal cell subclusters were annotated by SingleR (1.4.1). Unbiased cell type recognition from single-cell RNA sequencing data was performed by leveraging reference transcriptomic datasets of pure cell types to infer the origin of each single cell independently. For human, (HumanPrimaryCellAtlasData) was used.

To obtain the cell trajectories of stromal cells, the scVelo (0.2.0) package with default parameter was applied for studying the cellular differentiation status based on the normalized cell matrix from selected cells to solve the transcriptional dynamics of splicing kinetics using a likelihood-based dynamical model.

To enable a systematic analysis of cell-cell communication molecules, cell communication analysis based on the CellPhoneDB, a public repository of ligands, receptors, and their interactions, was applied. Significant mean and cell communication significance (p-value<0.05) was calculated based on the interaction and the normalized cell matrix achieved by Seurat Normalization. To identify potential interactions between and within immune cells and stromal cells, we used CellPhoneDB 2.0 with parameters threshold =0.25 and iterations=1000, which contains a curated repository of ligand-receptor interactions and a statistical framework for inferring lineage-specific interactions.

Samples were sectioned at a thickness of 5 µm. All the sections were incubated for 1 hour with primary antibody directed against IL1B (Cell Signaling Technology, #2022, 1:200), COL2A1 (Santa Cruz, sc-52658, 1:200), COL4A1 (Thermo Fisher Scientific, PA5-85634, 1:500), LGALS1 (Cell Signaling Technology, #5418, 1:200), MMP3 (Abcam, ab223666, 1:500) after blocking endogenous peroxidase using 3% hydrogen peroxide for 5 minutes at room temperature. After rinsing, the sections were incubated for 1 hour with biotinylated horseradish peroxidase-conjugated goat antirabbit IgG. Diaminobenzidine was used to develop peroxidase staining.

NP samples were placed in a 37° incubator for digestion for a half hour using 0.25% trypsin. After centrifugation at 1500 RPM/min, the supernatant was removed, 0.1% type II collagenase (Thermo Fisher Scientific, 17101015) was added, and the cells were collected in an incubator at 37°C overnight. Cultured cells were validated by immunofluorescence staining with aggrecan (Abcam, ab36861, 1:500), which was NP cell-specific markers ( Supplementary Figure S1 ).

After treated by the recombinant proteins such as TNF-α (MedChemExpress, HY-P7058), CD74 (MedChemExpress, HY-P72931) and CCL-3 (MedChemExpress, HY-P7256) for 24 hours, NP cells were further processed for immunofluorescence staining detection of the ossification markers including COL1A1, and RUNX2. Twenty thousand living cells were loaded onto each coverslip and incubated overnight. PFA (4%) was fixed at room temperature for 15 min. Then, blocking solution (containing 10% horse serum, 1% BSA and 0.1% Triton X-100 DPBS) was added and blocked for 1 hour. The blocking solution was discarded, and diluted primary antibodies including COL1A1 (ABclonal, A1352, 1:100), and RUNX2 (Thermo Fisher Scientific, 41-1400, 1:200) were added for incubation. After rinsing, the coverslips were placed in the dark and incubated in florescent-labeled goat-anti-rabbit IgG for 1 hour. Then, prolong gold antifade reagent with DAPI was added to the coverslips. All coverslips were randomly photographed, and at least 200 cells were counted for analysis (Leica DMI 6000B).

Research protocols involving animals were approved by the Ethics Committee of Shanghai Jiaotong University First People’s Hospital (2021AWS0022). Twenty-four rats’ coccyx 6-7 intervertebral discs were punctured by a 26-G syringe needle to establish the degeneration group (annulus fibrous puncture group: AFP group) (11). For twelve rats, the TNF-α specific inhibitor (MedChemExpress, Infliximab: HY-P9970) were injected during the puncture process. After 4 weeks, all the rats were anesthetized for the further analysis. X-rays and MRI were performed to determine the degeneration degree in different groups. Moreover, the rats were sacrificed at 12 weeks after puncture. X-rays and MRI were performed again to evaluate the degeneration of rat intervertebral disc. Additionally, HE staining and safranin O-fast green staining were further used to validate the degeneration of rat intervertebral disc. Moreover, we also detected the expression of TNF-α (Cell Signaling Technology, #3707, 1:500), Col-1 (Abcam, ab34710, 1:200), Runx-2 (Abcam, ab192256, 1:200) in vivo after harvesting disc tissues using immunohistochemical analysis.

SPSS V.18.0 software was utilized for statistical analysis. The measurement data were expressed as mean ± standard deviation. Differences between two groups were analyzed using two-tailed Student’s t-test. One-way analysis of variance was employed to compare all parameters among multiple groups, and LSD tests were performed for pairwise comparisons within groups. P<0.05 was considered statistically significant.

We dissociated NP tissues from four IDD patients with different degrees of degeneration and performed scRNA-seq ( Figure 1A , Table 1 ). After quality control, 18,855 nuclear pulposus cells (NPCs) with a median of 1053 genes per cell remained and were clustered into sixteen populations. Immune cell marker genes (CCL3L1, CCL3, and IL1B) were ubiquitously expressed in the blood clusters (Clusters 4, 11, 13, 14, and 15) (12). The specific genes (COL2A1, PLA2G2A, and SCRG1) of in Clusters 1, 3, and 6 were chondrocyte marker genes (13). COL4A1, IFI27, and COL4A2 were the most highly expressed in endothelial cells (ECs) (Cluster 10) (14). Pericytes (Clusters 9 and 12) were identified by marker gene expression (LGALS1, IBSP, and TAGLN) (15). Stromal cells (Clusters 0, 2, 5, 7, and 8) were mainly identified by the expression of genes including COL1A1, POSTN, and MMP3 ( Figure 1B ) (16). These marker genes are shown in feature plots, ridge plots and heatmaps ( Figures 1C–E ). Immunohistochemistry validated the expression of specific markers for each NP cell cluster at the protein level ( Figure 1F ). Interestingly, the UMAP shows the heterogeneity among these samples; therefore, the proportions of different cell clusters were calculated among the four patients. We found no pericytes or blood cells in Patient 4, which was probably associated with a reduction in vessel ingrowth in degenerative NP due to the young age and short course of disease in this patient ( Figure 1G ) (17). Overall, these results revealed the atlas of human NP during IDD pathogenesis at the single-cell level.

To investigate the potential functions of the differentially expressed genes (DEGs) identified in the five cell types, we used the Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) databases to explore the biological processes and pathways, respectively. The DEGs of blood cells were mainly enriched in processes related to extracellular matrix organization, cell chemotaxis and the immune response. The DEGs of chondrocytes were enriched in extracellular matrix organization, regulation of vasculature development and immune responses. The DEGs of ECs were enriched in extracellular matrix organization, the regulation of angiogenesis and cell migration. The DEGs of pericytes were enriched in extracellular matrix organization, ossification and aging. The DEGs of stromal cells were enriched in extracellular matrix organization, immune response, and cell adhesion ( Figure 2A ). Comprehensive analysis showed that the five cell types in degenerative NP were mainly enriched in extracellular matrix organization, the immune response, and ossification. Pathway analysis showed that the DEGs of blood cells were enriched in rheumatoid arthritis, phagosomes, and osteoclast differentiation. The DEGs of chondrocytes were enriched in rheumatoid arthritis, phagosome, and antigen processing and presentation. The DEGs of ECs were enriched in ECM-receptor interactions, rheumatoid arthritis, and antigen processing and presentation. The DEGs of pericytes were enriched in rheumatoid arthritis, ECM-receptor interaction, antigen processing and presentation. The DEGs of stromal cells were enriched in rheumatoid arthritis, ECM-receptor interactions, and osteoclast differentiation ( Figure 2B ). Comprehensive analysis showed that the five cell types in degenerative NP were mainly enriched in rheumatoid arthritis, antigen processing and presentation, and the ECM-receptor interaction pathway. As healthy intervertebral discs are immune-silenced, the appearance of blood cells in degenerative NP is an important event associated with IDD (5). Therefore, the activated status of pathways, including rheumatoid arthritis signaling and antigen processing and presentation signaling, in the blood cell cluster was further examined by gene set enrichment analysis (GSEA), which is an unbiased analysis that is independent of differentially expressed genes. The results showed that rheumatoid arthritis signaling and antigen processing and presentation signaling were activated in the blood cell cluster ( Figure 2C ). These results indicated that the immune response was activated in degenerative NP and may be associated with ECM remodeling, further inducing IDD.

To further investigate the effect of the immune response on IDD pathogenesis, we performed blood cell subpopulation analysis and revealed that blood cells could be further divided into MCs (CXCL8, IL-1B, and CCL3) (18), T cells (PTPRC, MT1X and CST7) (19) and CMPs (TPSB2, PRG2, TPSAB1) ( Figures 3A, C ). Statistical analysis of the proportion of blood cell subpopulations showed that the percentage of MCs was the highest, which was consistent with previous studies ( Figure 3B ) (20). Since CMPs accounted for 0.73% of all cells, we mainly focused on analyzing MCs and T cells. Interestingly, T cells were not located in normal NP, and this analysis showed the accumulation of T cells in degenerative NP; thus, we further examined blood samples from IDD patients. CD45, which is encoded by the PTPRC gene, is primarily expressed on the surface of T and B lymphocytes. Immunoassays showed that the proportion of CD45RO+ T cells (associated with activated and differentiated lymphocyte phenotypes) was significantly increased in the blood of LDH patients compared with that of healthy volunteers ( Figure 3D ), which indicated that CD45RO+ T cells may be an immunotherapy target in the future.

Stromal cells are the progenitor cells in the NP region and can contribute to regeneration of the adult intervertebral disc. Additionally, this cell type accounted for a high proportion in the scRNA-seq results. Therefore, we further classified stromal cells into 0-8 subpopulations based on their percentages ( Figures 4A, B ). PDGFRA and PROCR were universally expressed in all subpopulations, similar to the results of a previous study ( Figure 4C ) (3). We further analyzed the biological process and pathway enrichment in the subpopulations of stromal cells. The DEGs of Clusters 0 and 2 were enriched in connective tissue development, cartilage development, and extracellular matrix organization. The DEGs of Cluster 1 were enriched in the response to lipopolysaccharide, bacterial origin and interleukin-1. The DEGs of Clusters 3, 4 and 5 were enriched in extracellular matrix organization, ossification and collagen metabolic processes. The DEGs of Clusters 6 and 7 were enriched in extracellular matrix organization, endodermal cell differentiation, and cell−matrix adhesion. The DEGs of Cluster 8 were enriched in neutrophil degranulation, neutrophil activation involved in the immune response, and extracellular matrix organization ( Figure 5A ). Pathway analysis showed that the DEGs of Clusters 0 and 2 were enriched in osteoclast differentiation, pathogenic Escherichia coli infection, complement and coagulation cascades and phagosomes. The DEGs of Cluster 1 were enriched in cytokine−cytokine receptor interaction, rheumatoid arthritis and osteoclast differentiation. The DEGs of Clusters 3, 4, 5, 6, and 7 were enriched in ECM-receptor interactions, focal adhesions, protein digestion and absorption. The DEGs of Cluster 8 were enriched in rheumatoid arthritis, antigen processing and presentation, ECM-receptor interactions and osteoclast differentiation ( Figure 5B ). These results indicated that the functions of each subpopulation were different and that Subpopulations 3, 4, and 5 were mainly associated with NP ossification and extracellular matrix organization in degenerative NP. Additionally, some novel markers were identified in the subpopulations of stromal cells ( Supplemental Figure S2 ). These findings provide references for investigating the role of stromal cell subpopulations in the NP.

To further investigate the differentiation of stromal cells and corresponding gene expression, we conducted a velocity analysis of stromal cells. Subpopulation 3 was present at the start in the stromal cell differential trajectory ( Figure 6A ). Subpopulation 3 could differentiate into Subpopulation 5 and Subpopulation 8. Additionally, Subpopulation 3 could differentiate into Subpopulations 1, 2 and 4. These findings suggest that Subpopulation 3 possesses progenitor properties. The heatmap shows the specific genes at the differentiation stages of stromal cells. ZFP36L2, ADAM12, FHL1, FGF1, BTG2, FGF2, JARID2, LY6K, FOS, CDK6, ARID5B, GPRC5A, TWIST1, TCF4, GREM1, FOSB, CSF1, DKK3, SMOC1, and FGF18 were associated with cell proliferation, differentiation, development and regeneration. EFNB2, FAT1, ENAH, and POSTN were associated with cell connection, adhesion and motility. MMP2, MMP3, MMP10, MMP13, CEMIP, CEMIP2, COL3A1, COL5A1, COL5A2, COL6A3, COL12A1, COL14A1, COL15A1, TNFAIP6, ADAMTS12, ADAMTS1, ADAMTS5, ADAMTS14, THBS1, and THBS2 were associated with the extracellular matrix and extracellular matrix degradation ( Figure 6A ). Overall, we believed that stromal cells transferred from the characteristics of cell proliferation to the characteristics of cell adhesion and motility and then participated in the regulation of the extracellular matrix, which further resulted in the appearance of ossification characteristics. PDGFRA and THY1 (CD90) are signatures of progenitor cells that were notably expressed in Subpopulation 3; therefore, we further confirmed Subpopulation 3 as the progenitor cells in degenerative NP ( Figure 6B ). Importantly, ossification marker genes, such as OGN, CHSY1, ATP2B1, SNORC, ATF4, BMP2, COL11A1, ENPP1, SPP1, TNFRSF11B, COL1A1, IBSP, RUNX2, and LGALS3, were significantly expressed throughout the stromal cell differentiation process ( Figures 6A, B ) (21, 22). Activation of these ossification genes may promote stromal cell differentiation toward an ossification phenotype, further interfering with the repair functions of stromal cells.

As shown in the stromal cell subpopulation analysis, Clusters 3, 4, and 5 were enriched in extracellular matrix organization, ossification and collagen metabolic processes. The differentiation of stromal cells may be affected by immune cells, and so Cellphone DB analysis was performed. The results showed that MCs and T cells had more interactions with Clusters 3, 4, and 5 than with the other clusters ( Supplementary File 1 ). Therefore, we focused on the effect of immune cells (including macrophages and T cells) on these clusters. We found strong ligand−receptor interactions between MCs and T cells with stromal cells, including IL1B_ADRB2, CXCL2_DPP4, CCL3L1_DPP4, HLA-DRB1_OGN, TNF_RIPK1, TNF_TNFRSF1B, TNF_FAS, TNF_TNFRSF1A, TNF_SEMA4C, TNF_NOTCH1, TNF_DAG1, CCL3_IDE, CCL4_SLC7A1, CD74_APP, CD44_HBEGF, CD74_COPA, CD74_MIF, HLA-C_FAM3C, SPP1_CD44, and HLA-DAP1_TNFSF9. Further analysis revealed that these interactions were mainly associated with cell binding and adhesion, cell growth, and signal transduction and affected the biological processes of stromal cells, including ossification, cell proliferation and differentiation ( Figure 7A , Supplementary File 2 ). Among these interactions, TNF-α, CCL-3 and CD74 frequently acted as ligands, suggesting that TNF-α, CCL-3 and CD74 may play key roles in inducing the ossification of stromal cells. Therefore, immunofluorescence analysis was further performed for validation in vitro. NP cells were isolated from human tissues and treated with recombinant TNF-α, CD74 and CCL-3 proteins, and the expression of ossification markers, including COL1A1 and RUNX2, was measured. The results showed that these factors successfully induced the expression of COL1A1 and RUNX2 ( Figures 7B–D ). As shown previously, TNF-α interacted with TNFRSF1B and TNFRSF1A, which are receptors of the MAPK signaling pathway. To further clarify the signaling pathways in stromal cell ossification that are induced by immune factors, we focused on the effect of TNF-α on the MAPK signaling pathway. The results showed that the upregulated expression of ossification genes (COL1A1 and RUNX2) was significantly reversed after the administration of a MAPK inhibitor (U0126: HY-12031A), which indicated that TNF-α could induce NP cell ossification by activating the MAPK signaling pathway ( Figure 7E ). Additionally, we found strong ligand−receptor interactions between stromal cells, MCs and T cells, including CD99_PILRA MIF_TNFRSF14, MDK_SORL1, FAM3C_CLEC2D, CD44_HBEGF, TNFRSF1B_GRN, TNFRSF1A _GRN, PGRMC2_CCL4L2, and the THY1_axb2 complex, which may change immune cell functions such as cell proliferation, apoptosis and differentiation.

Furthermore, we examined the role of inhibiting immune factors, such as TNF-α, on NP ossification and IDD using a rat coccyx disc degeneration model. First, narrowing of the intervertebral space was found by X-ray in the degeneration group ( Supplemental Figure S3A ). Additionally, MRI showed that the NP signal switched to gray with decreased disc height in the degeneration group ( Supplemental Figure S3B ). HE staining and safranin O-fast green staining indicated that the NP had a long spindle shape with a decrease in cell number and disordered collagen fibers. The boundary between the AF and NP became blurred in the degeneration group ( Supplemental Figures S3C, D ). These results demonstrated the successful establishment of the disc degeneration animal model. In this animal model, immunohistochemical staining showed significantly increased expression of TNF-α in the AFP group ( Figure 8A ). Blocking TNF-α with a specific inhibitor decreased the expression of ossification markers, including Col-1 and Runx-2 ( Figure 8B ), and improved disc morphology ( Figures 8C, D ). Furthermore, X-ray and MRI showed improvements in disc imaging ( Figures 8E, F ). Therefore, we clearly indicated that inhibiting TNF-α in the NP could alleviate intervertebral disc ossification and degeneration.