In the Gene Expression Omnibus (GEO), with intervertebral disc degeneration and nucleus pulposus cells as the retrieval conditions, the species was set as human, and finally two datasets (GSE56081 and GSE70362) were selected for analysis. GSE56081 contains five nucleus pulposus samples from IDD patients and five normal nucleus pulposus samples. GSE70362 included samples from 10 IDD patients and 14 normal nucleus pulposus samples. The two datasets were normalized using the “preprocessCore” package. In addition, cell classification inference of immune cell composition was performed in both the initial and validation datasets, allowing the immunomodulatory effects of key IDD biomarkers to be compared across different IDD samples. The ferroptosis-related genes (FRGs) were derived from FerrDb database. The flowchart of the analysis process was summarized in Figure 1.

The “preprocessCore” packages in R software were used to avoid batch effects, and the normalization of both datasets was validated by boxplots with the ggpolt2 package. Differential genes (DEGs) were determined using the limma software package analysis, with |log2 (fold change)|>1 and p < 0.05 as DEGs with statistically significant differences. Volcano plots and heatmaps were generated using “ggplot2” and the “pheatmap” packages, showing statistically significant DEGs. Intersection of FRGs with DEGs using the Venn diagram package and showing the expression levels of DE-FRGs by drawing heatmap.

To further reveal the important functions and biological pathways of co-expressed differentially genes, we performed GO function and KEGG pathway enrichment analysis on them. The “clusterProfiler”, “org.Hs.e.g.db”, and “ggplot2” packages in the R software were used to analyze and visualize the GO and KEGG pathway enrichment of differential genes. Setting the minimum gene set to five and the maximum gene set to 5000, p values of <0.05 and FDR of <0.25 were considered statistically significant.

In order to further identify key genes from differentially expressed ferroptosis-related genes, the Lasso regression algorithm and the SVM-RFE support vector machine recursive feature elimination algorithm were used to screen eigengenes. The gene results were intersected, and the R software “venneuler” package was used to draw a Venn diagram to visualize the results.

For the eigengenes obtained by the above intersection, the value of the obtained genes as diagnostic markers was verified in two independent data sets, GSE56081 and GSE70362, and the diagnostic value was evaluated by drawing the receiver operating characteristic (ROC) curve. With p < 0.05 is the threshold to determine.

In order to evaluate the infiltration of immune cells in IDD and the relationship between the gene expression of the screened IDD diagnostic biomarkers and the infiltration of immune cells in the intervertebral disc tissue, the CIBERSORT algorithm (https://cibersort.stanford.edu/) was used to analyze the GSE56081 and GSE70362, and the appropriate samples were screened according to p < 0.05 and the percentage of each type of immune cell in the sample was calculated. Wilcox test was used to compare the differences between the two groups, and the correlation between the infiltration rates of various types of immune cells was determined through the corrplot package.

This study used NP tissue, which had the informed consent of the patient’s family and was approved by the Medical Ethics Committee of Fuyang People’s Hospital (Fuyang, China). Ethical approval number is Review of Medical Ethics [2022]2. During the period from January 2022 to August 2022, the degenerative samples were taken from the degenerative NP cell samples of patients with degenerative intervertebral disc disease who underwent discectomy in Fuyang People’s Hospital, and the control samples were taken from the patients who had undergone surgery for idiopathic scoliosis. After cell collection, NP cells were resuspended in DMEM/F12 medium containing 10% FBS for culture. Observe the cell growth and adherence every day, and pass it to the third generation for relevant experiments. The clinical characteristics of the patients are shown in Table 1.

The primer sequence of each gene was shown in Table 2. A total of 10 NP tissues were used, of which 5 were degenerative and five were normal. Nucleus pulposus cells are digested by trypsin and treated with the TriZol reagent (Invitrogen) to extract total RNA after ultrasonic fragmentation. The PrimeScript RT-PCR kit (Takara) was used to synthesize cDNA. The 7500 real-time fluorescent PCR system (Thermo Fisher) was used for RT-PCR, and GAPDH was used as the endogenous control.

The expression of DE-FRGs was evaluated by immunohistochemistry of sections of nucleus pulposus with different levels of degeneration. NP tissue was fixed with 4% paraformaldehyde, then paraffin embedded and sectioned. Histological analysis was performed. Sections were desaffinified, rehydrated, and stained with Eosin (HE), Masson, and Alcian blue, respectively. Immunohistochemistry was performed according to the aforementioned methods (Yi et al., 2020). Sections were incubated with a primary antibody (dilution 1:200) resistant to ENPP2 (14243-1-AP, Wuhan proteintech), NOX4 (14347-1-AP, Wuhan proteintech), FADS2 (680261-LG, Wuhan proteintech), and TFAP2A (67076-1-lg, Wuhan proteintech). Next, the slices were incubated with the secondary antibody (pv6000, Beijing Zhongshan). Buffer was used instead of primary antibody as negative control. Finally, three fields of each slide were randomly selected and observed with a microscope (Olympus, Tokyo, Japan).

Analyze the statistical data obtained from GEO database through R-3.6.1. Use GraphPad Prism 8 and SPAS software to process other data. p < 0.05 was considered statistically significant.

The data sets GSE56081 and GSE70362 were homogenized using the “preprocesscore” package (Supplementary Figures S1A, B). Next, we used the “limma” package to analyze the differences between the two datasets. A total of 339 differential genes were screened, including 182 upregulated genes and 157 downregulated genes. Figures 2A, C show the volcano map of differential genes analysis, while Figures 2B, D show the analysis heat map of differential genes. To explore the biological function of the differential genes related to IDD, we performed GO and KEGG enrichment analysis. The results of GO analysis showed that these differential genes were mainly related to neutrophil mediated immunity, neutrophil activation involved in immune response, mitotic nuclear division and so on (Figure 3C). In addition, the KEGG enrichment pathway analysis shows multiple important signaling pathways such as pathways neurodegeneration-multiple.

In order to explore the relationship between ferroptosis and intervertebral disc degeneration, we crossed the previously obtained differential genes with FRGs and obtained a total of 16 DE-FRGs (Figure 4A). Supplementary Figures S2 shows the volcanic map and heatmap of ferroptosis related differential genes in GSE56081 and GSE70362 separately. The enrichment results of GO and KEGG showed that the DE-FRGs were closely related to cellular aldehyde metabolic process, oxidoreductase activity and chemical carcinogenesis-reactive oxygen species (Figures 4B, C). In addition, we use node graph to show the relationship between KEGG results of top5 and related differential genes (Figure 4D).

Next, we analyzed the differential expression of 16 DE-FRGs in two data sets. As shown in Figures 5A, B genes were upregulated, and nine genes were downregulated in two databases. In order to further identify the key genes, we further screen the characteristic genes through machine learning. Seven genes were screened by LASSO analysis and 16 genes were screened by SVM algorithm. Seven feature genes were obtained through LASSO algorithm and SVM algorithm, and 16 feature genes are obtained through SVM-RFE algorithm (Figures 5C, D). The genes obtained by the two methods were intersected to obtain seven characteristic marker genes, and the Wayne map was drawn (Figure 5E). Finally, heat map is used to display the logFC values of seven genes in the difference analysis results of two data sets, including two highly expressed genes and five low expressed genes (Figure 5F).

Next, we use the “circle” package to analyze the correlation of seven core genes in the two data sets. As shown in Figures 6A, B, NOX4 and PIR were negatively correlated with other genes, while TFAP2A, ATF3, ENPP2, FADS2 and TIMM9 were positively correlated with other genes. We drawed ROC curves and evaluate the diagnostic value of key characteristic genes by calculating AUC values. ROC curve results showed that ENPP2 in two data sets AUC values of are 0.96 and 0.914 respectively, PIR in two data AUC values of the set are 1.0 and 0.8 respectively, TIMM9 in two data AUC values of the set were 1.0 and 0.85 respectively. The results show that they have high diagnostic value in IDD (Figures 6C, D).

In order to explore the connection between IDD and immune infiltration, we used the ssGSEA function of “GSVA” package to evaluate the degree of immune cell infiltration in GSE56081 dataset. Figure 7A showed the correlation between immune cells in IDD, macrophages and natural killer cells were positively correlated with many other immune cells including activated CD4 T cells and CD8 T cells. Figure 7B showed the difference in the infiltration of immune cells in normal intervertebral disc NP cells and degenerative intervertebral disc NP cells. The activated CD4 T cells, CD8 T cells, dendritic cells and so on were increased in degenerated NP cells, while eosinophil and type 2T helper cell were decreased in degenerated NP cells. Next, we use Cibersort algorithm to analyze the relationship between the expression of seven key DE-FRGs in GSE56081 database and immune cell infiltration (Figure 7C). The upregulated expression of NOX4 and PIR in IDD was positively correlated with multiple immune cells, while the downregulated expression of TIMM9, ATF3, ENPP2, FADS2 and TFAP2A were negatively correlated with multiple immune cells. These results suggest that high immune infiltration may play an important role in intervertebral disc degeneration.

In order to further clarify the role of DE-FRGs in IDD, we conducted gene set enrichment analysis (GSEA) function enrichment analysis. Use the data in GSE70362 to analyze the correlation between seven key genes and all genes and use the heat map to display the positive correlation top50 gene respectively (Figure 8A). Based on the results of correlation analysis, we used the “clusterProfiler” package for single gene GSEA analysis of Reactome (Figure 8B). The enrichment results indicate that ATF3 and PIR may be related to the regulation of mitochondrial related genes. ENPP2, NOX4 and TIMM9 may be related to RNA transcription and immune system, while FADS2 and TFAP2A may participate in the IDD process by regulating the cycle of NP cells.

Next, we used RT-PCR and immunohistochemical staining to confirm whether key DE-FRGs can be applied non-selectively to IDD patients. RT-PCR experiment results showed that differences in the expression levels of ENPP2, NOX4, FADS2 and TFAP2A genes were statistically significant. ATF3, PIR and TIMM9 showed no statistically significant difference. The expression of NOX4 was upregulated, whereas the expression of ENPP2, FADS2 and TFAP2A was downregulated (Figure 9). We also confirmed the above results by immunohistochemistry (Figure 10). These results were consistent with the bioinformatics analysis. In addition, the MRI images of patients with different levels of degeneration were shown in Supplementary Figure S3A, and HE, Masson and Alxin blue staining of nucleus pulposus were shown in Supplementary Figure S3B.