Introduction

Intervertebral disc degeneration (IVDD) is a prevalent age-related disease that affects a significant proportion of the adult population.¹ It is considered the main cause of lower back pain, a painful condition that imposes a substantial burden on individuals and healthcare systems worldwide.² The degeneration of the nucleus pulposus (NP) and annulus fibrosus (AF), the two main components of intervertebral disc, play crucial roles in the pathogenesis and progression of IVDD, with complex molecular mechanisms including inflammation and oxidative stress.^3,4 Understanding the intricate mechanisms underlying NP and AF degeneration is of great importance for the development of effective therapeutic strategies of IVDD.

Researchers have recognized that reduced intervertebral disc cell numbers or cell death and extracellular matrix (ECM) degradation are pivotal processes in NP and AF degeneration.^5,6 Programmed cell death (PCD), encompassing apoptosis, autophagy, and ferroptosis, were considered significant contributors to intervertebral disc cell death and IVDD.^7–9 Among them, ferroptosis has garnered increasing attention due to its involvement in various degenerative disorders.^10–12 Ferroptosis is characterized by iron-dependent lipid peroxidation, accumulation of reactive oxygen species (ROS), and mitochondrial morphological alteration.¹³ Previous research summarized that ferroptosis could be regulated by a set of genes, including drivers (promoting ferroptosis, eg nuclear receptor coactivator 4, NCOA4) and suppressors (inhibiting ferroptosis, eg ferroportin1, FPN1).¹⁴

Recently, growing evidence has elaborated the role of ferroptosis induced by oxidative stress or inflammation in the degeneration of human NP tissues and cells, highlighting its significance in IVDD.¹⁵ Targeting ferroptosis in NP cells has emerged as a promising therapeutic approach for IVDD treatment.¹⁶ The upregulated methylation of glutathione peroxidase 4 (GPX4), a vital regulator and marker of ferroptosis, contributes to the exacerbation of ferroptosis in NP cells by downregulating GPX4 expression.¹⁷ Additionally, our previous findings elucidated a critical mechanism involving FPN1 hypermethylation, which leads to ferroptosis and degeneration of NP cells.¹⁸ Notably, glutamine has been shown to suppress oxidative stress-induced ferroptosis and mitigate ECM degradation in NP cells by deubiquitinating NRF2 and inhibiting lipid oxidation.¹⁹ Nevertheless, a comprehensive understanding of the molecular mechanisms and regulatory networks of ferroptosis in NP degeneration necessitates further investigation.

On the other hand, researchers rarely dealt with the molecular regulatory network of ferroptosis in degenerated AF. Yang et al²⁰ demonstrated the occurrence of ferroptosis in tert-butyl hydroperoxide (TBHP)-treated AF cells and revealed that NCOA4-mediated ferritin selective autophagy, known as ferritinophagy, activated ferroptosis in AF cells. Furthermore, Han et al²¹ reported that removing intracellular ROS promotes the repair of AF tissues, preventing intervertebral disc herniation. Huo et al²² found that MGST1 is involved in regulating the ferroptosis phenotype of AF tissues. Revealing the key ferroptosis-related genes in degenerated AF holds significant promise for the identification of potential therapeutic targets or diagnostic markers for IVDD.

Despite the rising evidence emphasizing the correlation between ferroptosis and IVDD, the detailed role of ferroptosis in NP or AF degeneration remains elusive. Further evidence-based results are needed to reveal the regulatory network of ferroptosis during IVDD. In this study, we conducted a comprehensive analysis of gene expression profiles containing both normal and degenerated NP and AF tissues via bioinformatics analysis and experimental verification. This study aimed to identify ferroptosis-related biomarkers in IVDD, elucidate their important biological functions, perform immunoinfiltration analysis, and predict their associated drugs or ceRNAs to define their mechanisms in the degeneration process of NP and AF further. These findings could provide a foundation for further exploration of the pathological mechanisms underlying IVDD and offer potential drug targets for intervention. The workflow of the current study is illustrated in Figure 1, and the important abbreviations are shown in Supplementary Table S1.

Figure 1 Flowchart of the study. Abbreviations: GO, Gene Ontology; KEGG, Kyoto Encyclopedia of Genes and Genomes; LASSO, Least Absolute Shrinkage and Selection Operator; SVM-RFE, Support Vector Machine-Recursive Feature Elimination; ROC, receiver operator characteristic curve; GSEA, Gene Set Enrichment Analysis; GSVA, Gene Set Variation Analysis.

Materials and Methods

Data Selection

Six gene expression profiles were downloaded from the Gene Expression Omnibus (GEO) database (http://www.ncbi.nlm.nih.gov/geo).²³ We merged the GSE70362 and GSE34095 datasets and used the “Combat” algorithm in the R package “SVA” (version 3.29.1) to remove batch effects for identifying FRDEGs in NP. The GSE70362 dataset was considered as a training set for the analysis of AF tissue in this study. The GSE167199, GSE186542, and GSE244889 were used as validation analysis datasets for NP. The GSE147383 was used to verify the expression of the marker genes in AF. Detailed information about the six datasets was shown in Table 1. A total of 484 ferroptosis-related genes had been acquired from FerrDb V2 (http://www.zhounan.org/ferrdb/).²⁴

Table 1 Detailed Information for All the Relevant Datasets

Identification of FRDEGs in Degenerated NP and AF Tissues

All the data that was downloaded from GEO database were imported into R software version 4.2.1. Quality control (QC) analysis, such as principal component analysis (PCA), was used to determine the sample distribution. We first extracted the expression data of ferroptosis-related genes in non-degenerated and degenerated NP (345 genes) and AF (287 genes) tissues from the GSE70362 and GSE34095 datasets. The “limma” package was used to identify the FRDEGs, and the “heatmap” package was used to draw heatmaps, showing the genes with P value <0.05. The logFC value was set at 1.0 and 0.2 for NP and AF tissues, respectively. The correlations among the FRDEGs in degenerated NP and AF were analyzed and visualized via “corrplot” package.

Functional and Pathway Enrichment Analysis

Gene Ontology (GO) and the Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis were performed using the R “clusterProfilter” package to further investigate the potential biological functions of the FRDEGs in NP and AF degeneration, with a P value of <0.05 indicating statistically significant. We used “ggplot2” package to visualize the important results of function and pathway enrichment analysis. Three terms including biological process (BP), molecular function (MF), and cellular component (CC) were conducted for GO enrichment.

Identification of Feature FRDEGs via Algorithms

The least absolute shrinkage and selection operator (LASSO) algorithm was applied with the “glmnet” package to reduce the dimensions of the data.²⁵ The expression data of FRDEGs were retained for feature selection with 10-fold cross-verification. We utilized a novel machine learning technique, support vector machine-recursive feature elimination (SVM-RFE) algorithm, to rank features of FRDEGs via “e1071” package, which was compared by the average misjudgement rates of their 10‐fold cross‐validations.²⁶ Then, feature FRDEGs for NP and AF degeneration were identified by overlapping biomarkers derived from the two algorithms and were assessed by calculating the receiver operating characteristic (ROC) curve.

Gene Set Enrichment Analysis (GSEA) of Hub FRDEGs

The annotation gene set “c2.cp.kegg_legacy.v2023.2.Hs.symbols.gmt” was downloaded from the Molecular Signatures Database (MsigDB).²⁷ All genes were sorted according to their correlations from high to bottom after correlation analysis between the feature FRDEGs and all other genes, which were used for GSEA analysis via R “GSEA” package. Moreover, P < 0.05 was considered significantly enriched. R “enrichplot” package was used for visualization.

Gene Set Variation Analysis (GSVA) of Hub FRDEGs

The “c5.go.v2023.2.Hs.symbols.gmt” set downloaded from MsigDB was regarded as the reference set to perform GSVA analysis on each feature gene. Based on the expression data of feature FRDEGs, samples were classified into high- and low-expression groups. We applied the “limma” package to analyze the difference in GSVA scores, with |t| > 2 and P < 0.05 consider significant. Terms with t > 0 and t < 0 were considered to be activated in the high‐expression group and the low‐expression group, respectively.

Immune Infiltration Analysis

The relative contents of 22 types of infiltrating immune cells were calculated by cell-type identification by estimating relative subsets of RNA transcripts (CIBERSORT) algorithm.²⁸ The contents of 22 immune cells were compared between the normal and degenerated samples, showing in a violin diagram by “vioplot” package. Furthermore, Spearman correlation analysis between the expression levels of feature FRDEGs and the contents of the immune cells was conducted, which was visualized in a correlation heatmap via “ggplot2” package.

Drug Regulatory Network and Competitive Endogenous RNAs (ceRNA) Network Construction

Drugs targeting feature genes were screened using the Drug–Gene Interaction database (DGIdb, version 4.0, www.dgidb.org).²⁹ The Cytoscape (version 3.9.1) was applied to perform the interaction network between potential drug and feature genes. Three databases were used to predict mRNA‐microRNA (miRNA) interaction pairs based on feature genes, including miRanda, miRDB, and TargetScan. Then, we predicted lncRNA of the miRNAs in starBase to obtain the network of miRNA‐lncRNA. Cytoscape was used to construct a ceRNA network.

Validation of the Feature FRDEGs by Other GEO Datasets

The GSE167199, GSE186542, and GSE244889 were merged and used as validation analysis datasets for NP. The GSE147383 was used to verify the expression of the hub genes in AF. “limma” package in R software was used to identify the difference between the two groups, and “ggplot2” package was used for visualization. Statistical significance was set at P < 0.05.

Validation of the Feature FRDEGs by Single-Cell RNA Sequencing (scRNA-Seq) Data

GSE230808 raw data was downloaded. We used R “Seurat” package to conduct data process, including normalization, scaling, clustering, and identified the NPCs and AF cells from the clusters as previously described.^2,30,31 Clustering was performed with the FindClusters functions in Seurat using Clustree analysis (resolution from 0.1 to 1). The FindAllMarkers function was used to identify marker genes of each cluster, and then NP cells and AF cells were screened. Monocle pseudotime trajectory was used for identifying the NP cell and AF cells from early stage to late stage. To identify differentially expressed genes among clusters with interest, the function FindMarkers was used, with abs(logFC) > 0.5 and P value < 0.05 indicating significant.

Human Primary NP Cells and AF Cells Culture

Human primary NP cells (ScienCell, Catalog#4800, USA) and AF cells (ScienCell, Catalog#4810, USA) were purchased from ScienCell Research Laboratories according to our previous research.¹⁸ The cells were cultured with complete culture medium (DMEM, Gibco, Invitrogen, USA) containing 10% fetal bovine serum (FBS, Gibco, Invitrogen, USA) in a humidified atmosphere of 5% CO2 at 37°C. The degenerated cell models were constructed by using tert‐Butyl hydroperoxide (TBHP, 50μM) for 24 hours, then collected for subsequent experiments.

IVDD Rat Model Construction

The puncture rat model was developed as previously described.¹⁸ Briefly, 12 Male SD rats (3 months old) were randomly divided into control and puncture (degenerated) groups. After anesthesia by intraperitoneally injecting chloral hydrate (5%, 0.7mL/100g), and sterilized, a 31-gauge sterile needle was used to puncture the tail vertebral disc Co6/7 from the lateral side of the tail. All the SD rats were allowed unrestricted activity. After eight weeks, the rats were euthanized by intraperitoneal injection of excessive anesthetic, and the tail vertebra specimens were collected. Specimens were soaked in 4% paraformaldehyde for 1–2 d, followed by decalcified for 8 weeks via 10% Ethylenediaminetetraacetic acid (EDTA) at room temperature.

Cell Immunofluorescence

Cells were fixed with 4% paraformaldehyde solution, then permeabilized with 0.1% Triton-X100 solution for 30 min. After blocking with 1% goat serum, the cells were incubated overnight at 4°C with rabbit-derived MUC1 (Abcam, ab109185) and CREB3 (Proteintech, 11275-1-AP) antibody (diluted in 1% BSA-PBS at a dilution ratio of 1:200). Secondary antibody Alexa Fluor 594-conjugated goat anti-rabbit IgG was incubated for 1 hour. DAPI and phalloidin staining was performed for 20 minutes. Fluorescence microscope (Nikon, Japan) was used to obtain images, and the intensity was analyzed by Image J.

Western Blotting (WB)

Total proteins were extracted using RIPA lysis buffer with protease inhibitors and phosphatase inhibitors (Beyotime, Nantong, China). Tissues needed liquid nitrogen freezing and grinding. SDS–PAGE gel electrophoresis was used to separate proteins, followed by transferring to polyvinylidene fluoride (PVDF) membranes, and then blocked for 120 min. PVDF membranes were cut into ideal size according to the protein molecular weights and incubated overnight at 4°C in primary antibody, including GPX4 (Abcam, ab125066; 1:5000) and Tubulin (Proteintech, 10068-1-AP; 1:1000). After being washed by TBST, these membranes were put into secondary antibody Goat Anti-Rabbit IgG H&L (HRP) (Abcam, ab6721; 1:5000) for 60 min. Image J was used to quantify the density of the blots.

Ferrous Ion (Fe²⁺)

A FerroOrange fluorescent probe (DOJINDO, Japan) was used to detect intracellular Fe²⁺ according to the manufacturer’s instructions. In brief, AF cells were incubated with 1.0 µM FerroOrange for 30 min at 37 °C in a dark environment. After washing with PBS, fluorescence microscope (Nikon, Japan) was used to obtain images. The intensity was quantified by Image J.

Immunohistochemistry (IHC)

Intervertebral disc tissues (human and rats) were embedded in paraffin and sliced into 4 μm sections. For Immunohistochemical examination, the primary antibodies (MUC1 and CREB3) were added to the sections at a dilution of 1:200 in this assay and incubated overnight at 4°C, then incubated with secondary antibody. Disc sections were observed under an optical microscope (Olympus BX51) and analyzed with Image-Pro Plus software (Version 5.1.0.20. Media Cybernetics, Inc). The staining intensity was analyzed and used for comparison between groups.

Ros

ROS for AF cells was measured with 2′7′-dichlorodihydro-fluorescein diacetate (DCFH-DA) assay kit (Beyotime, China) according to the manufacturer’s instructions. Briefly, the cells in 6-well plates were incubated with DCFH-DA (10 μM/L) at 37°C for 20–30 min. Photos were taken using a fluorescence microscope (Nikon, Japan) and analyzed by Image J.

Transmission Electron Microscopy (TEM)

The TEM for AF cells was conducted as we previously described.¹⁸ Briefly, cells were fixed with 2.5% glutaraldehyde solution and 1% osmium tetroxide solution. Then, we use ethanol solutions to treat cells, followed by pure acetone. The embedded samples were sliced into 70–90 nm sections via an ultramicrotome (EM UC7, LEICA). The sections were stained with lead citrate and uranyl acetate and observed by a transmission electron microscope (JEM-1200EX, Japan).

Statistical Analysis

R of version 4.2.1 was used to conduct most statistical analysis in this study. The correlation analysis among the FRDEGs was performed by Pearson analysis. Data are presented as the mean ± Standard Deviation of triplicate-independent experiments, and GraphPad Software (Version 8.0.2, San Diego, CA, USA) was used for the statistical analysis via Student’s t-tests, one-way analysis of variance (ANOVA). P < 0.05 indicates statistical significance.

Results

Identification of the FRDEGs in Degenerated NP and AF Tissues

According to the QC results, two samples of NP tissue (GSM1725783 and GSM1725787) in GSE70362 were excluded, and two samples of AF tissue (GSM1725782 and GSM1725786) in GSE70362 were excluded from analysis (Supplementary Figure S1). The PCA results after exclusion were depicted in Figure 2A and B. 15 (eight upregulated genes and seven downregulated genes) and 18 (seven upregulated genes and 11 downregulated genes) FRDEGs in degenerated NP (Table 2) and AF (Table 3) tissues were identified, respectively. The heat map illustrated the expression level of the FRDEGs in NP samples (Figure 2C) and AF samples (Figure 2D). Figure 2E showed that MUC1 was positively correlated with MT1G and BEX1, MMD was negatively correlated with MT1G and BEX1 in degenerated NP samples. For degenerated AF samples, CREB3 was positively correlated with RPL8 and DPEP1, BECN1 was negatively correlated with IDO1 (Figure 2F).

Table 2 The Ferroptosis-Related Differentially Expressed Genes (FRDEGs) in Degenerated Nucleus Pulposus

Table 3 The Ferroptosis-Related Differentially Expressed Genes (FRDEGs) in Degenerated Annulus Fibrosus

Figure 2 The expression levels of ferroptosis-related differentially expressed genes (FRDEGs) in degenerated nucleus pulposus (NP) and annulus fibrosus (AF). (A) Principal component analysis (PCA) of the samples for NP. (B) PCA of the samples for AF. (C) Heatmap shows expression patterns of FRDEGs in NP samples from normal or patients with intervertebral disc degeneration (IDD). (D) Heatmap shows expression patterns of FRDEGs in AF samples. (E and F) The correlation of these genes in NP and AF, respectively. ***p < 0.001, **p < 0.01, *p < 0.05.

Enrichment Analysis for the FRDEGs

GO and KEGG pathway analysis were conducted to detect the important biological function of the FRDEGs. For NP, GO results showed that the FRDEGs were mainly enriched in oxidoreductase activity and peroxisome (Figure 3A and B). For AF, responses to oxidative stress and peroxisome were significantly enriched via GO analysis (Figure 3C and D). In addition, the KEGG analysis for NP depicted that the FRDEGs were significantly associated with apoptosis pathway, folate biosynthesis, and fatty acid degradation (Figure 3E). The results of KEGG for AF showed that the FRDEGs were significantly enriched in necroptosis, autophagy, apoptosis, and TNF signaling pathway (Figure 3F).

Figure 3 Gene Ontology (GO) enrichment and Kyoto Encyclopedia of Genes and Genomes (KEGG) analysis of ferroptosis-related differentially expressed genes (FRDEGs). (A and B) Bubble plot and circos diagram of enriched GO terms for nucleus pulposus (NP). (C and D) Bubble plot and circos diagram of enriched GO terms for annulus fibrosus (AF). (E) Bar plot of enriched KEGG terms for NP. (F) Bar plot of enriched KEGG terms for AF. BP, biological process; CC, cellular component; MF, molecular function.

Identification of Feature FRDEGs by LASSO and SVM-RFE

We combined LASSO analysis with SVM-RFE analysis to identify the feature genes for NP degeneration and AF degeneration from the 15 and 18 FRDEGs, respectively. For NP, four genes (AKR1C1, AKR1C3, MUC1, and ENPP2) were screened by LASSO analysis (Figure 4A and B) and 10 genes (AKR1C3, NQO1, AKR1C1, MUC1, GSTZ1, MMD, ENPP2, CA9, MT1G, and PGD) were recognized by SVM-RFE analysis (Figure 4C and D). Four feature FRDEGs were ultimately identified via a Venn diagram (Figure 4E), and the ROC curve of each gene was established for determining the diagnostic performance via R “glm” package (Figure 4F). As for AF, six genes (SCP2, ABCC1, TLR4, KLF2, IDO1, and CREB3) were selected by LASSO analysis (Figure 4G and H) and 11 genes (IDO1, KLF2, ABCC1, SCP2, PEX2, KDM5C, PTGS2, IFNA1, HMGB1, RPL8, and CREB3) were outputted by SVM-RFE analysis (Figure 4I and J). Five feature FRDEGs were ultimately identified (Figure 4K), and the ROC curve of each gene was established (Figure 4L).

Figure 4 Identification of feature genes for nucleus pulposus (NP) and annulus fibrosus (AF) degeneration by algorithm. (A and B) By Least Absolute Shrinkage and Selection Operator (LASSO) algorithm, with penalty parameter tuning conducted by 10‐fold cross‐validation for NP degeneration. (C and D) Support Vector Machine-Recursive Feature Elimination (SVM‐RFE) algorithm to identify the optimal combination of feature genes. (E) The four hub genes obtained from the LASSO and SVM‐RFE by venn diagram. (F) Receiver operator characteristic curves and the area under the curve (AUC) was calculated for the hub genes. (G-L) Similar to NP, the hub genes identification for AF degeneration by LASSO and SVM-RFE.

GSEA and GSVA Enrichment Analysis for the Feature FRDEGs

GSEA enrichment analysis was conducted to further explore the biological function of the feature FRDEGs (Supplementary Figures S2 and S3). AKR1C1, ENPP2, and MUC1 in NP were mainly enriched in apoptosis, DNA replication, proteasome, arginine and proline metabolism, and AKR1C3 was associated with steroid hormone biosynthesis. CREB3, IDO1, and SCP2 in AF were significantly correlated with arginine and proline metabolism, notch signaling pathway, steroid hormone biosynthesis, glutathione metabolism.

GSVA enrichment analysis was performed to identify differentially activated biological functions between the high‐and low‐expression groups based on the expression levels of each feature gene (Supplementary Figures S4 and S5). The feature genes in NP were mainly enriched in cellular response to hydrogen peroxide, DNA strand elongation involved in DNA replication, DNA cytosine deamination, and DNA deamination. The feature genes in AF were significantly correlated with autophagosome-lysosome fusion, cell recognition, cellular response to hydroperoxide, cell matrix adhesion, cell cycle DNA replication initiation, and negative regulation of cytokine production involved in immune response.

Immune Infiltration Analysis

Enrichment analysis demonstrated that the feature FRDEGs of NP and AF were closely related to the immune response. Based on single-cell sequencing data, we have also found an inseparable connection between immune microenvironment and IVDD.³² CIBERSORT algorithm was used for identifying different immune microenvironment between controlled and degenerated NP or AF samples. Immune microenvironment was not statistically changed in both degenerated NP tissue (Figure 5A) and AF (Figure 5C) tissue in contrast to controls. Furthermore, we analyzed the correlation between these feature FRDEGs and different types of infiltrating immune cells in degenerated samples by spearman. AKR1C1 in NP had a positive correlation with Macrophages M0 (r = 0.78, p = 0.007) and Mast cells activated (r = 0.61, p = 0.046) (Figure 5B). CREB3 in AF was correlated with Monocytes (r = 0.88, p = 0.007), Dendritic cells activated (r = 0.90, p = 0.004), and Eosinophils (r = −0.85, p = 0.008) (Figure 5D).

Figure 5 Immune landscape analysis. The immune microenvironment in degenerated nucleus pulposus (NP) (A) and annulus fibrosus (AF) (C) by CIBERSORT algorithm. Pearson correlation analysis revealed the correlation between the hub genes and immune cells in degenerated NP (B) and AF (D). **p < 0.01, *p < 0.05. Abbreviation: IDD, intervertebral disc degeneration.

Established Drug and ceRNA Networks for the Feature FRDEGs

The DGIdb database was utilized for predicting potential pharmacological targets with feature FRDEGs, and Cytoscape software was used to visualize the results (Supplementary Figure S6). A possible inhibitor INDOMETHACIN targeting AKR1C3 in degenerated NP tissue was found. A possible inhibitor (SULFINPYRAZONE) for ABCC1 and three inhibitors (EPACADOSTAT, PF-06840003, and LINRODOSTAT) for IDO1 were identified in AF tissue.

We utilized the miRanda, miRDB, and TargetScan databases to build a ceRNA network based on the feature genes. The network of the feature genes in NP comprised 71 nodes (3 feature genes, 38 miRNAs, and 30 lncRNA). The ceRNA network of ENPP2 was relatively complex, including 31 miRNAs and 27 lncRNAs. Among them, has-miR-1184 and has-miR-1238-3p could be regulated by a variety of lncRNAs (Supplementary Figure S6). The network of the feature genes in AF comprised 221 nodes (5 feature genes, 87 miRNAs, and 129 lncRNA). We found that 19 lncRNAs and 11 lncRNAs could regulate the expression of ABCC1 through competitive binding with hsa‐miR‐149-3p and hsa‐miR‐185-5p, respectively. For CREB3, six miRNAs including miR-4270, miR-2355-5p, miR-542-3p, miR-34b-3p, miR-595, and miR-1225-5p exert regulatory roles in this gene (Supplementary Figure S7).

Expression of the Feature FRDEGs in the Validation Sets

GSE167199, GSE186542, and GSE244889 (containing nine control and nine degenerated NP tissues) were used for validating the expression level of these feature FRDEGs. We merged these datasets and used R package “sva” to remove batch effects. Figure 6A showed the PCA results of all the samples. Consistently, in the merged dataset, the expression level of ferroptosis-suppressor MUC1 in the degenerated NP tissue was lower than that in the control samples (P = 0.035), which was visualized by “ggpubr” package (Figure 6B-E). We used GSE147383 (containing two normal and two degenerated AF tissues) dataset to assess the five feature genes expression (method=Wilcox.test) and found that none of these genes were significantly changed in degenerated AF (Figure 6F-J). Based on the function of these five genes in ferroptosis and their expression levels in the GSE70362 dataset, the upregulated ferroptosis driver SCP2 and IDO1, as well as the downregulated suppressor CREB3 were selected as candidate genes for verification.

Figure 6 The expression level of each hub gene in the validation set. (A) Principal component analysis (PCA) of all the nucleus pulposus (NP) samples in GSE167199, GSE186542, and GSE244889. The expression of AKR1C1 (B), AKR1C3 (C), MUC1 (D), and ENPP2 (E) in NP samples. The expression of SCP2 (F), ABCC1 (G), KLF2 (H), IDO1 (I), and CREB3 (J) in annulus fibrosus samples. *p < 0.05; ns, nonstatistical significance. Abbreviations: IDD, intervertebral disc degeneration; ns, nonstatistically significant.

Expression of the Feature FRDEGs in NP by scRNA-Seq Analysis

Two NP samples with mild degeneration (stage II–III) were used for the analysis (GSM7235341 and GSM7235344). Nine clusters (resolution = 0.6) were identified (Figure 7A and B), and after finding marker genes, the NP cells were selected as high expression of ACAN and SOX9 (cluster 1, 4, and 7) (Figure 7C). We then conducted a subcluster analysis with a resolution of 0.5 (Figure 7D and E), using the following criterion: progenitor NP cells (UBE2C and TOP2A); fibrotic NP cells (FBLN1); metabolic NP cells (DKK1); adhesive NP cells (MSMO1); and stress-responsive NP cells (CP). Fibrotic NP cells and progenitor NP cells were considered as late stage and prior stage cells, respectively. Monocle pseudotime trajectory showed the same results (Figure 7F-H). By using FindMarkers function, we found that MUC1 was downregulated in Fibrotic NP cells in contrast to progenitor NP cells (Figure 7I-K).

Figure 7 Single-cell RNC-seq analysis of the nucleus pulposus (NP) samples. (A) Cluster analysis showing clustering patterns from resolutions 0 through 1. Resolution 0.6 was used. (B) Visualization of clustering by TSNE plot using 0.6 resolution. (C) Bubble plot showing the expression of the ACAN and SOX9 related different clusters. (D) Five subclusters for the NP cells were identified using 0.5 resolution. (E) Five cell types were annotated including progenitor NP cell (Pro-NPC), fibrotic NP cell (Fibro- NPC), metabolic NP cell (Met-NPC), metabolic NP cell (Met-NPC), stress-responsive NP cell (SR-NPC). (F-H) Monocle pseudotime trajectory showing the progression of Pro- and Fibro-NPC. (I) Bubble plot showing the expression of MUC1. (J) volcano plot showing that MUC1 was downregulated in Fibro-NPC. (K) TSNE image of MUC1 in subcluster 3 and 4.

Expression of the Feature FRDEGs in AF by scRNA-Seq Analysis

Two AF samples with mild degeneration (stage II–III, GSM7235336 and GSM7235334) and six AF samples with severe degeneration (stage III–IV) were analyzed. Eight clusters (resolution = 0.5) were identified (Figure 8A), and after finding marker genes, the AF cells were selected as high expression of ACAN and COL1A1 (cluster 2, 4, and 8) (Figure 8B and C). We then conducted subcluster analysis with a resolution of 0.4, and seven clusters were identified (Figure 8D). Monocle pseudotime trajectory showed that cluster 3 was in the early stage and cluster 5 was in the late stage (Figure 8E-G). The expression of CREB3 in cluster 3 was higher than that in cluster 5, but the expression of SCP2 and IDO1 was not significantly changed (Figure 8H).

Figure 8 Single-cell RNC-seq analysis of the annulus fibrosus (AF) samples. (A) Visualization of clustering by TSNE plot using 0.5 resolution. (B) Bubble plot showing the expression of the ACAN, SOX9, and COL1A1 related different clusters. (C) Six cell types were identified, including AF cell (AFC), fibroblast, macrophage, muscle stem cell (MSC), nucleus pulposus cell (NPC), and T cell. (D) Seven subclusters for the AF cells were identified using 0.4 resolution. (E-G) Monocle pseudotime trajectory showing the progression of seven subclusters of AF cells. (H) TSNE image of CREB3 in subcluster 3 and 5, showing the high expression in cluster 3.

Experimental Verification of MUC1 and CREB3 Expression

WB showed that the protein expression of GPX4 in degenerated human NP tissues and AF tissues was significantly downregulated, suggesting the existence of ferroptosis (Supplementary Figure S8). Based on the results of our previous study, we further verified the expression of MUC1 in degenerated NP. Immunohistochemical results showed that the protein expression of MUC1 in severely degenerated human NP tissues were significantly downregulated in contrast to mild NP samples (Figure 9A and B). And the demographic data of the participants is shown in Supplementary Table S2. An IVDD model of SD rats was established. The results of immunohistochemistry (Figure 9C and D) demonstrated the lower expression of MUC1 in the degenerated rat intervertebral disc tissues. In vitro, we used TBHP to treat human NP cells to establish an oxidative stress-related degenerated model. Immunofluorescence results showed that MUC1 was significantly downregulated in TBHP (50μM)-induced degenerated cells (Figure 9E).

Figure 9 Experimental validations of MUC1 and CREB3 in nucleus pulposus (NP) and annulus fibrosus (AF). (A and B) Immunohistochemical results demonstrated that the expression of MUC1 was decreased in degenerated NP tissues. (C and D) Immunohistochemical results showed that the expression of MUC1 decreased in degenerated rat NP tissues. (E) Immunofluorescence detection of protein expression of MUC1 in NP cells, scale bar=50 μm. (F and G) Immunohistochemical results of the expression of CREB3 in degenerated rat AF tissues. The Fe2+ (H and J, scale bar = 20μm) and ROS (H and K, scale bar = 100μm) intensity was increased in TBHP-treated AF cells, and mitochondrial morphology detection by transmission electron microscopy (TEM) (I, scale bar = 1000 nm, the red arrows indicate mitochondria). (L) Immunofluorescence detection of protein expression of CREB3 in AF cells, scale bar=50 μm. ***p < 0.001.

We detected the protein expression of CREB3 in AF tissues. The protein expression of CREB3 in degenerated rat AF tissues was significantly downregulated in contrast to controls (Figure 9F and G). TBHP was used to establish degenerated AF cell model. Previous study has shown that TBHP can induce oxidative stress and ferroptosis in AF cells.²⁰ We further verified that TBHP can significantly induce intracellular ferrous ion and ROS accumulation, and mitochondria morphology alteration (Figure 9H-K). Immunofluorescence results showed that CREB3 was significantly downregulated in TBHP (50μM)-induced degenerative AF cells (Figure 9L).

Discussion

The degeneration of NP or AF tissues could contribute to IVDD, a complex condition results in chronic back pain and disability.^1,2 Despite extensive efforts to explore the mechanisms underlying NP and AF degeneration, the outcomes of current treatment strategies for disc degeneration have been less than satisfactory.³³ Ferroptosis, a form of PCD, has emerged as a significant factor in the development and progression of various degenerative diseases.^10,11 Our previous literature review discussed the involvement of ferroptosis in IVDD and concluded that there was still considerable work to be done to uncover the underlying molecular mechanisms.³⁴ Bioinformatics analysis provides a valuable opportunity to identify hub ferroptosis-related genes in degenerative tissues and find out the novel molecular mechanisms of ferroptosis in disc degeneration.⁹ In this study, by integrated bioinformatics analysis and experimental verification, we aimed to identify and analyze the FRDEGs that associated with NP and AF degeneration.

Most studies on ferroptosis in intervertebral discs have predominantly focused on the NP, but research on ferroptosis in the AF remains limited. Jia et al³⁵ concluded that the selenium-SelK-GPX4 axis was responsible for the ferroptosis in mechanically overloaded-induced IVDD. Zhang et al³⁶ found that Cynarin could protect NP cells from ferroptosis and alleviate IVDD progression. Lu et al³⁷ elucidated that the m6A reader YTHDF1 could upregulate SLC7A11 expression, inhibiting ferroptosis in NP cells. Interestingly, we discovered that the DNA methyltransferase DNMT3B downregulated SLC40A1, promoting ferroptosis in NP cells.¹⁸ Xiang et al³⁸ revealed that the enhancement of Piezo1 channel-induced iron influx was implicated in the exacerbation of ferroptosis in NP cells. Zhu et al³⁹ demonstrated that the deubiquitinase USP11 stabilized SIRT3 to inhibit oxidative stress-induced ferroptosis in human NP cells. Regarding ferroptosis in the AF, Yang et al²⁰ reported that NCOA4-mediated ferritin selective autophagy is closely correlated with TBHP-induced ferroptosis in AF cells. Further research is necessary to investigate the molecular mechanisms of ferroptosis in AF degeneration and its impact on IVDD.

The identification of FRDEGs in degenerated NP and AF tissues constitutes a significant aspect of the current study. We identified a total of 15 FRDEGs in degenerated NP tissues and 18 FRDEGs in degenerated AF tissues. Functional enrichment analysis revealed that these FRDEGs were associated with important biological processes and pathways. In degenerated NP tissues, these genes were significantly enriched in oxidoreductase activity and the apoptosis pathway, which were contributors to NP degeneration.^7,40 Similarly, in degenerated AF tissues, the FRDEGs were enriched in response to oxidative stress, necroptosis, autophagy, apoptosis, and TNF signaling pathway, all of which play crucial roles in the degenerative process.^41–43 These genes exhibited dysregulated expression patterns, with some being upregulated and others downregulated, suggesting their potential involvement in the pathogenesis of IVDD. To further identify key FRDEGs contributing to NP and AF degeneration, we employed LASSO and SVM-RFE analysis. In NP degeneration, four genes were identified as feature FRDEGs. The area under the curve (AUC) of MUC1 was greater than 0.9, indicating that it has a certain level of accuracy and specificity in distinguishing degenerated NP tissues from normal samples. Similarly, in AF degeneration, five genes were identified as feature FRDEGs, with the AUC value of CREB3 greater than 0.9.

We examined the expression levels of the feature FRDEGs in the validation datasets and scRNA-seq data. The expression of the ferroptosis suppressor MUC1 was significantly downregulated in degenerated NP tissues or cells. Since the GSE147383 dataset included only two AF tissue samples in each group, the expression levels of SCP2, IDO1, and CREB3 did not show significant changes. The results of scRNA-seq showed that CREB3 was significantly downregulated in the late stage of AF cells. The results of IHC and immunofluorescence confirmed the low expression of MUC1 in degenerated NP tissues and CREB3 in degenerated AF tissues.

MUC1 is a macromolecular transmembrane protein belonging to the mucin family, and its involvement in ferroptosis has been studied in the context of cancer.^44–46 Liu et al⁴⁶ reported that apolipoprotein M holds a great importance for hepatocellular carcinoma treatment by downregulating MUC1 in cancer cells, contributing to the facilitation of ferroptosis. Liang et al⁴⁵ found that MUC1-related ferroptosis was regulated by lncRNA MALAT1/miR-145-5p axis in endometriosis. Treatment with Mesenchymal Stem Cells suppresses inflammatory responses in dextran sulfate sodium-induced colitis in mice by upregulating the expression of MUC1 and inhibiting ferroptosis.⁴⁷ However, the role of MUC1 low expression-induced ferroptosis in NP degeneration remains unknown. CREB3, a ferroptosis inhibitor, is a transcription factor localized to the endoplasmic reticulum⁴⁸ and is implicated in the survival of hepatitis B virus‐related hepatocellular carcinoma patients.⁴⁹ But its role in ferroptosis and AF degeneration is not fully studied.

GSVA enrichment analysis provided further insight into the biological functions of the hub FRDEGs. The GSVA analysis of MUC1 confirmed the activation of DNA cytosine deamination in the low-expression group of NP tissues. Cytosine deamination has been reported to cause mitochondrial DNA damage.⁵⁰ Cytidine deamination, as a form of DNA modification, plays a significant biological role in diseases.⁵¹ The GSVA results showed that autophagosome-lysosome fusion was activated in the low-expression group of CREB3 for AF tissues. Autophagosome-lysosome fusion represents the final step of autophagy, and abnormal activation of autophagy can promote cell death.⁵² These findings suggest that the dysregulation of these pathways may contribute to the degenerative processes in NP and AF tissues.

Moreover, immune infiltration analysis revealed a connection between the feature FRDEGs and the immune microenvironment in degenerated NP and AF tissues. MUC1 in NP showed no correlation with immune cells, indicating that MUC1-related ferroptosis may not contribute to IVDD through immune responses. Previous literature has found that MUC1 can regulate immune infiltrates in the clear cell renal cell carcinoma microenvironment by activating the complement system.⁵³ In the case of AF, CREB3 exhibited a negative correlation with Eosinophils, suggesting the involvement of this immune cell in AF degeneration. Eosinophils exert their functions by releasing various pro-inflammatory mediators.⁵⁴ This study is the first to reveal a correlation between CREB3 and Eosinophils, but its specific role in IVDD needs further investigation.

Furthermore, we conducted drug and ceRNA network analysis to explore potential therapeutic targets and regulatory mechanisms. In degenerated NP tissues, four drugs targeting MUC1 were selected, but most of them were antibodies used to target MUC1 in cancer treatment.^55–57 Their specific roles in NP degeneration require further investigation. However, in degenerated AF tissues, we were unable to identify underlying drugs for the hub gene CREB3, indicating the need for further research. The ceRNA networks in NP tissues revealed that miR-1291 and miR-4252 may regulate ferroptosis by targeting MUC1. In AF tissues, CREB3 was predicted to be regulated by six miRNAs, including miR-4270, miR-2355-5p, miR-542-3p, miR-34b-3p, miR-595, and miR-1225-5p. Focusing on the post-transcriptional level and epigenetic mechanisms in the miRNA-mRNA network provides a novel direction to discover new mechanisms of IVDD.⁵⁸ These findings provide a foundation for further exploration of the pathological mechanisms underlying IVDD and offer potential drug targets for intervention.

Our study has some limitations. Firstly, the mechanisms underlying the downregulation of MUC1 and CREB3 in degenerated NP and AF tissues, as well as their roles in ferroptosis, remain elusive. In vivo and in vitro experiments are necessary to elucidate these aspects. Secondly, the exact mechanisms of immune reactions induced by CREB3 downregulation in AF tissues need to be further investigated. Thirdly, although we established ceRNAs networks, their molecular mechanisms have not been comprehensively studied, which will be the focus of future research. Finally, the validation set for AF only contained two samples in each group, which may lead to inaccurate results.

Conclusion

In this study, we identified ferroptosis-related biomarkers for NP (MUC1) or AF (CREB3) degeneration through comprehensive bioinformatic and experimental verification. Then, we elucidated their important biological functions and relation to immunoinfiltration, and predicted their associated drugs or ceRNAs to define their mechanisms in the degeneration process. Our findings provide a foundation for further exploration of pathological mechanisms and novel therapeutic targets for IVDD.