Identification and Experimental Verification of Potential Immune Cell-Associated Gene Biomarkers in Human Intervertebral Disc Degeneration

Introduction

Low back pain has a high incidence in patients with skeletal and muscular diseases, and is the main cause of dysfunction, worldwide.^1–3 Intervertebral disc degeneration (IDD) is a major cause of low back pain; however, its specific pathogenesis is unclear.⁴ IDD is a prevalent degenerative disease that is usually accompanied by lower back pain, lumbar spinal stenosis, lumbar disc herniation, and sciatica.⁵ The progression of intervertebral disc (IVD) degeneration involves many etiological factors including heredity, aging, and biomechanical overload. However, the progression of degenerative IVD has not yet been fully elucidated.^6,7 At present, common IDD treatments include symptomatic and surgical treatments; however, these treatments can only relieve the symptoms of IDD and cannot fundamentally improve the condition.^8–10 Therefore, it is necessary to identify potential diagnostic markers and molecular therapeutic targets for IDD.

Starkweather et al suggested that the neuroimmune system is triggered during disc herniation-caused nerve root pain, which triggers immune cells to discharge proinflammatory cytokines, which signal the brain via humoral and neural pathways, causing ache and functional changes in neural action.¹¹ Immunity abnormalities are a crucial contributor to the development of IDD.^4,12 Several studies have examined the characteristics of immune-related inflammatory response and molecular immunology in low back pain.¹³ Increasing evidence suggests that immune infiltration plays an important role in the occurrence and development of IDD. IDD shows different patterns of immune cell infiltration in IVD tissues, and different cytokines secreted by immune cells participate in the initiation and development of IDD by regulating the inflammatory response. Many studies have examined immune cells in IVD tissues. For example, the main function of macrophages in herniated discs is to promote the reabsorption of prominent tissues, participate in the process of blood vessel ingrowth, and release a variety of inflammatory cytokines.^14,15 Dynamic polarization of macrophage M1 and M2 cell subtypes was found in the progression of IDD and played a significant role in regulating cell metabolism, especially in Pro nucleus pulposus cells.¹⁶ In addition, T cells promote the autoimmune response in the disc through the release of interleukin (IL) and TNF-α, leading to the occurrence of pain.¹⁷ The presence of B cells¹⁸ and natural killer (NK) cells¹⁹ in human discs has been reported. Furthermore, some studies have involved bioinformatic analysis of immune cell-associated gene,²⁰ but there have been few reports on the expression and role of immune cell-related genes in IDD.

This study conducted weighted gene co-expression network analysis (WGCNA) to identify immune cell-related hub genes. Next, three machine learning methods were used to screen candidate signature genes from IDD immune-related hub genes. Furthermore, a network containing 5 diagnostic markers and 24 long non-coding RNAs (lncRNAs) was constructed. Through the above analyses, five diagnostic markers were identified, which provide a theoretical foundation and have reference value in the study of IDD from an immunological perspective.

Methods

Data Source

The Gene Expression Omnibus Series (GSE)150408 and GSE124272 datasets were from the Gene Expression Omnibus (GEO) database (https://www.ncbi.nlm.nih.gov/geo/). The GSE150408 dataset includes gene expression profiles of whole-blood samples of 17 IDD and 17 healthy control samples. Messenger RNA (mRNA) and lncRNA expression profiles were primarily used in this study. Whole blood samples from eight IDD and eight healthy control samples were obtained from the GSE124272 dataset.

Differential Expression (DE) Analysis

DE analysis of the GSE150408 and GSE124272 datasets was performed using the R package limma. |log₂ fold change (FC)| >0.5 and P value <0.05 were used as the significance thresholds. We obtained differentially expressed (DE)-mRNAs and DE-lncRNAs from the GSE150408 dataset and DE-mRNAs from the GSE124272 dataset. We performed an intersection analysis to identify robust DE-mRNAs. In brief, an intersection analysis was performed on the upregulated DE-mRNAs in both datasets to identify the commonly upregulated and downregulated mRNAs. These common mRNAs were defined as IDD-related DE mRNAs. Furthermore, intersection analysis was performed using the Venn Diagram online tool.

Functional Enrichment Analysis

Functional enrichment analyses were performed to explore potential biological functions of the target genes. In this study, gene ontology (GO) analysis was performed on IDD-related DE-mRNAs using the R package cluster Profiler. The GO system consists of three main categories: cellular component (CC), biological process (BP), and molecular function (MF). P value <0.05 indicated a significantly enriched term. However, no enriched terms for IDD-related DE-mRNAs were identified in the BP category.

Immune Infiltration Analyses

Immune infiltration was estimated using both the marker gene expression-based microenvironment cell population-counter (MCP-counter) and single-sample gene set enrichment analysis (ssGSEA) algorithms for each sample in each of the two datasets (GSE150408 and GSE124272) in GEO database. MCP-counter was implemented in the R package MCP-counter to quantify the abundance of different immune cells based on specific molecular markers.²¹ The Wilcoxon rank-sum test was used to examine the differences between immune cells from healthy controls and IDD samples. P <0.05.

Construction of the Co-Expression Network by a WGCNA

Co-expression networks and hub genes targeting immune cells (the DE and consistently expressed genes in both the GSE150408 and GSE124272 datasets) were constructed and analyzed using the R package WGCNA based on the mRNA expression profiles of the GSE150408 and GSE124272 datasets. In the GSE150408 and GSE124272 cohorts, sample aggregation was first assessed by cluster analysis. Next, a scale-free co-expression network was created using an adjacency matrix based on the soft threshold parameter, β. Hierarchical clustering analysis was then performed to identify the modules of genes with similar expression patterns. Similar modules in each cohort were then merged based on MEDissThres of 0.4. Finally, we calculated the Pearson correlations between the modules and specific features (immune cells) and identified the modules of interest based on a P value ≤0.05, and a |correlation coefficient (cor)| value >0.3. After selecting the modules of interest, we defined a |ModuleMembership (MM)| value >0.6 (the correlation between a phenotype and a gene) and a |GeneTraitSignificance (GS)| value >0.4 (the correlation between a gene expression profile and module eigengene) as the thresholds for screening the immune-related hub genes in the modules of interest.

To identify the IDD- and immune-related hub genes, we conducted an overlap analysis to identify common genes among the IDD-related DE-mRNAs, GSE150408-immune-related hub genes, and GSE124272-immune-related hub genes, which were then defined as IDD- and immune-related hub genes.

Protein-Protein Interaction Networks

Protein-protein interaction networks (confidence =0.15) were constructed to predict the protein-protein interactions of the IDD- and immune-related hub genes using the STRING database (https://cn.string-db.org/). Next, these data were uploaded into Cytoscape software to visualize the networks of the IDD- and immune-related hub genes.

Identification of Signature Genes

This study used machine learning, a novel algorithmic analysis tool, to identify signature genes that could be used to distinguish between healthy control samples and IDD samples. We extracted the expression profile data of IDD- and immune-related hub genes in the GSE150408 and GSE124272 datasets as feature subsets. The signature genes were selected by three machine-learning algorithms (Support vector machine-recursive feature elimination (SVM-RFE), last absolute convergence and selection operator (LASSO), and extreme gradient boosting decision trees (XGBoost)) in the GSE150408 data set and validated in the GSE124272 data set. LASSO is a regularized regression algorithm performed using the R package glmnet. SVM-RFE is a supervised learning technique that applies recursion to the rank features. We applied the SVM algorithm using the R package e1071. XGBoost is a decision-tree-based ensemble of machine learning algorithms implemented using the XGBoost package. The real feature genes were the common genes filtered by the three machine-learning algorithms. The ability of the identified signature gene/linear regression model (constructed based on the signature gene) to distinguish between IDD and healthy control samples was assessed by the area under the curve (AUC) of the receiver operating characteristic (ROC) curve. The ROC curves were analyzed using the pROC in R software.

Gene Set Enrichment Analysis (GSEA) of Signature Genes

GSEA was conducted using the GSE150408 dataset. We downloaded the C5 (c5.all.v7.5.1. entrez.gmt) and C2 (c2.cp.kegg.v7.5.1. entrez.gmt) sub-collections from the MSigDB as the reference gene sets to perform the GSEA analysis for the signature genes. Gene sets were considered significant if they had a normalized enrichment score| >1, adjusted P value <0.05, and q value <0.25.

Construction of lncRNA-mRNA Co-Expression Network

In this study, Pearson correlation analysis was performed to investigate the immune-related lncRNA targeting the signature genes based on the expression profiles of the signature genes and DE-lncRNAs extracted from the GSE150408 dataset. According to lncRNA-mRNA co-expression relationship pair theory, the relationship between lncRNAs and mRNAs should be positively correlated. Thus, we integrated and uploaded DE-lncRNA-featured gene relationship pairs with P values <0.05, and cor values >0.3 into the Cytoscape software for co-expression network visualization.

Real-Time Quantitative PCR (RT-qPCR)

Peripheral blood samples were collected from 20 healthy individuals and 20 IDD patients. Patients with intervertebral disc degeneration (IDD) were identified based on the presence of specific clinical symptoms and a Pfirrmann grade ranging from IV to V. In contrast, healthy controls exhibited no clinical symptoms and had a Pfirrmann grade of I to II.²² The recruitment period for participants was from December 1, 2023, to April 1, 2024, and the participants signed an informed consent form (No. E20231128083), and complied with the Declaration of Helsinki. The RNAiso Blood Kit (Takara Beijing, China) was used for total RNA extraction. A PrimeScript ™ II 1st Strand complementary DNA (cDNA) Synthesis Kit (Takara, Shiga, Japan) was used for cDNA synthesis. RT-qPCR was performed on a LightCycler 96 (Roche Life Sciences, Swiss Basel, Switzerland) using Real-Time PCR Mix (Cowin Biotech, Taizhou, Jiangsu, China). The 2^−ΔΔCt method was used to evaluate relative mRNA expression. GAPDH was used as the internal control. The independent experiments were repeated three times. Detailed primer sequences are shown in Supplementary Table 1.

Statistical Analysis

All statistical analyses were performed using R software. The expression levels of the signature genes were compared between patients with IDD and healthy controls using Wilcoxon’s test. P <0.05.

Results

Analysis of IDD-Related DE-mRNAs

The GSE150408 and GSE124272 datasets were obtained from the GEO database. After matching the mRNA expression profiles of the GEO database, differential expression (DE) analysis was performed. In total, 520 DE-mRNAs were identified in the GSE150408 dataset. Compared to the healthy control samples, 290 and 230 mRNAs were significantly upregulated and downregulated, respectively, in the IDD samples (Figure 1A). We also identified 2502 DE-mRNAs between IDD and healthy control samples in the GSE124272 dataset. Compared to the healthy control samples, 1296 and 1206 mRNAs were significantly upregulated and downregulated, respectively, in the IDD samples (Figure 1B). To identify robust DE-mRNAs, we performed overlap analysis of GSE150408-DE-mRNAs and GSE124272-DE-mRNAs. The results are shown in Figure 1C-D. 67 mRNAs were upregulated in both GEO datasets (IDD vs healthy control samples) and 33 mRNAs were downregulated in IDD in both GEO datasets (IDD vs healthy control samples), which were defined as IDD-related DE-mRNAs.

Figure 1 Identification of DE IDD-related genes. (A-B) Volcano plot of the DEGs in the IDD and control samples in the GSE150408 (A) and GSE124272 (B) datasets. (C-D) Venn diagram of the differentially upregulated (C) and differentially downregulated (D) key common genes. (D) Venn diagram of the key common genes. (E-F) GO terms in the CC (E) and MF (F) categories enriched by the differentially expressed IDD-related DE-mRNAs. Abbreviations: DE, Differential expression; IDD, intervertebral disc degeneration; DEGs, differentially expressed genes; GSE, Gene Expression Omnibus Series; GO, gene ontology; CC, cellular component; MF, molecular function.

Furthermore, we performed GO enrichment analysis to determine the potential biological functions of IDD-related DE-mRNAs. Five CC terms and 12 MF terms were found to be enriched. In the CC category, IDD-related DE-mRNAs were significantly associated with endocytic vesicles and the nicotinamide adenine dinucleotide phosphate (NADPH) oxidase complex (Figure 1E). Figure 1F shows the top 10 terms in the MF category, and these genes were found to be significantly associated with NADPH. Moreover, Immunoglobulin G (IgG) binding, immune receptor activity, and immunoglobulin binding were significantly enriched.

Immune Landscape in IDD and Healthy Controls

In the GSE150408 and GSE124272 datasets, the MCP-counter was first used to impute the contents of eight immune cells, one fibroblast, and one epithelial cell in the IDD and healthy control samples, and the differences between the two groups were assessed using the Wilcoxon rank-sum test. The results are shown in Figure 2A. In the GSE150408 dataset, natural killer cells were significantly downregulated in the IDD samples compared to the healthy control samples, while the content of neutrophils was significantly upregulated in the IDD samples (all P <0.05; Figure 2A). In the GSE124272 dataset, T cells, cytotoxic lymphocytes, and natural killer cells were reduced more in the IDD samples than in the healthy control samples, while neutrophils and endothelial cells were more abundant in the IDD samples (all P <0.05; Figure 2B). Furthermore, the ssGSEA approach was used for immune infiltration analysis of 28 immune cells in both the datasets. As shown in Figure 2C, neutrophil and MDSC expression in the GSE150408 cohort was higher in IDD samples than in controls (all P <0.05). For the GSE124272 dataset, the expression levels of activated CD4 T cells and natural killer CD56dim cells in the case samples were lower, while activated dendritic cells and mast cells were upregulated (all P <0.05; Figure 2D). The above results suggest that natural killer cells and neutrophils are robust differentially abundant immune cells in IDD and healthy control samples.

Figure 2 Immune cells infiltration was analyzed between the IDD and healthy control groups. (A-B) MCP-counter calculated in the GSE150408 (A) and GSE124272 (B) datasets. (C-D) Box plot of 28 immune cells infiltration between the IDD and healthy control groups through ssGSEA algorithm in GSE150408 (C) and GSE124272 (D). ns represents no significant difference, * P <0.05, ** P <0.01. Abbreviations: IDD, intervertebral disc degeneration; MCP, microenvironment cell population; GSE, Gene Expression Omnibus Series; ssGSEA, single sample gene set enrichment analysis.

Identification of Immune Cell-Related Hub Genes by WGCNA

We performed WGCNA to identify immune cell-related hub genes. Natural killer cells and neutrophils were used as phenotypic data. To construct a scale-free network, the soft threshold was set to two in the GSE150408 dataset (Figure 3A) and three in the GSE124272 dataset (Figure 3B). Subsequently, 20 and 22 modules were obtained from the GSE150408 (Figure 3C) and GSE124272 (Figure 3D) datasets, respectively, based on the dynamic shear tree algorithm and a MEDissThres of 0.4. In the GSE150408 dataset, the royal-blue module (cor =0.38, P=0.03) was significantly and positively correlated with natural killer cells, and the black (cor =−0.41, P=0.02), green (cor =−0.41, P=0.02), and dark grey (cor =−0.36, P=0.04) modules were significantly and negatively correlated with neutrophils and were considered to be interesting modules (Figure 3E). In the GSE124272 dataset, the black (correlated with neutrophils; cor =0.77, P=5e-04), blue (correlated with neutrophils; cor =−0.86, P=5e-02), coral1 (correlated with natural killer cells; cor =0.6, P=0.02), and light steel-blue1 (correlated with neutrophils; cor =0.86, P=3e-05) modules were considered interesting modules (Figure 3F). We then calculated the values between MM and GS in the modules of interest and obtained hub genes. In the GSE150408 dataset, 22, 100, 100, and 70 hub genes were identified in the royal-blue, black, green, and dark-gray modules, respectively (Figure 4A). In the GSE124272 dataset (Figure 4B), 9802 hub genes were identified in the four modules of interest (black: 2396; blue: 6963; coral1:363; light steel-blue1:80). The genes identified in the modules of interest were referred to as immune cell-related hub genes and used in subsequent analyses.

Figure 3 Identification of the module and construction of the weighted gene co-expression network involved in IDD. (A-B) Construction of a scale-free network using the GSE150408 (A) and GSE124272 (B) datasets. (C-D) Dendrogram of the genes clustered on the basis of the highly correlated eigengenes in the GSE150408 (C) and GSE124272 (D) datasets. (E-F) Heatmap of the correlations between the module eigengenes and different cells in the GSE150408 (E) and GSE124272 (F) datasets. Abbreviations: IDD, intervertebral disc degeneration; GSE, Gene Expression Omnibus Series.

Figure 4 The hub gene and reciprocal network. (A-B) Scatter plot of module eigengenes in the royal-blue, black, green, and dark-grey modules of the GSE150408 (A) and GSE124272 (B) datasets. (C-D) Overlap analysis between the IDD-related DE-mRNAs and GSE150408 (C) or GSE124272 (D) -immune cell-related hub genes. (E) The reciprocal network between the IDD- and immune-related hub genes. Abbreviations: IDD, intervertebral disc degeneration; GSE, Gene Expression Omnibus Series.

Identification and Evaluation of Signature Genes

To identify IDD- and immune-related hub genes, we performed an overlap analysis. As shown in Figure 4C, there were five common genes between IDD-related DE-mRNAs and the GSE150408-immune cell-related hub genes. Additionally, there were 68 common genes between IDD-related DE-mRNAs and the GSE124272-immune cell-related hub genes (Figure 4D). The integration of the above common genes resulted in 68 co-genes, which were defined as IDD- and immune-related hub genes. We also visualized the reciprocal network between these IDD- and immune-related hub genes, which contained 20 nodes and 36 edges, using STRING and Cytoscape software (Supplementary Figure S1). Furthermore, based on the GSE150408 dataset, candidate signature genes from 68 IDD and immune-related hub genes were screened. LASSO regression analysis identified the following 12 candidate signature genes based on a λ _min value of 0.042 (Figure 5A): ASAP1-IT1, IKZF2, KLHL14, KRTAP5-7, lnc-C10orf131-1, lnc-TBK1-1, lnc-TNFAIP6-2, lnc-ZNF438-1, LOC101927805, NRCAM, PDGFD, and RGL3. ROC curve analysis showed that the LASSO regression model, which comprised the above-mentioned 12 genes, was able to accurately distinguish between the IDD samples and healthy control samples in the GSE150408 dataset (AUC =1; Figure 5B). The SVM-RFE algorithm identified 59 candidate signature genes using 10-fold cross-validation, at a minimum error rate of 0.138 (Figure 5C). Fourteen candidate signature genes were identified from the 68 IDD and immune-related hub genes by XGBoost analysis (Supplementary Table 2). The ROC curve of the XGboost model had an AUC of 1 (Figure 5D), indicating that this 14-gene based model had a satisfactory ability to identify the sample types. We then performed an overlap analysis to identify the intersection between the LASSO-candidate signature genes, SVM-RFE-candidate signature genes, and XGboost-candidate signature genes. As shown in Figure 5D, five overlapping genes were identified: ASAP1-IT1, IKZF2, KLHL14, lnc-C10orf131-1, and LOC101927805. These five genes were considered as real signature genes.

Figure 5 Identification of the signature genes. (A) LASSO regression analysis. (B) ROC analysis of 12 candidate signature genes in the LASSO regression model. (C) SVM-RFE algorithm. (D) ROC analysis of 68 candidate signature genes in XGBoost model. (E) Overlap analysis between the LASSO-candidate, SVM-RFE-candidate, and XGBoost-candidate. Abbreviations: LASSO, last absolute convergence and selection operator; SVM-RFE, support vector machine-recursive feature elimination, XGBoost, extreme gradient boosting decision trees; ROC, receiver operating characteristic.

Next, we evaluated the ability of the five signature genes to distinguish between the IDD samples and healthy control samples in the GSE124272 and GSE150408 datasets. In the GSE150408 dataset, ASAP1-IT1 had the highest accuracy (AUC =0.9), followed by IKZF2 (AUC =0.789), LOC101927805 (AUC =0.751), KLHL14 (AUC =0.74), and lnc-C10orf131-1 (AUC =0.706) (Figure 6A). In the GSE124272 dataset, the ROC curves of all five signature genes had AUC >0.8 (Figure 6B). We then constructed a linear regression model based on the five signature genes, and the ROC curves showed that the 5-gene model had an AUC of 1 in the GSE150408 dataset (Figure 6C) and an AUC of 0.953 in the GSE124272 dataset (Figure 6D). Furthermore, we assessed the differences in expression of the five signature genes between the IDD samples and healthy control samples. Specifically, IKZF2, KLHL14, and LOC101927805 were significantly downregulated in the IDD samples compared to the healthy control samples, whereas lnc-C10orf131-1 and ASAP1-IT1 were significantly upregulated (all P <0.05; Figure 6E-F). This suggested that the five signature genes (ASAP1-IT1, IKZF2, KLHL14, lnc-C10orf131-1, and LOC101927805) could serve as potential diagnostic markers for IDD.

Figure 6 Evaluation of the signature genes. (A-B) ROC curve for the diagnostic markers in the GSE150408 (A) and GSE124272 (B) datasets. (C-D) Logistic regression of ROC curve in the GSE150408 (C) and GSE124272 (D) datasets. (E-F) Wilcoxon rank-sum test of diagnostic markers in the GSE150408 (E) and GSE124272 (F) datasets. * P <0.05, ** P <0.01, **** P <0.0001. Abbreviations: ROC, receiver operating characteristic; GSE, Gene Expression Omnibus Series.

Exploring the Potential Molecular Mechanisms of the 5 Signature Genes

We conducted to examine the potential biological pathways of the five signature genes. Collectively, we found that the five signature genes were significantly associated with the immune response, inflammatory response, and cell cycle-related signaling pathways, including the adaptive immune response, regulation of inflammatory response, and negative regulation of the cell cycle. Supplementary Figure S2 shows the top 10 GO terms enriched by ASAP1-IT1 (Figure S2A), IKZF2 (Figure S2B), KLHL14 (Figure S2C), lnc-C10orf131-1 (Figure S2D), and LOC101927805 (Figure S2E). Kyoto Encyclopedia of Genes and Genomes (KEGG) analysis was also performed using GSEA to reveal the pathways in which the individual signature genes were involved (Supplementary Figure S3). A comprehensive analysis revealed that these signature genes correlated with the immune-inflammatory response, amino acid degradation and metabolism, cell cycle, and multiple disease-related signaling pathways, including antigen processing and presentation, valine leucine and isoleucine degradation, DNA replication, and acute myeloid leukemia-related signaling pathways. Supplementary Figure S2 shows the top 10 pathways, including ASAP1-IT1 (Supplementary Figure S3A), IKZF2 (Supplementary Figure S3B), KLHL14 (Supplementary Figure S3C), lnc-C10orf131-1 (Supplementary Figure S3D), and LOC101927805 (Supplementary Figure S3E).

Co-Expression Network of DE-lncRNA-Signature Genes

After matching the lncRNA expression profiles in the GSE150408 dataset, DE analysis was performed to identify the DE-lncRNAs. We identified 37 DE-lncRNAs and calculated their Pearson’s correlations with the five signature genes. Based on the lncRNA-mRNA regulatory relationship theory, we identified lncRNA-mRNA relationship pairs that had a cor value >0.3 and a P value <0.05. Ultimately, we identified 38 lncRNA-mRNA relationship pairs and visualized the DE-lncRNA-signature gene co-expression network using the Cytoscape software, which contained 29 nodes (5 signature genes and 24 lncRNAs) and 38 edges (Supplementary Figure S4). KLHL14, LOC101927805, and lnc-C10orf131-1 were simultaneously regulated by LINC01560. IFNG-AS1 and PLCG1-AS1 were regulated by IKZF2, KLHL14, and LOC101927805. lnc-C10orf131-1 and ASAP1-IT1 are simultaneously regulated by six lncRNAs. Both FUT8-AS1 and LINC00920 regulate IKZF2 and KLHL14 expression.

Expression Level of Key Differentially Expressed Genes (DEGs)

To further validate the correctness of the analyzed outcomes, the expression of the three DEGs was identified using RT-qPCR. The results showed that three DEGs associated with IDD were upregulated in patients with IDD. Furthermore, the results demonstrated that IKZF2 and KLHL14 were significantly downregulated in IDD patients, while ASAP1-IT1 was significantly upregulated (Figure 7A-C).

Figure 7 Expression level of 3 hub DEGs in healthy control and IDD patients. (A-C) Expression levels of ASAP1-IT1 (A), KLHL14 (B) and IKZF2 (C). N=20, ** P <0.01, *** P <0.001. Abbreviations: DEGs, differentially expressed genes; IDD, intervertebral disc degeneration.

Discussion

Lumbar disc herniation is often caused by lumbar disc degeneration,²³ mechanical overload,²⁴ or trauma.²⁵ When the IVD protrudes, the physical barrier between the IVD and the immune system is broken, and the nucleus pulposus is exposed to the immune system, which in turn causes the expression of chemical and inflammatory factors and infiltration of inflammatory cells.⁷ Inflammatory cytokines are secreted by the nucleus pulposus, annulus fibrosus, and infiltrating immune cells around herniated IVD. By activating matrix metalloproteinases and other proteases, cytokines cause the decomposition and destruction of the extracellular matrix and increase the recruitment of immune cells, thus forming an immune and inflammatory microenvironment around the IVD.²⁶

There are few reports about natural killer cells in IVD herniations. It has been reported that the activity of natural killer cells is more downregulated in patients with IDD than in healthy individuals.²⁷ In addition to being recognized by macrophages, natural killer cells also have positive cytotoxic effects on nucleus pulposus cells, which may have an immunological function in the early phase of IDD.¹⁹ Another bioinformatics analysis reported that the DEGs of IDD patients were involved in natural killer cell-mediated cytotoxicity, which is consistent with our findings.²⁸

Neutrophils have been shown to release proteolytic enzymes to accelerate the degenerative process of the lumbar disc, especially in IDD patients who are cigarette smokers.²⁹ Gorth et al verified that cytokine expression induced by IDD in a mouse model could result in neutrophil infiltration into the disc, which accelerates the process of IDD.^30,31 Two recently published WGCNA co-expression analyses revealed that neutrophil activation, regulation of the cytokine-mediated signaling pathway, and apoptosis play important roles in the IDD process, which is consistent with our findings.^32,33

Several studies have contended that ncRNAs, including microRNAs and lncRNAs, play significant roles in IDD.^34–36 In this study, we found that some lncRNAs were closely related to signature genes. Previous studies have shown that these genes and lncRNAs are involved in the development and progression of several cancers. For example, lncRNA FUT8-AS1 has been shown to inhibit melanoma malignancy by promoting miR-145-5p biogenesis and inhibiting NRAS/MAPK signaling.³⁷ INC00920, which is regulated by ERG, has been shown to increase prostate cancer cell survival via the 14-3-3ε-FOXO pathway.³⁸ LINC00174 has been shown to increase glycolysis and tumor progression by regulating the miR-152-3p/SLC2A1 axis in glioma.³⁹ SGMS1-AS1 has been shown to restrain the development of lung adenocarcinoma via the miR-106a-5p/MYLI9 axis.⁴⁰ TMEM92-AS1 stimulates gastric cancer progression by binding to YBX1 to mediate CCL5.⁴¹ TPRG1-AS1 suppresses liver cancer progression by sponging miR-4691-5p and miR-3659.⁴² FAM201A has been shown to have a protective effect against interleukin (IL)-1 β-induced chondrocyte injury.⁴³ However, the roles of FUT8-AS1, INC00920, LINC00174, TMEM92-AS1, TPRG1-AS1, and FAM201A in IDD have not yet been reported. As the first well-described lncRNA that enhances interferon (IFN)-γ production in T cells, IFNG-AS1 plays an essential role in TH1 polarization and T-cell effector functions.⁴⁴ We speculate that IKZF2, KLHL14, and LOC101927805 regulate IDD by regulating T-cell immunity through IFNG-AS1. However, comprehensive and integrated analyses of lncRNAs associated with IDD are lacking.

This study identified signature genes IKZF2, LOC101927805, KLHL14, lnc-C10orf131-1, and ASAP1-IT1. In this study, GSEA results suggested that these five signature genes were mainly involved in NADPH, IgG binding, and immune receptor activity. NADPH oxidase activation and IgG binding have been shown to be associated with IDD.^45,46 Collectively, we found that these five signature genes were significantly associated with the immune system, immune-inflammatory responses, skeletal system, cell cycle, amino acid degradation, and metabolism. Thus, we speculated that these five signature genes play crucial roles in IDD by regulating immune feedback and the skeletal system.

Among the hematopoiesis-specific transcription factors that regulate lymphocyte development, IKZF2 belongs to the Ikaros family of zinc-finger proteins. IKZF2 has been widely studied in leukemia and lymphoma, but has seldom been studied in IDD. ZFP36L2 reduces IKZF2 expression and inhibits regulatory T-cell function, indicating that IKZF2 is involved in the post-transcriptional regulation of immunological responses and inflammatory diseases.⁴⁷ In GSEA, IKZF2 was significantly associated with mitochondrial gene expression, oxidative phosphorylation, and the cytokine-mediated signaling pathway, which participates in IDD.^32,48,49

A previous study examined the role of Klhl14 in the IVD and axial skeleton development processes for the first time.⁵⁰ According to our study, Klhl14 is significantly correlated with ribosome and B cell activation, which are implicated in the process of IDD.⁵¹ ASAP1-IT1 expression was associated with the proportion of neutrophils and gamma delta T cells, which was beneficial for regulating the infiltration of immune cells and contributing to the diagnosis and treatment of IDD.³³ As for LOC101927805 and lnc-C10orf131-1, no relevant studies have been conducted. Additionally, ASAP1-IT1 and LOC101927805 were found to be associated with the adaptive immune response and DNA repair in GSEA analysis, which is involved in the process of IDD.^52,53 Thus, we are the first to find that LOC101927805 and lnc-C10orf131-1 may influence the occurrence and progression of IDD and could potentially be useful in diagnosing IDD.

Molecular markers serve as indicators of biological alterations. The shifts observed in the biomarker spectra within degenerated intervertebral discs and peripheral blood could potentially hold a significant sway in the early identification of IDD).⁵⁴ Diagnostic value analyses demonstrated that the identified biomarkers exhibit high precision and stability in discerning IDD. Nevertheless, clinical validation is essential to establish these genes as viable therapeutic targets or diagnostic markers. Once clinical verification is complete, the expression patterns of the central genes can be harnessed to formulate pertinent diagnostic markers or biomarkers, which will be instrumental in diagnosing the disease, forecasting its course, and gauging treatment efficacy. On the other hand, by capitalizing on the specific expression patterns of these central genes, drugs can be developed to directly act on them, thereby modulating the immune response and impeding the progression of IDD.

There are also some limitations in the current study, including the small sample size and incomplete clinical parameters of IDD patients. In addition, several studies have focused on biomarkers of IDD involved in autophagy, aging, and lymphatics.^55–57. However, we examined the expression and diagnostic value of these hub genes using different machine-learning algorithms and intend to apply more in vitro experiments as well as clinical work to further determine the findings of this study in the future. In addition, the predicted lncRNA-mRNA interactions should be further examined to determine the underlying regulatory mechanism of IDD through subsequent basic experiments. In addition, we will predict potential therapeutic drugs through a drug screening system and verify them at the animal level.

Conclusions

We identified five signature genes and constructed an lncRNA signature gene network to reveal the regulatory mechanisms of the biomarkers. Furthermore, we verified the three signature genes using biological experiments. These verified hub genes may provide important insights into the potential molecular mechanisms of IDD and may reveal new biomarkers for diagnosis.