Single-Cell Transcriptome Profiling Reveals Multicellular Ecosystem of Nucleus Pulposus during Degeneration Progression.

Although degeneration of the nucleus pulposus (NP) is a major contributor to intervertebral disc degeneration (IVDD) and low back pain, the underlying molecular complexity and cellular heterogeneity remain poorly understood. Here, a comprehensive single-cell resolution transcript landscape of human NP is reported. Six novel human NP cells (NPCs) populations are identified by their distinct molecular signatures. The potential functional differences among NPC subpopulations are analyzed. Predictive transcripts, transcriptional factors, and signal pathways with respect to degeneration grades are explored. It is reported that fibroNPCs is the subpopulation for end-stage degeneration. CD90+NPCs are observed to be progenitor cells in degenerative NP tissues. NP-infiltrating immune cells comprise a previously unrecognized diversity of cell types, including granulocytic myeloid-derived suppressor cells (G-MDSCs). Integrin αM (CD11b) and oxidized low density lipoprotein receptor 1 (OLR1) as surface markers of NP-derived G-MDSCs are uncovered. The G-MDSCs are found to be enriched in mildly degenerated (grade II and III) NP tissues compared to severely degenerated (grade IV and V) NP tissues. Their immunosuppressive function and alleviation effects on NPCs' matrix degradation are revealed in vitro. Collectively, this study reveals the NPC-type complexity and phenotypic characteristics in NP, thereby providing new insights and clues for IVDD treatment.

Introduction

Low back pain is a major disabling health condition in humans, with a lifetime prevalence of as high as 84%.[ ] The socioeconomic burden of this rheumatologic disorder of the spine is enormous. It was estimated at $85 billion in 2008, and the economic cost has still increased in the last decade.[ ] Intervertebral disc degeneration (IVDD) is a widely recognized contributor to low back pain.[ ] The degeneration of the nucleus pulposus (NP), the central gel‐like part of the intervertebral disc, is a significant mechanism of IVDD.[ ] The current treatments are limited to relieving back or leg symptoms. They do not focus on replenishing the NP loss and restoring the native disk structure. The failure leads to unsatisfactory outcomes, such as recrudescence or degeneration of adjacent motion segments.[ ] New therapeutic targets are therefore needed.

NP cells (NPCs) are the main cell type residing in the NP, and they are responsible for maintaining tissue homeostasis.[ ] NPCs proliferate slowly and lack self‐regeneration capacity adding to the intractability of the disease. Current studies on the pathophysiology of NPCs are usually supported by transcriptomic and epigenomic analyses. However, bulk‐tissue level resolution masks the complexity of alterations across cells and within cell types. The uncharacterized cell types and markers residing in the NP raise interest in terms of unexplored cellular heterogeneity.

IVD is the largest avascular organ of the body. Studies addressing immune reactions against the NP have been limited and focused on their detrimental aspects. However, given the complexity of immunity in other immune‐privileged sites,[ ] immune cell subpopulations may also help restore IVD structure and lessen degeneration. Moreover, intrinsic properties of NPCs, including expression of immunomodulatory factors, and extrinsic microenvironmental changes to immune compartments, remain largely unknown. These issues highlight the importance of understanding the immune panorama of the NP during IVDD pathogenesis.

Single‐cell RNA sequencing (scRNA‐seq) provides a powerful alternative to study the cellular heterogeneity of NP tissues. Here, we aimed to provide a single‐cell view of IVDD pathology, profiling 39 732 cells from NP tissues across eight individuals with different grades of progressive degeneration. Notably, we comprehensively characterized the transcriptome feature of NPCs and immune cells, and we decoded the cell percentage, the heterogeneity of cell subtypes during degeneration, providing a unique cellular‐level insight into transcriptional alterations associated with IVDD pathology.

Results

Single‐Cell Profiling of NPC Atlas in Human Subjects with IVDD Pathology

We dissociated NP tissues from eight IVDD patients (Table ) with different degrees of degeneration (see the Experimental Section) and performed scRNA‐seq on the BD Rhapsody system (Figure and Figure S1A, Supporting Information). After quality control and doublet exclusion filtering to remove cells with low gene detection (<600 genes) and high mitochondrial gene content (>8%) (Figure S1E–H, Supporting Information), NPCs were identified based on their levels of the transcripts that encode different proteins [e.g., aggrecan proteoglycan (ACAN) and SRY‐box transcription factor 9 (SOX9)]. Consistent with previous studies,[ ] we found a lack of distinctive clusters in the tSNE map (Figure 1B) which may suggest the heterogeneity exhibited from these cell populations. Six subpopulations were identified based on their highly expressed genes and published single‐cell studies (Figure 1C):

Table 1

Basic information and characteristics for participants

Patient ID	Age [years]	Gender	Weight [kg]	Reason for surgery	Pfirrmann grading	CRP [mg dL−1]	WBC (109 L−1)	Lymphocyte (109 L−1)	Monocyte (109 L−1)
S1	63	Male	60	Burst fracture	II	0.047	6.46	2.16	0.64
S2	41	Male	73.5	Burst fracture	II	0.706	11.37	2.35	1.02
S3	56	Female	62	Lumbar disc herniation	III	0.079	7.25	1.63	0.47
S4	65	Female	76	Lumbar disc herniation	III	0.154	6.78	2.15	0.54
S5	64	Female	50	Lumbar disc herniation	IV	0.027	4.83	1.55	0.19
S6	53	Female	60	Lumbar disc herniation	IV	0.184	5.53	2.28	0.49
S7	54	Male	68	Lumbar disc herniation	V	0.011	10.79	3.55	0.83
S8	56	Male	55	Lumbar disc herniation	V	0.022	5.75	1.86	0.61

Figure 1

Identification of human NPC atlas and transcriptional changes correlated with IVDD severity. A) Graphical representation of the experimental workflow. B) UMAP visualization of human NP cells identified six different clusters after unsupervised clustering. Each dot corresponds to one single cell colored according to cell cluster. C) Heatmap revealing the scaled expression of differentially expressed genes for each cluster. D) Dot plots showing the grade distribution in each NP cell subsets. E) RT‐qPCR for the representative genes of NPC atlas in different degenerative grades discs (n = 3 with mean ± SD shown). F) Representative immunohistochemistry assay of indicated genes in NP tissues. G) Heatmap showing grade‐related transcription factors. H) Enriched signal pathways related with degeneration grades.

hypertrophy chondrocyte‐like NPCs (HT‐CLNPs; Groups 0, 1; expressing FRZB, DKK);[ ]

effector NPCs (Groups 2 and 4, expressing mRNAs that encode proteins that participate in genes cellular metabolism, e.g., MSMO1[ ] and HMGCS1[ ]);

homeostatic NPCs (Group 5; expressing RPS29 and RPS21);[ ]

regulatory NPCs (Group 6; expressing CHI3L1,[ ] CXCL2, and NFKB[ ]);

fibroNPCs (Groups 8, 9, and 10, expressing mRNAs that encode related with fibrosis, COL1A1, COL3A1, and COL6A1); and[ ]

adhesion NPCs (Groups 3 and 11; expressing mRNAs that related to cell adhesion and migration, such as FN1[ ] and CRTAC1[ ]).

We next analyzed the relationship between the degeneration grades and distributions of the cell populations. HT‐CLNPs and regulatory NPCs were the main cells for grade II discs. Effector NPCs and HT‐CLNPs were major subpopulations in grades III and IV discs. For grade V discs, fibroNPCs and adhesion NPs were the main NPCs (Figure 1D). Moreover, the proportion of adhesion NPCs and fibroNPCs increased with the severity of IVDD (from grade II to V), while the proportion of homeostatic NPCs showed opposite trends.

Real‐time quantitative PCR (qPCR) was used to validate the indicated gene expression in NPCs from different degeneration grades. The results revealed that the expressions of FN1 and CRTAC1, COL3A1, and MMP2 (markers of adhesion NPCs and fibroNPCs) were significantly elevated in NPCs from late‐stage degenerative discs (grade IV and V). The expressions of homeostatic NPC markers, MSMO1 and HMGS1, significantly decreased in NPCs from late‐stage degenerative discs. (Figure 1E). We then performed an immunohistochemistry assay to validate the expressions of markers for each NPC subpopulation at the protein level (Figure 1F).

Identification of NPC Transcriptional Changes Correlated with IVDD Severity

A group of transcriptional factors (TFs) were identified to be related to grade of degeneration (Figure 1G). Some TFs, like AR, REL, LMX1A, and PRDM1, were first revealed. Enrichment for TFs of REL in grade IV and V discs may be related to Rel/NF‐κB signal transduction pathway in IVDD.[ ]

Besides, we compared differentially expressed genes (DEGs) between severe (grades IV and V) and mild (grades II and III) degeneration and identified candidate genes related to IVDD progression, including HSPH1, CTGF, MMP13, and HSP90AA1 (Table ). There are significantly elevated expressions of a list of heat shock protein (HSP) genes in severe degenerative discs, including HSPA1B, HSPH1, HSP90AA1, HSPA8, HSPA1A, HSPB8, and HSPD1, indicating the importance of stress‐related mechanisms in IVDD. Enrichment for TF, SOX4, and INHBA genes is consistent with the importance of TGF‐β signaling in IVDD.[ ] ANXA5 was also elevated, which may be related to mitochondrial dysfunction induced cell apoptosis.[ ]

Table 2

DEGs by comparing severe degeneration grade (Grade III and IV) to mild degeneration grade (Grade I and II)

Upregulated genes
Gene	P value	avg_logFC	pct.1	pct.2	p_val_adj
INHBA	0	1.542242	0.787	0.491	0
FHL2	0	1.355454	0.65	0.274	0
KLHL21	0	1.150942	0.633	0.334	0
CYR61	0	1.149118	0.947	0.887	0
HSPA1B	0	1.059508	0.868	0.733	0
NFATC2	0	1.027662	0.497	0.173	0
HSPH1	0	0.983483	0.813	0.634	0
CTGF	0	0.979719	0.957	0.877	0
MMP13	0	0.948365	0.822	0.626	0
HSP90AA1	0	0.899521	0.964	0.896	0
RTN4	0	0.759251	0.836	0.694	0
HSPA8	0	0.650494	0.882	0.8	0
ACTG1	0	0.627876	0.979	0.96	0
ANXA5	0	0.608453	0.908	0.839	0
LUM	0	0.58595	0.984	0.97	0
H3F3B	0	0.496407	0.963	0.94	0
HSPA1A	5.2 × 10−284	0.735349	0.91	0.866	1.1 × 10−279
GFPT2	6.9 × 10−261	0.734121	0.682	0.489	1.5 × 10−256
KLF4	7 × 10−258	0.702163	0.905	0.846	1.5 × 10−253
ATP1B1	6.7 × 10−252	0.710599	0.614	0.392	1.5 × 10−247
PCOLCE2	1.2 × 10−250	0.563587	0.869	0.752	2.6 × 10−246
EMP1	7.9 × 10−244	0.770794	0.907	0.84	1.7 × 10−239
SLC39A14	9.2 × 10−242	0.606197	0.801	0.68	2 × 10−236
HERPUD1	1.3 × 10−237	0.515521	0.873	0.801	2.8 × 10−233
PLOD2	5.6 × 10−227	0.471544	0.866	0.805	1.2 × 10−222
HSPB8	1.8 × 10−224	0.780165	0.687	0.502	4 × 10−220
ZSWIM6	1.8 × 10−219	0.548961	0.309	0.088	3.8 × 10−215
BTF3	3.4 × 10−219	0.457349	0.863	0.785	7.4 × 10−215
ANXA1	2.2 × 10−216	0.58282	0.928	0.894	4.7 × 10−212
ACKR3	4.6 × 10−215	0.909398	0.595	0.397	1 × 10−210
MYADM	1.2 × 10−214	0.682027	0.739	0.614	2.5 × 10−210
ACTB	1.6 × 10−214	0.357681	0.965	0.945	3.4 × 10−210
AMOTL2	1.1 × 10−212	0.883341	0.388	0.169	2.5 × 10−208
FOXC2	1.9 × 10−210	0.758533	0.43	0.208	4.2 × 10−206
TXN	1.2 × 10−201	0.630013	0.708	0.571	2.6 × 10−197
UAP1	7.7 × 10−201	0.602079	0.708	0.556	1.7 × 10−196
TM4SF1	4.7 × 10−200	0.936974	0.488	0.282	1 × 10−195
RYBP	1.2 × 10−199	0.648708	0.654	0.479	2.7 × 10−195
FMOD	1.1 × 10−195	0.535043	0.908	0.89	2.3 × 10−191
RPL30	1.8 × 10−186	0.341223	0.949	0.941	3.9 × 10−182
GPRC5A	7.1 × 10−186	0.670452	0.596	0.403	1.5 × 10−181
HSPD1	1.2 × 10−184	0.574874	0.775	0.666	2.7 × 10−180
FN1	4.3 × 10−180	0.740568	0.968	0.95	9.4 × 10−176
TMED2	2.4 × 10−177	0.583225	0.672	0.564	5.2 × 10−173
RAN	1.5 × 10−174	0.462583	0.799	0.715	3.3 × 10−170
SERPINE2	2 × 10−174	0.721224	0.845	0.773	4.4 × 10−170
CTNNB1	1.9 × 10−167	0.511848	0.66	0.526	4.2 × 10−163
IL11	7.6 × 10−165	0.928692	0.287	0.1	1.6 × 10−160
OAT	1 × 10−163	0.473536	0.728	0.608	2.3 × 10−159
CD55	6.3 × 10−161	0.717146	0.707	0.576	1.4 × 10−156

Based on a IVDD grade‐related gene set, a group of signaling pathways were found to potentially promote degeneration, including antigen processing and presentation, the TNF pathway, MAPK, and Hippo pathways (Figure 1H).

ScRNA‐seq Reveals Transcriptional Features of NPCs Subpopulations

To analyze the functional differences among subpopulations, we used the Quantitative Set Analysis for Gene Expression (QuSAGE) (Figure ) and Gene Ontology (GO) database (Table )[ ] to investigate the biological process of our identified NPCs subpopulations. Effector NPCs were enriched with metabolic process‐related genes (e.g., sterol biosynthesis and glycosaminoglycan metabolism) and positive regulation of extracellular matrix (ECM) assembly. Regulatory NPCs were enriched with highly expressed genes responsible for cellular responses to inflammation and endogenous stimuli, such as CXCL3, IL‐6, and CHI3L1, indicating that these cells might potentially regulate immune functions. Homeostatic NPCs were enriched for processes related to cellular homeostasis, including translational regulation and protein/RNA metabolism. Adhesion NPCs were enriched for cell migration and cell–matrix adhesion.

Figure 2

ScRNA‐seq reveals transcriptional features of NPCs subpopulations. A) QuSAGE analysis of cell subpopulation specific differential expression colored by statistically significant normalized enrichment scores. B–E) Violin plots of steroid biosynthesis, immune response, innervation, and angiogenesis score for each cluster. F) Correlation of scRNA‐seq defined NPCs subpopulations with cell senescence. G) Heatmap showing the scaled expression of the differentially expressed genes (DEGs) for HT‐CLNP‐I and HT‐CLNP‐II subsets. H) Pseudotime trajectory axis revealing the progression of HT‐CLNP‐I and HT‐CLNP‐II.

Table 3

GO Terms for NP subpopulations

Efector NPCs GO terms
GOID	GOTerm	P‐value
GO:0016126	sterol biosynthetic process	2.432 × 10−11
GO:0006695	cholesterol biosynthetic process	5.795 × 10−10
GO:0030198	extracellular matrix organization	3.136 × 10−8
GO:0044281	small molecule metabolic process	4.528 × 10−8
GO:0001558	regulation of cell growth	2.285 × 10−6
GO:0001501	skeletal system development	3.133 × 10−6
GO:0006694	steroid biosynthetic process	5.308 × 10−6
GO:0008203	cholesterol metabolic process	9.283 × 10−6
GO:0008202	steroid metabolic process	1.061 × 10−5
GO:0006048	UDP‐N‐acetylglucosamine biosynthetic process	1.904 × 10−5
GO:0008299	isoprenoid biosynthetic process	2.17 × 10−5
GO:0010628	positive regulation of gene expression	3.007 × 10−5
GO:0008285	negative regulation of cell proliferation	3.96 × 10−5
GO:0009405	pathogenesis	5.295 × 10−5
GO:0072593	reactive oxygen species metabolic process	6.08 × 10−5
GO:0006629	lipid metabolic process	9.184 × 10−5
GO:0048661	positive regulation of smooth muscle cell proliferation	9.297 × 10−5
GO:0033173	calcineurin‐NFAT signaling cascade	0.0001037
GO:0030199	collagen fibril organization	0.0001053
GO:0042127	regulation of cell proliferation	0.0001366
GO:0001525	angiogenesis	0.0001439
GO:0014070	response to organic cyclic compound	0.0001497
GO:0010955	negative regulation of protein processing	0.0001541
GO:0042340	keratan sulfate catabolic process	0.0001541
GO:0006011	UDP‐glucose metabolic process	0.0001563
GO:0007568	aging	0.0001871
GO:0006633	fatty acid biosynthetic process	0.0002058
GO:0043434	response to peptide hormone	0.000223
GO:0055114	oxidation‐reduction process	0.0002321
GO:0005975	carbohydrate metabolic process	0.0003237
GO:0043065	positive regulation of apoptotic process	0.0003622
GO:0030203	glycosaminoglycan metabolic process	0.0003728
GO:0043066	negative regulation of apoptotic process	0.000445
GO:0003331	positive regulation of extracellular matrix constituent secretion	0.0004652
GO:0006065	UDP‐glucuronate biosynthetic process	0.0004652
GO:2000544	regulation of endothelial cell chemotaxis to fibroblast growth factor	0.0004652
GO:0042542	response to hydrogen peroxide	0.0004732
GO:0032270	positive regulation of cellular protein metabolic process	0.0005055
GO:0006915	apoptotic process	0.0006357
GO:0010035	response to inorganic substance	0.0006375
GO:0008284	positive regulation of cell proliferation	0.0006802
GO:0044255	cellular lipid metabolic process	0.0007222
GO:0060591	chondroblast differentiation	0.0009226
GO:1903169	regulation of calcium ion transmembrane transport	0.0009226
GO:0006933	negative regulation of cell adhesion involved in substrate‐bound cell migration	0.0009226
GO:0051387	negative regulation of neurotrophin TRK receptor signaling pathway	0.0009226
GO:0009314	response to radiation	0.0009242
GO:0060291	long‐term synaptic potentiation	0.0010297
GO:0045444	fat cell differentiation	0.0013192
GO:0031394	positive regulation of prostaglandin biosynthetic process	0.001525

We specifically scored for innervation, angiogenesis, and cell senescence, which are important events related to IVDD,[ ] based on GO terms and published literature (see the Experimental Section). It was suggested that there is no significant difference for innervation among NPC subpopulations (Figure 2D). For angiogenesis, fibroNPCs showed the highest score, and homeostatic NPCs showed the lowest score. (Figure 2E) Via SCENIC analysis, we identified the specific signature transcription factor motifs for each subpopulation, and their predictive transcriptional control (Figure S2B, Supporting Information).

Also, it was indicated that senescence‐associated secretory phenotype (SASP) was almost not expressed in HT‐CLNPs or homeostatic NPCs. However, HT‐CLNPs were active for cell cycle arrest (DNA damage) and SASP‐related epigenetics changes. FibroNPCs were highly active for macromolecular damage, indicating that HT‐CLNPs showed different cellular senescence mechanisms to FNPs (Figure 2F).

Overall, our findings suggest differences in functions and biological processes among NPC subpopulations. Moreover, multiple potential mechanisms for cell senescence exist in different subpopulations.

Two transcriptionally distinct clusters were identified: HT‐CLNP‐I and HT‐CLNP‐II. The DEGs of these two clusters are shown in Figure 2G. The HT‐CLNP‐I was enriched in genes related to programmed cell death and negative regulation of cell proliferation. In contrast, the HT‐CLNP‐II subset was enriched in genes related to ECM organization and ECM disassembly (Table ). Additionally, the Monocle pseudotime trajectory revealed progression of the HT‐CLNP‐I and HT‐CLNP‐II clusters (Figure 2H). These findings contribute to a deeper understanding of the hydrophilic mechanisms in NPCs during IVDD pathogenesis.

Table 4

Comparison of HT‐CLNPs‐I and HT‐CLNPs‐II

HT‐CLNPs‐I GO terms
GOID	GOTerm	P‐value
GO:0006048	UDP‐N‐acetylglucosamine biosynthetic process	3.524 × 10−7
GO:0012501	programmed cell death	1.952 × 10−6
GO:0008285	negative regulation of cell proliferation	2.502 × 10−6
GO:0042026	protein refolding	3.032 × 10−6
GO:0010467	gene expression	1.394 × 10−5
GO:0043066	negative regulation of apoptotic process	1.964 × 10−5
GO:0031397	negative regulation of protein ubiquitination	2.025 × 10−5
GO:0010941	regulation of cell death	2.151 × 10−5
GO:0007596	blood coagulation	2.385 × 10−5
GO:1900740	positive regulation of protein insertion into mitochondrial membrane involved in apoptotic signaling pathway	4.689 × 10−5
GO:0045648	positive regulation of erythrocyte differentiation	4.689 × 10−5
GO:0002042	cell migration involved in sprouting angiogenesis	6.276 × 10−5
GO:1900034	regulation of cellular response to heat	7.543 × 10−5
GO:0006047	UDP‐N‐acetylglucosamine metabolic process	8.473 × 10−5
GO:0048010	vascular endothelial growth factor receptor signaling pathway	9.284 × 10−5
GO:0045944	positive regulation of transcription from RNA polymerase II promoter	0.0001026
GO:0006928	cellular component movement	0.0001158
GO:0006915	apoptotic process	0.0001216
GO:0034605	cellular response to heat	0.0001572
GO:0007623	circadian rhythm	0.0001957
GO:0032092	positive regulation of protein binding	0.0002012
GO:0042267	natural killer cell mediated cytotoxicity	0.0002278
GO:0034629	cellular protein complex localization	0.0002315
GO:0033173	calcineurin‐NFAT signaling cascade	0.0002315
GO:0045893	positive regulation of transcription, DNA‐templated	0.0002562
GO:0007264	small GTPase mediated signal transduction	0.0002668
GO:0006986	response to unfolded protein	0.0003037
GO:0032060	bleb assembly	0.000343
GO:0006357	regulation of transcription from RNA polymerase II promoter	0.0003793
GO:0097193	intrinsic apoptotic signaling pathway	0.0004046
GO:0006367	transcription initiation from RNA polymerase II promoter	0.0004339
GO:0050821	protein stabilization	0.0004345
GO:0010628	positive regulation of gene expression	0.0004551
GO:0061045	negative regulation of wound healing	0.0004841
GO:0070374	positive regulation of ERK1 and ERK2 cascade	0.000569
GO:0051085	chaperone mediated protein folding requiring cofactor	0.0006576
GO:0006366	transcription from RNA polymerase II promoter	0.0006606
GO:0045597	positive regulation of cell differentiation	0.0007432
GO:0060128	corticotropin hormone secreting cell differentiation	0.0007994
GO:0070370	cellular heat acclimation	0.0007994
GO:2000544	regulation of endothelial cell chemotaxis to fibroblast growth factor	0.0007994
GO:0043123	positive regulation of I‐kappaB kinase/NF‐kappaB signaling	0.0008046
GO:0033138	positive regulation of peptidyl‐serine phosphorylation	0.0010132
GO:0060129	thyroid‐stimulating hormone‐secreting cell differentiation	0.0015814
GO:0060591	chondroblast differentiation	0.0015814
GO:0010664	negative regulation of striated muscle cell apoptotic process	0.0015814
GO:0022614	membrane to membrane docking	0.0015814
GO:0015936	coenzyme A metabolic process	0.0015814
GO:0070434	positive regulation of nucleotide‐binding oligomerization domain containing 2 signaling pathway	0.0015814
GO:0045765	regulation of angiogenesis	0.0015841

CD90+NPCs Is the Progenitor within FibroNPCs, the End‐Stage Subpopulation

Next, we investigate the relationship between the different cell clusters. Cellular (Cyto) Trajectory Reconstruction Analysis using gene Counts and Expression (CytoTRACE) can accurately uncover the direction of differentiation and predict cell lineage trajectories.[ ] The results indicated the order of the differentiation states as fibroNPCs, adhesion NPCs, effector NPCs, regulatory NPCs, HT‐CLNPs, and homeostatic NPCs (Figure ).

Figure 3

CD90+NPCs is the progenitor within FibroNPCs, the end‐stage subpopulation. A) Plot of the cytoTRACE pseudotime order for the NP subpopulations. The value of cytoTRACE represents the predicted order. B) Visualization for dynamic velocities projected into the UMAP‐based embedding. C) The expression of CD90, CD44, CD73, and CD29 in NPCs, the red box represents the region of fibroNPCs on UMAP. D) Histogram to evaluate the relative expression of CD90 after cell sorting. E) Left: oil red staining for CD90+ NPC‐induced adipogenic differentiation, respectively (n = 3). Scale bar, 100 µm. Right: alizarin red staining for CD90+ NPC‐induced osteogenic differentiation (n = 3). Scale bar, 100 µm. F) RT‐qPCR of fibroNPCs phenotype mRNA levels between CD90+/− NPCs. (n = 3 with mean ± SD shown, *P < 0.05, **P < 0.005). G) RT‐qPCR of adhesion NPCs phenotype at different time points after cultured CD90+NPCs in chondrogenesis induced medium (n = 3 with mean ± SD shown). H) Safranin O staining of CD90+NPCs from mild (grade II and III) or severely (grade IV and V) degenerative individuals after culturing in chondrogenesis medium for 21 d. I) Correlation of scRNA‐seq defined NPC subpopulations with cell death and inflammasome. J) Western blot analysis with representative blots including Bax, Bcl‐2m NLRP3 levels in the CD90+/− NPCs. Densitometric analysis is shown as mean ± SD, n = 3; *P < 0.05, **P < 0.005. K) Immunofluorescence (IF) visualization of CD90 (red) and nuclei (blue) in degenerative disc tissues induced in rats.

In addition, we applied the scVelo trajectory algorithms, a novel method developed for recovering the position of each cell in the underlying differentiation processed based on inferred gene‐specific rates of transcription, splicing, and degradation.[ ] The arrows in scVelo indicate the estimated sequence of transcriptomic events (Figure 3B). The results predict that the fibrosis NPCs may be differentiated into adhesion NPCs and other NPCs subsets. Monocle analysis also revealed that FNPs existed at the start of the pseudospace trajectory. Adhesion NPCs were distributed along the trajectory, and homeostatic NPCs were mainly distributed at the end (Figure S3A, Supporting Information). Together, these findings imply that fibroNPCs possess progenitor properties.

The scVelo analysis helped us to systematically identify putative driver genes as genes characterized by high likelihoods in fibroNPC populations (Figure S3B, Supporting Information). In other words, these genes may work as candidates for important drivers of the main process in firboNPCs. These genes have been associated with matrix remodeling (CoL14A1, coL12a1, CHAD, CRTAC1, TNC, Lamb1),[ ] antioxidation (FTH1),[ ] and inflammation (S100a8).[ ]

Because fibroNPCs may have progenitor properties, we analyzed the expression of markers of mesenchymal stem cells (MSCs) in our single‐cell data and gated the fibroNPC region with MSCs markers (red box in Figure 3C) on visualized UMAPs. The results show that CD90 was expressed mainly in the fibroNPC region compared to the other MSC markers, CD44, CD73, and CD29. We, therefore, hypothesize that CD90+ NPCs may be the progenitor cells in degenerative NP tissues.

To confirm whether CD90+ NPCs are the progenitor cells in NP, we isolated CD90+NPCs from human degenerated NPCs using microbeads (Figure 3D). The CD90+NPCs were positive for CD44 and CD29 and negative for CD34 and HLA‐DR (Figure S5A, Supporting Information). Oil red staining and Alizarin red staining revealed that CD90+ cells had multipotent capabilities (differentiated into various cell lineages, including osteoblasts and adipocytes) (Figure 3E). RT‐qPCR showed that fibroNPCs phenotype genes (including CLO1A1, COL3A1, MMP2) were significantly elevated in CD90+NPCs compared to CD90‐NPCs, however, the expression of TGFb β showed no significant difference between two groups (Figure 3F). We also cultured CD90+NPCs in the chondrogenesis differentiation medium and found that adhesion NPCs genes, including FN1, CRTAC1, and FMOD, were elevated from day 1 and peaked at day 6 but decreased after that (Figure 3G). This result indicates that CD90+NPCs can differentiate into cells with adhesion NPCs phenotypes. Safranin O staining showed a higher chondrogenesis potential for CD90+ NPCs isolated from mild degenerative (grade II and III) individuals compared to those from severe degenerative (grade IV and V)individuals (Figure 3H). These findings demonstrate that CD90+NPCs expressed phenotypes of fibroNPCs and can serve as progenitor cells in degenerative NP tissues.

Both cell repair and cell death are involved in tissue regeneration.[ ] To explore whether different cell death events occurred in our predictive NPC subsets, we correlated our single‐cell data with previous publications.[ ] The correlation analysis showed that fibroNPCs were active for all four cell death types: apoptosis, pyroptosis, necroptosis, and ferroptosis. Inflammasome may be highly involved in cell death in fibroNPCs but not in other cell types. Ferroptosis was evident in HT‐CLNPs, adhesion NPCs, and regulatory NPCs. Necroptosis occurred in fibroNPCs, adhesion NPCs, and regulatory NPCs (Figure 3I). Next, we compared apoptosis and NLRP3‐related proteins expressed between CD90+ NPCs and CD90− NPCs; CD90+ cells showed decreased levels of the Bax and NLRP3 transcripts, as well as increased Bcl‐2 expression (Figure 3J). Furthermore, we used immunofluorescence to explore the location of CD90+ NPCs in puncture‐induced degenerative rat IVD.[ ] It was found that CD90+ NPCs were enriched in the boundary between NP and AF; only a few positive cells could be detected in the inner area of the discs (Figure 3K). Therefore, the combined results showing fibroNPCs showed profound cell death activity. Nevertheless, CD90+NPCs may help NP regeneration in this subpopulation. As it was shown in Figure 1D, the proportion of fibroNPCs is elevated with degeneration progression, which indicating fibroNPCs could be the end‐stage subpopulation. Overall, fibroNPCs are the major cell population for end‐stage degenerative NPCs with the persistence of both cell death and progenitors. The CD90+NPCs as the progenitor NPCs can be a potential therapeutic target in the future study.

Identification of NP‐Derived G‐MDSCs

Apart from NPCs, immune cells were identified in NP tissues. Using the canonical correlation analysis (CCA) method,[ ] natural killer (NK) cells were identified by CD94.[ ] Macrophage/monocytes were identified by CD163,[ ] and T cells were identified by T Cell Receptor Alpha Constant (TRAC) (Figure 1B,C). G‐GMPs (expressing MS4A3, MPO, and ELANE),[ ] neutrophils (FCGR3B (CD16b) + HLADR),[ ] and granulocytic myeloid‐derived suppressor cell (G‐MDSC, ITGAM (CD11b), OLR1, and ARG1)[ ] were also identified (Figure ). Monocle trajectory analysis revealed GMPs were distributed at the start of the trajectory, G‐MDSCs at the intermediate segment, and neutrophils were dispersed along the trajectory (Figure 4D).

Figure 4

Identification of NP‐derived G‐MDSCs. A) Uniform manifold approximation and projection (UMAP) visualization showing immune cells inside NP tissues. B) Heatmap showing the typically expressed genes in each cell type. C) Dot plot showing scaled expression of selected signature genes for GMP, G‐MDSCs, and neutrophils, by average expression of each gene in each cluster scaled across all clusters. Dot size represents the percentage of cells in each cluster with more than one read of the corresponding gene. D) Monocle method reconstructed pseudospace trajectory for GMP, G‐MDSC, and neutrophils. E) FACS isolation for CD45+ CD11b+OLR+ cells from mild/severe degenerative NP tisues. F) Safranin O/Fast Green staining of the intervertebral discs sham and experimental rats. Scale bar, 1 mm in left and 100 µm in right. G) Merged immunofluorescence staining of DAPI, CD11b, OLR1 in the intervertebral discs of sham and experimental rat. Scale bar, 1 mm in merge images and 100 µm in others.

Our data show that there is an enrichment of G‐MDSCs in severely degenerated NP tissues (grade IV and V) compared to mildly degenerative tissues (grade II and III) at the single‐cell level. To confirm this finding, we isolated G‐MDSCs from human NP tissues via Fluorescence‐activated cell sorting (FACS) (Figure S4A, Supporting Information). CD45, the marker of hematopoietic cells, CD11b, OLR1 were markers used for G‐MDSCs sorting. G‐MDSCs (CD45+CD11b+OLR1+ cells) decreased by almost threefolds in severe degeneration discs (6.5% vs 2.2%) (Figure 4E). To further validate this finding in vivo, we used rat models of IVDD with the needle puncture method.[ ] Histological staining showed notable changes including cell cloning (cell clusters), loss of demarcation between NP and AF, and proteoglycan loss in punctured IVDs, demonstrating progressive IVDD at two and six weeks (Figure 4F). Immunofluorescence staining for CD11b, OLR, and CD24 showed the existence of G‐MDSCs in rat NPs, and G‐MDSCs decreased in the six‐week group relative to the two‐week group (Figure 4G).

Functional Validation NP‐Derived MDSCs

The most typical and important function of G‐MDSCs is immunosuppression.[ ] The hallmark of this immunosuppressive activity is the capability of suppressing T cell activation and ROS production.[ ] Next, we validate the functions of NP‐derived MDSCs via immunosuppression and ROS production (Figure ).

Figure 5

Validation of functions of NP‐derived G‐MDSCs. A) Schematic workflow of the experimental strategy. B) FACS isolation for CD45+ CD11b+OLR+CD24+ and CD24− cells. C) T cell suppression analysis NP‐derived G‐MDSC identification. N = 3, *P < 0.05, **P < 0.005, ***P <0 .001. D,E) CD11b+OLR+CD24+ show increased reactive oxygen species (ROS) formation compared to CD24− cells. Rosup treated cells were used as positive control. N = 3, *P < 0.05, **P < 0.005. F) RT‐qPCR of degeneration related genes, aggrecan, collagen II, ADAMTS4,5, and MMP13 in untreated, IL‐1b+G‐MDSCs, and IL‐1b alone NPC groups. N = 3, *P < 0.05, **P < 0.005, ***P < 0.001. G) Western blot analysis with representative blots including aggrecan and MMP13 in untreated, IL‐1b+G‐MDSCs, and IL‐1b alone NPC groups. N = 3, *P < 0.05, **P < 0.005, ***P < 0.001.

As we found that CD24 also highly expressed in MDSCs at single‐cell level, we tested this finding by FACS. As it is shown in Figure 5B, the NP‐derived MDSCs were divided into CD24+ MDSCs and CD24‐MDSCs. To functionally validate whether CD45+CD11b+OLR1+CD24+/− cells inhibit immune cell activation, we performed co‐cultures with activated T cells (Figure 5C).[ ] Human NP tissues, both CD45+CD11b+OLR1+CD24+ and CD24−, exhibited T cell suppressive capacity; however, CD24+ cells showed significantly stronger suppressive capacity compared to CD24− cells (Figure 5C). Similar to T‐cell proliferation, interleukin 2 (IL2) production was significantly suppressed by G‐MDSC co‐cultured (Figure S4B, Supporting Information). We also observed that CD45+CD11b+OLR1+CD24+ cells produced significantly higher amounts of ROS compared to CD24− cells (Figure 5D,E). These findings indicate that our identified G‐MDSCs in human NP tissues have T‐cell suppression and ROS production capabilities.

Because G‐MDSCs decreased in severely degenerated NP tissues, we hypothesized that G‐MDSCs might protect against IVDD. To explore the effects of G‐MDSCs on NPCs, G‐MDSCs were isolated from mildly degenerated (grade II and III) human NP tissues and co‐cultured with non‐degenerative NPCs. IL‐1β was used to induce NPC degeneration.[ ] Consequently, compared with IL‐1β treatment alone, G‐MDSCs co‐cultured with NPCs showed increased aggrecan expression and significantly lower MMP13 expression (Figure 5F,G). These in vitro findings indicate that G‐MDSCs could alleviate IL‐1β induced NPCs matrix degradation.

Characterization of Cell‐to‐Cell Interactions Involved in IVDD

Our study showed that NP tissues consist of NPCs and immune cells. We next predicted the cell–cell interaction network among the cell types using CellPhoneDB 2.0.[ ] Macrophages showed the most interactions with other cell types. Interactions among NPCs were the most active (Figure ).

Figure 6

Predicted immune‐NPCs regulatory network in IVDD. A) Heatmap showing the number of potential ligand–receptor pairs between cell groups. B) Bubble plots showing ligand–receptor pairs of immunomodulation, growth factors, angiogenesis, and adhesion between NPCs and other cell groups. C) Predicted regulatory network centered on NPCs.

CD74/APP and CD44/HBEGF widely existed in the interaction between immune cells and NPCs. NPCs also produced VEGFB and VEGFA, which bind to VEGF receptors (FLT1, ADRB2, NRP1, and NRP2) in both NPCs and immune cells, especially macrophages. This implies that the interaction between macrophages and NPCs is potentially proangiogenic. Notably, a healthy IVD is largely avascular. However, angiogenesis occurs during degeneration.[ ]

It has also been suggested that NPCs produce FN1 and FGF for their receptors. Immunomodulation‐related cytokines, such as pro‐inflammatory cytokine, TNF, anti‐inflammatory cytokines, and IL‐6, are expressed by macrophages, neutrophils, and G‐MDSCs, and can bind to NPC receptors. Neutrophils and macrophages expressed growth factors, including TGFB, PDGF, and OSM, and showed pro‐proliferation effects on NPCs (Figure 6B,C). Overall, our results predicted a network between immune cells and degenerative NPCs, which may play a role in inflammation balance, cell proliferation, and angiogenesis in NP. However, the exact mechanisms should be confirmed by further experiments.

Discussion

Our identified subpopulations showed transcriptomic heterogeneity and potential functional differences at the single‐cell level. Effector NPCs were identified and indicated to be active in cellular anabolism and catabolism activities, as they showed a high metabolic rate, including sterol biosynthesis, glycosaminoglycan metabolism, and protein metabolism. During IVDD, the ECM and cells undergo profound metabolic processes.[ ] The QuSAGE analysis indicated effector NPCs are active for “positive regulation of ECM assembly,” which is opposite to fibroNPCs with “negative regulation of ECM assembly.” Thus, effector NPCs may play a role in maintaining ECM homeostasis. It was also suggested that with IVDD progressing, the proportion of effector NPCs decreased, while fibroNPCs proportion increased. Regulatory NPCs are characterized by inflammation and endogenous stimuli responses. The Toll‐like signaling pathway, a crucial pathway related to the innate immune system, is active in regulatory NPCs. Moreover, regulatory NPCs have been shown to be one of the main subsets in grade II NP degeneration, suggesting that they may serve as a cell population for immunity maintenance. Calió et al. found 15 distinct cell clusters in healthy bovine NP tissues. Their study demonstrated that NP is heterogeneous, they also identified notochordal‐like cell cluster and a progenitor stem cell cluster inside the bovine NP.[ ]

During cartilage development, chondrocytes undergo terminal differentiation when they become hypertrophic. It has been widely acknowledged that chondrocyte hypertrophy‐like changes play a role in osteoarthritis (OA).[ ] It has been shown that hypertrophic differentiation occurs in NP tissues during IVDD progression.[ ] Here, we annotated a cell population as “HT‐CLNPs” and identified HT‐CLNPs‐I and HT‐CLNPs‐II subpopulations possessing different ontologies. The HT‐CLNP‐I subpopulation was enriched in genes related to negative regulation of apoptosis, response to unfolded protein, and circadian rhythm, while the HT‐CLNP‐II subpopulation expressed genes involved in ECM organization and disassembly (Table 3). This indicates that HT‐CLNP‐I may serve as a regulatory cell cluster that may play a potential protective role in chondrocyte hypertrophy‐like events. In support, our data showed decreased HT‐CLNP‐II proportion in mild degenerative discs (grade II and III). Moreover, a distinct cell marker of HT‐CLNP‐II, CHRDL2, was found to negatively regulate chondrocyte development by competitively inhibiting BMPs.[ ] BMP/TGFβ is a vital pathway promoting chondrocyte hypertrophy. Here, functional analysis (Figure 3I) showed different cell death events occurring in HT‐CLNPs compared to other subpopulations. In HT‐CLNPs, ferroptosis, instead of apoptosis, was the main cell death type. HT‐CLNPs also showed decreased cell senescence and a low level of SASP. These results revealed potential cell dysfunction events happened in NPCs chondrocyte hypertrophy‐like changes.

It has been reported that NP originates from the notochord; however, at about age 4 in humans, notochordal cells disappear, replaced by smaller non‐vacuolated cells. More importantly, the origin of these cells is unknown. Previous studies indicated that Tie2 positive (Tie2+) and disialoganglioside 2 positive (GD2+) populations serve as progenitor NPCs.[ ] However, few Tie2+ and GD2+ NPCs can be found in patients aged over 50.[ ] Interestingly, the single‐cell data showed that there are few Tie2+ or GD2+ NPCs in degenerative NP tissues (Figure S5B, Supporting Information). IVDD is an age‐related condition; the average age is 59.13 ± 12.5 for males and 61.02 ± 10.6 for females.[ ] Thus, finding out progenitor cells in aged and degenerative NP tissues is vitally important for understanding tissue homeostasis and regeneration. Blanco et al. identified the presence of nucleus pulposus mesenchymal stem cells (NP‐MSCs) within the IVD. And their study showed the NP‐MSCs do not have adipogenic differentiation ability.[ ] Our scRNAseq showed that fibroNPCs demonstrated the highest score for progenitor properties and stemness. Increasing evidence shows that apoptotic cells can induce compensatory proliferation and promote regeneration in invertebrates and vertebrates.[ ] The relationship among apoptosis, proliferation, and tissue regeneration has been linked in some animal models, like drosophila, planarians, newts, and mice.[ ] Cell apoptosis is an important cell death mechanism in IVDD.[ ] Our study showed that apoptosis and inflammasome activity occurred mainly in fibroNPC populations. Moreover, the score for angiogenesis is highest in fibroNPCs and angiogenesis is essential for the growth and regeneration of tissues.[ ] These imply that fibroNPCs are a significant subpopulation for cell regeneration in NP tissues.

We selected CD90+ NPCs as candidate progenitors, which transcriptionally co‐express genes for fibroNPCs and MSCs. CD90+NPCs showed decreased apoptosis compared to CD90‐ NPCs and could differentiate into adipocytes, osteoblasts, and other NPCs atlas (like adhesion NPCs). Furthermore, the CD90+NPCs isolated from mild degeneration individuals (grade II and III) demonstrated increased cell chondrogenesis ability. This study revealed that CD90+ is important for maintaining NP homeostasis and can be used in cell‐based scaffolding for NP repair and regeneration. Moreover, CD90+ was found most abundant in the NP/AF boundary regions in rat degenerative disc tissues. Exploring the origin of progenitor cells is of researchers’ interest. Previous descriptive studies found that progenitor and stem cell niche patterns in healthy IVDs are detected in certain AF areas. It was possible that MSCs from the bone marrow niche were involved in NP origination.[ ] Here, we hypothesized that a potential niche with CD90+NPCs enriched exists in the NP/AF boundary (Figure S5C, Supporting Information). CD90+NPCs may originate from the AF region, perichondrium region, or from bone marrow‐MSCs differentiation and migration. Detailed experiments should be conducted to trace the origin and differentiation of CD90+NPCs in the future.

The IVDs have been identified as immune privilege organs; the steady state of immune privilege is fundamental to organ homeostasis.[ ] Although the NP is often described in general terms in the literature as being avascular, previous studies have reported blood vessels can penetrate the NP, especially when degeneration. Binch et al. reported the “abnormal” localization of blood vessels within the NP associated with degeneration.[ ] Freemont et al. indicated the blood vessels occurs into normally aneural and avascular regions of the IVD, notably the NP in painful IVDD patients.[ ] Considering the angiogenesis might link to disturbed cell biology, we explored the correlation between angiogenesis and NPCs heterogeneity. We found that the angiogenesis scores are highest in firboNPCs, but lowest in homestatic NPCs (Figure 2E). We also founded the homeostatic NPCs showed negative correlations with degeneration severity and fibroNPCs showed positive correlations with severity (Figure 1D). The angiogenesis in NP is also associated with in‐growth of nerves, which contribute to chronic back pain symptoms.[ ] David et al. showed the process of angiogenesis at the IVD affects the patient's quality of life both pre and postoperatively.[ ] The angiogenetic NPCs populations might contribute to angiogenesis by express angiogenic factors like vascular endothelial growth factor (VEGF). The neovascularization can bring cytokines and immune cells. The endothelial cells of these blood vessels produce nerve growth factor (NGF) which previously has been shown to be responsible for the in‐growth of nociceptive nerve fibers into the normally aneural IVD.[ ] The angiogenetic NPCs populations therefore might accelerate degeneration and lead to chronic back pain. However, more experiments should be conducted to elucidate their significance in the future.

This single‐cell study showed there are multiple immune cell lines inside the NP, which may play a role in IVDD progression. Based on correlation analyses of previously reported single‐cell data, we showed that G‐GMPs aligned with bone marrow GMPs, proNeu, and preNeu. G‐MDSCs aligned with CXCR2low immNeu and CXCR2high mNeu, and neutrophils aligned with polymorphonuclear leukocytes (PMNs).[ ] G‐MDSC represents a group of immature broadly defined neutrophils with immunosuppressive functions characterized as CD14−/CD11b+/CD15+/CD66b+ in human.[ ] MDSCs have been widely studied in cancer, acting as a suppressor of antitumor immune responses.[ ] Accumulation of G‐MDSCs has also been reported in cancers and other inflammatory diseases, including inflammatory bowel diseases, rheumatoid arthritis, autoimmune arthritis, and autoimmune hepatitis.[ ] Nevertheless, the involvement of MDSCs in IVDD has not been revealed. MDSCs have a close correlation with neutrophils. Our findings indicated that in NP tissues, G‐MDSCs might emerge from GMPs through a differentiation trajectory.

Although we identified CD24, which can be used in combination with CD11b and OLR1 to detect the presence of G‐MDSCs in degenerative NP tissues, this finding does not deny the CD24 as a marker for NPCs. The expression of CD24 is decreased with IVDD severity.[ ] Previously, CD24+ NPCs were identified as progenitors in NP.[ ] In humans, the percentage of CD24+NPCs decreased significantly with aging or degeneration, which dropped to less than 10% in elderly adults or severely degenerated discs.[ ] In our study, scRNAseq showed that there are few CD24+ NPCs in degenerative discs. That is because there is only one participant below 50 years old in this study. Here, we identified CD24 as a cell marker for detecting NP‐derived G‐MDSCs when combined with CD11b and OLR1.[ , ] The results of FACS showed that most CD11b+OLR1+ were CD24 positive. Thus, we can conclude that CD24 is mainly expressed in G‐MDSCs, but not in NPCs, especially in elderly populations. The immunofluorescence for CD24 in discs of rat models displayed positivity for NPCs, and the percentage of positive cells decreased with IVDD.[ ]This might be a result of age and species difference as the rats we used were 2–3 months old. It is currently known that CD24 serves as a costimulatory factor of T cells, regulating their homeostasis and proliferation, and is involved in B cell activation and differentiation.[ ] Our findings showed CD24 can regulate suppression of T cell responses, suggesting CD24 might allow G‐MDSCs to directly regulate immune checkpoints in NP tissues. Our single‐cell data and animal experiments showed that G‐MDSCs substantially decreased in mildly degenerated NP tissues. Our co‐culture experiment also suggested G‐MDSCs might be a potential treatment option for IVDD.

Although our study provides some previously unrevealed insights into NP tissue biology, there are potential limitations to our study here. First, we need to acknowledge that part of the functional analysis was based on scRNAseq prediction. For instance, the functional difference was among six identified NPC subsets. Other alternative methods should be used to estimate the functional biology behind transcriptional changes in future research. Second, the dissociation bias when processing tissues for scRNAseq may lead to spurious changes of cell distribution, which can limit our ability to provide an exact numeric description of cell population changes. Several strategies were undertaken to minimize the bias. First, all samples were handled following the same protocol. Second, sample acquisition was carefully made during the surgery to avoid contamination of AF tissues or other tissues. Third, we used multiple methods to show the presence of identified novel cell types, like G‐MDSCs. In summary, we have comprehensively decoded the multicellular ecosystem of NP during degeneration. These findings may be leveraged to improve diagnostics and develop preventative strategies for degenerative spinal diseases.

Experimental Section

NP Tissues Specimens

This study protocol was approved by the Ethics Committee of the First Affiliated hospital of Soochow University. All participants signed a written informed consent form. The enrolled subjects were patients who required discectomy and/or interbody fusion for burst fracture (n = 2) lumbar disc herniation or spondylolysis (n = 6). The patients’ characteristics were listed in Table 1. In this study, the grading of IVDD based on Pfirrmann grading system. It should be noted that it is difficult to obtain Pfirrmann grade I disc tissues, especially in elderly participants. Thus, only grade II to V discs were analyzed in this study (Figure S1A, Supporting Information).

A standard surgical protocol assured correct sample acquisition and thus accurate selection of AF and NP tissue by experienced surgeons. For burst fracture cases, the NP tissue from the disc attached to the intact bottom endplate was harvested. For severe degenerative discs, the tissue from the central region of the NP was harvested to assure no AF tissues were harvested. If local bleeding happened during the tissue harvested, the tissues were discarded to ensure no contaminated blood in collected tissues.

Tissues immediately following surgical procurement were processed, and single‐cell suspensions were generated within ≈45 min with an experimental protocol optimized to reduce artifactual transcriptional changes introduced by disaggregation, temperature, or time. Human NP tissue samples were transported in MACS Tissue Storage Solution (Miltenyi Biotec) immediately after surgical resection. After tissues were harvested, the specimens were washed with PBS for three times and further checked under a dissecting microscope to guarantee there was no contamination of blood, AF tissues, or other tissues. For few cases which were difficult to distinguish between the NP and inner AF, histological staining was used to confirm the absence of AF tissues or other tissues (Figure S1B–D, Supporting Information).

Single‐Cell Dissociation of Human Nucleus Pulposus

The tissues were surgically removed and kept in MACS Tissue Storage Solution (Miltenyi Biotec) until processing. The tissue samples were processed as described below. Briefly, samples were first washed with phosphate‐buffered saline (PBS), minced into small pieces (≈1 mm3) on ice and enzymatically digested with 500 U mL−1 collagenase I, 150 U mL−1 collagenase II, 50 U mL−1 collagenase IV, 0.1 mg mL−1 hyaluronidase, 30 U mL−1 DNaseI, and 5% Fetal Bovine Serum Oringin South America (Yeasen) for 95 min at 37 °C, with agitation. After digestion, samples were sieved through a 70 µm cell strainer, and centrifuged at 300×g for 5 min. After washing with PBS containing 0.04% BSA, the cell pellets were re‐suspended in PBS containing 0.04% BSA and re‐filtered through a 35 µm cell strainer. Dissociated single cells were then stained with AO/PI for viability assessment using Countstar Fluorescence Cell Analyzer. The single‐cell suspension was further enriched with a MACS dead cell removal kit (Miltenyi Biotec).

Single‐Cell RNA Sequencing

BD Rhapsody system was used to capture the transcriptomic information of the (8 sample‐derived) single cells. Single‐cell capture was achieved by random distribution of a single‐cell suspension across >200 000 microwells through a limited dilution approach. Beads with oligonucleotide barcodes were added to saturation so that a bead was paired with a cell in a microwell. Cell lysis buffer was added so that poly‐adenylated RNA molecules hybridized to the beads. Beads were collected into a single tube for reverse transcription. Upon cDNA synthesis, each cDNA molecule was tagged on the 5′ end (i.e., the 3′ end of an mRNA transcript) with a unique molecular identifier (UMI) and cell label indicating its cell of origin. Whole transcriptome libraries were prepared using the BD Rhapsody single‐cell whole transcriptome amplification workflow. In brief, second strand cDNA was synthesized, followed by ligation of the WTA adaptor for universal amplification. Eighteen cycles of PCR were used to amplify the adaptor‐ligated cDNA products. Sequencing libraries were prepared using random priming PCR of the whole‐transcriptome amplification products to enrich the 3′ end of the transcripts linked with the cell label and UMI. Sequencing libraries were quantified using a High Sensitivity DNA chip (Agilent) on a Bioanalyzer 2200 and the Qubit High Sensitivity DNA assay (Thermo Fisher Scientific). The library for each sample was sequenced by illumina sequencer (Illumina, San Diego, CA) on a 150 bp paired‐end run.

Single‐Cell RNA Statistical Analysis

fastp[ ] with default parameter was applied filtering the adaptor sequence, and the low quality reads were removed to achieve the clean data. UMI tools were applied for Single Cell Transcriptome Analysis to identify the cell barcode whitelist. The UMI‐based clean data were mapped to human genome (Ensemble version 91) using STAR mapping with customized parameter from UMI tools standard pipeline to obtain the UMIs counts of each sample. Cells contained over 200 expressed genes and mitochondria UMI rate below 20% passed the cell quality filtering and mitochondria genes were removed in the expression table. Seurat packageSeurat package (version: 2.3.4, https://satijalab.org/seurat/) was used for cell normalization and regression based on the expression table according to the UMI counts of each sample and percent of mitochondria rate to obtain the scaled data. To correct the batch effect, fastMNN function from scater package (https://github.com/Alanocallaghan/scater/tree/master/R) was applied with k value equals 5 and UMAP as well as tSNE dimension reduction construction was calculated following the batch effect correction result. PCA was constructed based on the scaled data with top 2000 high variable genes and top 10 principals were used for tSNE construction and UMAP construction. Using graph‐based cluster method (resolution = 0.8), the unsupervised cell cluster result based on the PCA top 10 principals acquired, and the marker genes were calculated by FindAllMarkers function with Wilcox rank sum test algorithm under following criteria: lnFC > 0.25; p‐value < 0.05; and min.pct > 0.1. To identify the cell type detail, the clusters of same cell type were selected for re‐tSNE analysis, graph‐based clustering, and marker analysis. To identify differentially expressed genes among samples, the function FindMarkers with Wilcox rank sum test algorithm was used under following criteria: lnFC > 0.25; p value < 0.05; and min.pct > 0.1.

Pseudo‐Time Analysis

The single‐cell trajectories analysis was applied using Monocle2 (http://cole‐trapnell‐lab.github.io/monocle‐release) using DDR‐Tree and default parameter. Before Monocle analysis, marker genes of the Seurat clustering result and raw expression counts of the cell passed filtering were selected. Based on the pseudo‐time analysis, branch expression analysis modeling (BEAM Analysis) was applied for branch fate determined gene analysis.

More single‐cell trajectories analyses were applied by using PAGA in scanpy package (https://scanpy.readthedocs.io/en/latest/index.html, version 1.6.0) and Slingshot (https://bioconductor.org/packages/release/bioc/html/slingshot.html, version 1.4.0). Before analysis, marker genes of the Seurat clustering result and raw expression counts of the cell passed filtering were selected.

Cell Communication Analysis and SCENIC Analysis

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. Membrane, secreted and peripheral proteins of the cluster of different time point were annotated. 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 assess transcription factor regulation strength, the single‐cell regulatory network inference and clustering (pySCENIC, v0.9.5) workflow was applied, using the 20 000 motifs database for RcisTarget and GRNboost.

QuSAGE Analysis (Gene Enrichment Analysis)

To characterize the relative activation of a given gene set such as pathway activation, “CellDeath_Ferrdb_soring” and “CellDeath_Inflammasome_review” as described before, QuSAGE (2.16.1)[ ] analysis.

Cell Fate Analysis, Including CytoTRACE and Velocity

The scVelo package with default parameter was applied for studying the cellular differentiation status based on the bam mapping file from UMI tools STAR mapping steps to solve the transcriptional dynamics of splicing kinetics using a likelihood‐based dynamical model. CytoTRACE Analysis was applied for cell stem stage analysis with default parameter.[ ]

Pathway Analysis

Pathway analysis was used to find out the significant pathway of the marker genes and differentially expressed genes according to KEGG database. Fisher's exact test was used to select the significant pathway, and the threshold of significance was defined by P‐value and FDR.

Rat Model of IVDD

Animal experiments were conducted following the guidelines of the Ethical Committee of The First Affiliated Hospital of Soochow University. A total of 15 male Sprague‐Dawley rats, aged three months, were used for the experiments in vivo. Ten rats underwent the surgeries, and the remaining five rats underwent no surgical intervention as negative controls. On the day of surgery, the animals were anesthetized in the induction chamber with oxygen flow at 1 L min−1 and up to 4% isoflurane. The isoflurane was reduced to between 2% and 2.5% once the animal was asleep. Once anesthetized, the animal was placed on a heating pad in supine position. Buprenorphine (0.1 mg kg−1) in saline were administered by subcutaneous injection for post operation pain relief. The intervertebral space was located by digital palpation. Next, a needle was affixed using a clamp so that 5 mm tip sticks out. After cleaned injection site with ethanol, the needle head was then inserted into the intervertebral space and held in place for 20 s to ensure puncture. After surgeries, the animals were returned to a warm and clean cage and monitored until wake.

Quantitative Real‐Time PCR

The samples used for RT‐qPCR were independent sample pools from the sampled used for sequencing. The total RNA was harvested from the NPCs, and it was reverse transcribed into cDNA with SuperScript II reverse transcriptase (Invitrogen, Cat No. 18064014). The sequence of primers for RT‐PCR are as below: MSMO1 (Forward 5′‐3′: TGCTTTGGTTGTGCAGTCATT; Reverse 5′‐3′: GGATGTGCATATTCAGCTTCCA) HMGS1 (Forward primer: GAGCCCATACTCATCAAGTACCG; Reverse primer: CCTCGGGAGAGATGCACAC); DKK (Forward primer: ACGAGTGCATCATCGACGAG; Reverse primer: GCAGTCCCTCTGGTTGTCAC); FN1 (Forward primer: AGGAAGCCGAGGTTTTAACTG; Reverse primer: AGGACGCTCATAAGTGTCACC); CRTAC1 (Forward primer: TGTCCAGGATGTTACCGTTCC; Reverse primer: AGCTGGGTGGGATTACTGTCA); COL1A1 (Forward primer: ATCAACCGGAGGAATTTCCGT; Reverse primer: CACCAGGACGACCAGGTTTTC); COL3A1 (Forward primer: GGAGCTGGCTACTTCTCGC; Reverse primer: GGGAACATCCTCCTTCAACAG); MMP2 (Forward primer: CCCACTGCGGTTTTCTCGAAT; Reverse primer: CAAAGGGGTATCCATCGCCAT); RPS29 (Forward primer: CGCTCTTGTCGTGTCTGTTCA; Reverse primer: CCTTCGCGTACTGACGGAAA); RPS21 (Forward primer: AGCAATCGCATCATCGGTG; Reverse primer: CCCCGCAGATAGCATAAGTTTTA); CHI3L1 (Forward primer: AAGCAACGATCACATCGACAC; Reverse primer: TCAGGTTGGGGTTCCTGTTCT); NF‐kB (Forward primer: AACAGAGAGGATTTCGTTTCCG; Reverse primer: TTTGACCTGAGGGTAAGACTTCT); CXCL2 (Forward primer: CTTGTCTCAACCCCGCATC; Reverse primer: CAGGAACAGCCACCAATAAGC); FMOD (Forward primer: ATTGGTGGTTCCACTACCTCC; Reverse primer: GGTAAGGCTCGTAGGTCTCATA).

Scoring of Biological Processes

Individual cells were scored for their expression of gene signatures representing certain biological functions. The functional signatures were derived from the Gene Ontology database. The immune response and steroid biosynthetic process were measured by GO:0006955 and GO:0006694, respectively. The innervation score was measured by the calculating the average expression of genes in the Gene Ontology term “innervation” (GO:0060384). The angiogenesis was measured by term “positive regulation of angiogenesis” (GO:0045766). The indicated scores were calculated by scaling the normalized expression of a gene across all cells. Gene weights were set to either 1 or −1 to reflect positive or negative relationships.

Immunohistochemical (IHC) assays

The samples used for IHC were independent sample pools from the sampled used for sequencing. NP tissues were fixed for 48 h in 4% buffered paraformaldehyde. The sections were pre‐treated for 10 min with trypsin (0.05%) and then treated with 3% (vol/vol) H2O2 for 15 min. Then, the sections were blocked at room temperature for 1 h with 10% goat serum. After washing with PBS, sections were incubated with anti‐CRTAC1 (1:50 dilution, ab25469, Abcam), anti‐CHI3L1 (1:250 dilution, ab255297, Abcam), anti‐MSMO1 (1:100 dilution, ab203587, Abcam), anti‐RPS14 (1:50 dilution, Abcam), anti‐FRZB (1:100 dilution, ab205284, Abcam), and anti‐MMP2 (1:200 dilution, ab86607, Abcam) antibodies overnight at 4 °C. The sections were then washed with PBS and incubated with a biotinylated secondary antibody for 15 min from a Histostain Plus kit (Invitrogen, CA, USA). The sections were then washed and incubated with 3, 3′‐diaminobenzidine for 2 min. Finally, using light microscopy to observe the section.

Cell Sorting

Tissue samples were harvested from patients with IVDD and mechanically dissociated to generate single‐cell suspensions as described above. Cells were blocked with FcR Blocking Reagent, human (Miltenyi, 130‐059‐901) on ice for at least 10 min. Cells were then centrifuged at 300g for 5 min at 4 °C and washed once with BSA running buffer (0.5% BSA). Cells were incubated for 30 min at 4 °C with pre‐conjugated fluorescent labeled antibodies with the following combinations: CD45 (BioLegend, 368532 (PE/Cyanine7)), CD11b (BioLegend, 301330 (FITC)), OLR1 (BioLegend, 358604 (PE)), and CD24 (BioLegend, 311118 (APC)). Cells sorted by Beckman Moflo XDP and desired populations were isolated for different experiments. For CD90+ NPCs isolation, the primary cells were used within one passage. Cells were incubated with CD90 microbeads (Miltenyi Biotec, 130‐096‐253; 1 µL per 107 cells) to enrich CD90+ NPCs.

Multilineage Differentiation Assays In Vitro

For osteogenic differentiation, the CD90+ NP cells were cultured in a humidified incubator (37 °C, 5% CO2) with renewal of the culture medium every 3 d. The cells were incubated in a differentiation medium for 2–4 weeks once the cells had reached 100% confluence, during which time the medium was changed every 2–3 d. The differentiation medium was as follows: DMEM‐LG supplemented with 10% FBS, 1 × 10−6 m dexamethasone (Sigma), 50 µg mL−1 ascorbic acid (Sigma), 10 × 10−3 m sodium β‐glycerophosphate (Sigma), and 1% penicillin/streptomycin (Sigma). The cells were fixed with ice‐cold 70% ethanol and stained with Alizarin Red S (Amresco, Solon, OH), as well as the von Kossa stain, to detect mineralization (calcium deposits). The alkaline phosphatase (ALP) activity was also tested. For adipogenic differentiation, the CD90+ cells were first grown to 100% confluence and then incubated for 3 d in an induction medium consisting of DMEM‐LG supplemented with 10% FBS, 100 × 10−6 m indomethacin (Sigma), 0.1 × 10−6 m dexamethasone, 0.5 × 10−3 m IBMX, 10 µg mL−1 human insulin, and 1% penicillin/streptomycin. The cells were incubated in the induction and maintenance media for >2 weeks and then fixed with 4% paraformaldehyde for 30 min at room temperature and stained with Oil Red O, as well as Sudan Black B, to detect fat deposition. For chondrogenic differentiation, cells were maintained as pellet cultures (2.5–5 × 105 cells/pellet) in DMEM high‐glucose medium supplemented with insulin transferin selenium (with albumin), sodium pyruvate (100 µg mL−1), l‐proline (40 µg mL−1), ascorbate 2‐phosphate (50 µg mL−1), and TGF‐β3 (10 ng mL−1).

SDS‐PAGE Immunoblotting

SDS‐PAGE Immunoblotting was performed to analyze protein levels. The cells were lysed using RIPA lysis buffer, and protein concentrations were determined using the BCA assay. Nuclear and cytoplasmic proteins were lysed using the Nuclear/Cytosolic Fractionation assay kit from BioVision (Mountain View, CA) following the manufacturer's protocol. The membrane was blocked with 5% non‐fat milk and then incubated with the following primary antibodies at 4 °C overnight: anti‐BAX (1:1000), anti‐Bcl‐2 (1:1000), anti‐NLRP3 (1:1000), anti‐aggrecan (1:1000), and anti‐MMP‐13 (1:4000). After this incubation, the membranes were incubated for 2 h on a shaker at 37 °C with horseradish peroxidase (HRP)‐conjugated secondary antibodies (Boster, Wuhan, China) used at a dilution of 1: 2000 and then washed. Finally, the protein bands were visualized and detected using the enhanced chemiluminescence (ECL) system, and immunoreactive bands were quantified using the ImageQuant LAS 400 software (GE Healthcare Life Sciences) and calculated by normalization to the reference bands of GAPDH or lamin B.

Histopathology and Immunofluorescence

For histopathological analysis, tissues were fixed in 10% formalin for 96 h, decalcified in 10% ethylenediaminetetraacetic acid (EDTA) for at least one month, embedded in paraffin, and 4–5 µm sections were stained with Safranin O/Fast Green Stain. For immunofluorescence, the caudal disc sections were blocked in 5% normal serum (Thermo Fisher Scientific, 10000 C) in PBS‐T (0.4% Triton X‐100 in PBS) and incubated with the primary antibody. The primary antibodies used were CD11b (1:500, 120772, Absin), OLR1(1:500, Absin, 123947), CD24 (1:200, Proteintech, 10600‐1‐AP). After washed with PBS/1% BSA, the sections were incubated with Alexa‐conjugated antibodies and washed with PBS/1% BSA followed by PBS. The DAPI was used for counterstaining.

T‐Cell Suppression Assay

T‐lymphocytes were isolated from healthy donor's PBMCs via Pan T Cell Isolation Kit (Miltenyi Biotec, 130‐096‐535) following the manufacturer's instructions. Isolated T‐cells were then cultured in 96‐well round‐bottom plates in complete culture medium containing soluble, or plate‐bound, anti‐CD3 (1 µg mL−1) and soluble anti‐CD28 (5 µg mL−1). Sorted CD45+CD11b+OLR1+CD24+ and CD45+CD11b+OLR1+CD24− cells were added to T‐cells in 1:1 ratio. After 4 d of culture, T‐cell proliferation was assessed by MTT assay. T‐cells without sorted cells co‐cultured were used as positive control. Supernatants of the co‐culture were collected, the level of interleukin 2 (IL2) were measured by ELISA (Abcam) according to the manufacturer's protocol.

ROS Production Assay

After isolating NP‐derived MDSCs by fluorescence‐activated cell sorting, the cellular ROS level was determined with the ROS assay kit. The cell suspension was used at a concentration of 1 × 106 cells mL−1 and was diluted with DCFH‐DA solution (10 × 10−6 m) and incubated at 37 °C for 20 min with mixing by inversion every 5 min. After washing the cells with serum‐free medium, the sample was analyzed by flow cytometry (FAC Scan, Becton Dickenson, USA). Rosup was used as positive control.

Conflict of Interest

The authors declare no conflict of interest.

Author Contributions

J.T. and W.L. and S.D.Y. contributed equally to this work. J.T. conceptualized the project and planned, executed, and prepared the work for publication. J.Z. provided supports, supervised and initiated this study. J.T. performed most experiments, analyzed the scRNAseq data, and produced all figures and tables. W.T.L. and S.D.Y. performed part of the experiments and helped in data analysis. J.T. and W.T.L. wrote the manuscript with inputs from all authors. P.Y.Y. provided comments and feedback to the scRNA‐seq and manuscript. Q.Y. and C.H.W. collected samples and clinical data. K.T.L. and X.P.B. edited the manuscript. W.Y.D., J.C.W., A.D.D., C.Y., and H.L.Y. revised the manuscript, provided comments, and coordinated the collaboration.

Supporting Information

Click here for additional data file.