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

Abstract

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.

Keywords: intervertebral disc degeneration; low back pain; nucleus pulposus; single-cell RNA sequencing.

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.

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.

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.

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.

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.

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.