The microarray data and corresponding clinical information, including 30 human disc tissues, were downloaded and filtered from the GEO dataset (GSE124272, GSE147383, and GSE153761). The single-cell transcriptome data from intervertebral discs were obtained from GEO (GSE154884). Multiple datasets of gene expression matrices were merged and the inter-batch differences were removed for further processing.

The edgeR package was used to screen for DEGs in stage IV lumbar disc herniation (LDH) tissues and controls, including spondylolisthesis and LDH stages I to III. The selection criteria were |log2 FC| > 1.5, and false discovery rate (FDR) < 0.05. The volcano plot and heatmap of DEGs were generated by the pheatmap R package.

The CIBERSORT deconvolution algorithm was applied to quantify the degree of infiltration of 22 types of immune cells through the transcriptome data. Investigated immune cells included plasma cells, resting memory CD4+ T cells, CD8+ T cells, naive CD4+ T cells, T follicular helper cells, regulatory T cells (Tregs), activated memory CD4+ T cells, gamma delta T cells, naive B cells, memory B cells, monocytes, M0 macrophages, M1 macrophages, M2 macrophages, resting natural killer (NK) cells, activated NK cells, activated mast cells, eosinophils, neutrophils, resting dendritic cells, activated dendritic cells, and resting mast cells. The differences between the two groups were compared using the wilcox test and the results were visualized by applying the vioplot package. Finally, the correlation between infiltration rate of each types of immune cells was determined by the corrplot package.

Immune cell infiltration-related genes were identified using the WGCNA R package. First, the gene expression matrices outliers were filtered by hierarchical cluster analysis. Then, the correlation coefficient of genes was constructed and transformed into a weighted adjacency matrix. Next, these genes were allocated into minimum-sized modules and a cluster dendrogram was drawn, and then merged with a height cutoff (cutoff < 0.3). The correlation between gene expression and sample trait (immune cell infiltration score) was determined by the criterion of gene significance (GS) > 0.5 and module membership (MM) > 0.8. The relevant genes in the module were then tested for correlation with all other genes, the screening conditions were set at cor > 0.7 and p < 0.01, the regulatory network was visualized using cytoscape. Finally, matascape was used to conducted functional analysis to the selected genes.

Single-cell RNA sequencing data sets from IVD tissue (GSE154884) were obtained. The Seurat pipeline was used for data reprocessing and to classify the cell groups, SingleR identified the cell type, and Monocle was used to analyze the cell differentiation trajectory.

The clusterProfiler package (9) was used for performing gene ontology (GO) enrichment analyses on cell cluster marker genes, where p < 0.05 and FDR < 0.25 were considered significantly enriched.

Nine patients who underwent posterior open discectomy for radiating pain due to single-level LDH (classified as Pfirrmann grade IV by MRI) were included. Six patients were female and three were male (10). The mean duration between onset of symptoms and operation was 13.5 weeks. Patients with degenerative spondylolisthesis or a history of diabetes mellitus or renal disease were excluded. As a control, 14 patients, including 8 females and 6 males, who underwent posterior open discectomy for neurological symptoms due to single-level lumbar disc herniation (classified as Pfirrmann grades I to III by MRI) were selected for this study. Patient immune responses were evaluated as either “high” or “low” and the results are shown in Table 1 , to be specific, the immune status of each individual is evaluated by whether the percentage of total lymphocytes higher than the median level of 32 or the percentage of CD4-CD8- lymphocytes is lower than median level of 6. Tissues were excised and transferred to liquid nitrogen for RNA and protein extraction. Written informed consent was obtained from all patients. Clinical data were acquired from hospital records and pathology reports. The study protocol (approval number: 2020–171) was approved by The Ethics Committee of the Affiliated Hospital of Chongqing Medical University.

IVD specimens obtained from patients were first homogenized in phosphate-buffered saline (Tissue Tearor kits, Racine, model 985–370). The supernatant was obtained by centrifugation at 15,000 rpm at 4°C for 30 min. The Bradford protein assay method (#500–0006; Biorad, Hercules, CA, USA) was used for protein quantification and the results were used for measurement of TNF-α and TGF-β. IVD specimens were also used for RNA extraction using the UNIQ-10 Column Total RNA Purification Kit (Sangon Biotech, China). The quality and concentration of RNA were evaluated using a SMA4000 microspectrophotometer (Merinton Instrument, Inc. MI, USA).

IVD extracts were quantitatively analyzed for TNF-α and TGF-β using an enzyme-linked immunosorbent assay (ELISA) kit (Beijing Jingmei Biological Engineering Co, Ltd, China) and measurement at 450 nm. TNF-α and TGF-β concentrations were extrapolated from the standard curve.

The RR047A cDNA synthesis kit (TaKaRa, China) was used to perform the reverse-transcription of the extracted RNA and the 2X SG Fast qPCR Master Mix (High Rox, B639273, BBI) was used for quantitative PCR of hub genes on an ABI PRISM 3700 instrument (Foster, CA, USA). GAPDH was used as an internal control and primers are as follows:

A flow diagram showing the research design is shown in Figure 1 . In brief, DEGs in the LDH were screened from microarray data of samples in the GEO database. CIBERSORT was then applied to DEGs for identifying LDH-related immune cells. Next, WGCNA and correlation test were used to find hub genes associated with the identified immune cells. The expression pattern of these hub genes was then examined at the single-cell level and several LDH-related NP cells were recognized. Finally, qRT-PCR and ELISAs were performed to verify the relationship between hub gene expression levels and clinical characteristics in LDH patients.

Batch correction, normalization, and difference analysis of RNA-seq data from GSE124272, GSE147383, and GSE153761 were performed to screen for DEGs in IVD samples. A total of 410 DEGs, including 195 downregulated and 215 upregulated genes were identified. The results were visualized using a volcano plot ( Figure 2A ), which identifies important genes.

In order to further reveal the immune microenvironment in degenerated IVDs, the CIBERSORT method was used to analyze specific immune cell types that infiltrated into IVD tissue. Among the 22 types of immune cells investigated, the results showed that Tregs and macrophage levels were significantly higher in degenerated IVD (p<0.05) ( Figure 2B ), As is shown in Figure 2C , by further analyzing the CIBERSORT scores, there was a positive strong correlation between Mast cells, neutrophils and M1 macrophages, On the other hand, there was a negative correlation between the infiltration of plasma cells, regulatory T cells and M1 macrophages.

WGCNA was applied to identify differentially expressed immune cell-related genes from 410 DEGs, as shown in Figures 3A–D . These DEGs were then divided into modules and merged with different colors ( Figures 3A, B ). Five merged modules were analyzed and three gene modules were significantly correlated with immune cells ( Figure 3C ). Among them, we selected the green module, which was the most significant module that was related to macrophages M0, with a positive correlation of 0.58 and p<0.001 ( Figure 3D ). 15 genes were screened according to the criteria, of which 12 were selected by correlation test. The expression pattern of these key genes was then analyzed accordingly ( Figures 3E, F ). Cytoscape constructed the interaction network between these 12 genes and their target genes ( Figure 3G ). Next, to acquire more information about these genes we performed functional analysis using Matascape, as shown in Figure 3H , genes were enriched in critical biological processes, such as interleukin signaling and the regulation of cytokine production.

To determine the single-cell level transcriptomic landscape of IVD compartments, we employed scRNA-sequencing data from healthy, non-degenerated rat IVDs. We first conducted quality control of the gene expression matrix ( Figure 4A ). Then, normalization of scRNA-seq data was performed and 20 principal components (p<0.05) were screened for subsequent analysis ( Figure 4B ). Reduced dimension process analysis was achieved by using Discriminative Dimensionality Reduction Tree ( Figure 4E ). Unsupervised analysis was then conducted for cell clustering using the t-distributed stochastic neighbor embedding (t-SNE) method ( Figures 4C, D ). The result showed high cell heterogeneity, in which IVD cells were segregated into four major distinct clusters, including chondrocytes (NP cells), fibroblasts (AF cells), adipocytes, and epithelial cells, which were determined using singleR and cell markers. Next, we tested the expression pattern of the hub genes in these cell clusters. As expected, most of the key genes were highly expressed in NP cells, with the key gene ID1 found to be highly expressed in the AF cells ( Figure 4F ).

To uncover the detailed heterogeneity of NP cells, we re-subclustered the scRNA-seq data by upregulating the resolution value. Four cell clusters, including clusters 1, 2, 4, and 5, were identified according to their gene expression pattern ( Figure 5A ). Trajectory analysis was then performed to illustrate the degree of cell differentiation. All cells were projected onto one root and two branches. Interestingly, cells in cluster 5 were mainly located in the root, whereas cells in clusters 1, 4 were mostly located in the left, and most cells in cluster 2 were located in the right ( Figure 5C ). Based on the hub gene markers, clusters 1, 2, 4, and 5 were defined as RBB7+Pdlim5-, RBB7+ Pdlim5+, RBB7- ID1+, and RBB7- Pdlim5+ cell groups ( Figure 5B ). Next, we defined the molecular features of these NP cell groups by tracking their gene markers ( Figure 5D ). Among them, TOMM20, TOMM22, and CXCL2 were recognized as specific marker genes in cluster 1 ( Figures 6A–C ). TOMM20 has been shown to induce tissue inflammation responses in adipose and muscle tissue (11). CXCL2 is reported to promote tumor cell migration through the induction of M2-like macrophage polarization. Therefore, we named cluster 1 “NP inflammatory response cells”. Genes, including TOMM7 and placenta expressed transcript 1 (PLET1), are specifically expressed in cluster 2. PLET1 encodes for a cell surface protein that is specifically expressed in trophoblast stem cells. Of interest, functions enriched in cluster 2 included regulation of extrinsic apoptotic pathway and responses to metal ions, suggesting that cluster 2 cells are involved in multidirectional maintenance and differentiation of chondroid cells ( Figures 6B–D ). Therefore, cluster 2 was named “NP repair cells”. BNIP3, SOD2, and other genes are highly expressed in cluster 4. Functional enrichment results from cluster 4 marker genes showed that these cells function to response to oxygen levels. Cluster 4 cells provide extracellular matrix structural constituents, which confer resistance to compression ( Figures 7A, C ). Thus, we named cluster 4 “extracellular matrix NP cells”. Cluster 5 cells have highly expressed genes, including CLU, MEF2A, MCAM57, and FAM162A, which are known markers of mesenchymal stem cells and NP progenitor cells ( Figures 7B–D ). The trajectory result showed that cluster 4 cells were mostly populated in the root of two branches. We speculate that this cell group of the IVD contains NP progenitor cells, which other researchers have suggested. Thus, we refer to this group as the “stem-like NP cells”.

The Pfirrmann grade for each study group is represented by a set of lumbar intervertebral disc MRIs, as shown in Figures 8A, B . qRT-PCR results showed that the expression levels of ID1, RAP2C, and PTPRK were all significantly higher in NP tissues of the grade IV disc degeneration and high immunity groups ( Figures 8C–E ). ELISA results also showed that the expression of these genes was positively correlated with TNF-α ( Figures 8G–I ) but not with TGF-β (p<0.05, r=0.492). We believe that the expression of hub genes may play a role in chronic disorders such as IVD degeneration ( Figure 8F ).