Bioinformatics and Experimental Validation of FLVCR1 and SOX4 in Regulating Mitochondria-Macrophage Crosstalk in Disc Degeneration

Introduction

Intervertebral disc degeneration (IDD) is a prevalent degenerative musculoskeletal disorder and the main pathological cause of spinal degenerative diseases, such as intervertebral disc herniation and cervical spondylosis. IDD is characterized by excessive mechanical stress, apoptosis of nucleus pulposus (NP) cells, abnormal extracellular matrix (ECM) degradation, inflammatory responses, dysregulated autophagy, oxidative stress, and genetic susceptibility factors, and it is closely associated with symptoms such as low back pain.¹ It is particularly common among older adults and represents a major global public health problem due to its high prevalence and socioeconomic burden.² IDD is the principal pathological basis for degenerative spinal diseases such as disc herniation and cervical spondylosis, and is closely associated with chronic low back pain and functional disability.² Current treatment strategies, including conservative measures and surgical interventions, mainly aim to alleviate symptoms or correct structural abnormalities, but remain limited in halting disease progression or reversing biomechanical changes of the spine.^3,4 Therefore, it is critical to identify reliable biomarkers and elucidate molecular mechanisms to develop effective disease-modifying therapies.

Mitochondria, as essential organelles for energy metabolism and cellular homeostasis, play a pivotal role in the pathophysiology of IDD.⁵ Dysfunctional mitochondria impair respiratory chain activity, increase reactive oxygen species (ROS) accumulation, and induce apoptosis of NP cells.^5,6 Mitochondrial DNA mutations further compromise ATP synthesis and accelerate ECM degradation via MMP-3 activation, which correlates with the severity of disc degeneration.⁷ In parallel, macrophage polarization has emerged as another critical factor, with pro-inflammatory M1-like macrophages amplifying tissue damage, while M2-like macrophages contribute to tissue repair.⁸ Recent studies suggest that mitochondrial dysfunction may promote M1 polarization through multiple mechanisms, including mitochondrial ROS-mediated NF-κB activation, impaired mitophagy-induced cGAS-STING signaling via mtDNA leakage, and dysregulated mitochondrial dynamics such as DRP1-mediated fission, thereby sustaining a pro-inflammatory microenvironment.^8,9 Transcriptomic approaches offer a powerful tool to dissect the molecular basis of complex diseases by enabling systematic profiling of gene expression and identification of key regulators.¹⁰ Previous studies have successfully applied transcriptomics to reveal mechanisms by which sIL-13Rα2-Fc alleviates IDD in rat discs and to identify mediators of mannose’s protective effects on disc degeneration.^11,12

Previous studies have mainly focused on the independent roles of mitochondria-related genes (MRGs) and macrophage polarization-related genes (MPRGs) in the pathological progression of IDD. Whether MRGs and MPRGs exhibit direct interactions and share signaling pathways in disc tissue, and how they might synergistically drive IDD pathogenesis, has not yet been reported. Based on this, we propose for the first time a novel hypothesis that “MRGs and MPRGs may directly interact in disc tissue, share key signaling pathways, and jointly regulate the progression of IDD.” To test this hypothesis, datasets of IDD, MRGs, and MPRGs were obtained from public databases. A series of bioinformatics approaches were employed to identify key genes associated with mitochondria and macrophage polarization in IDD, and to further explore their potential biological functions, immunological characteristics, regulatory networks, and candidate therapeutic drugs. The expression of the identified key genes was further validated through animal experiments, with the overall workflow illustrated in Figure 1. This study addresses the current knowledge gap regarding the interplay between MRGs and MPRGs in the molecular mechanisms of IDD, providing new theoretical insights and experimental evidence for its diagnosis and treatment.

Figure 1 Workflow of the study. The study integrated differential expression analysis, WGCNA, and machine learning to identify candidate genes related to mitochondria and macrophage polarization in IDD. Key genes were validated in a rabbit IDD model and further analyzed for diagnostic value, biological functions, immune characteristics, regulatory networks, and potential therapeutic drugs.

Materials and Methods

Data Collection

The GSE70362 (16 IDD and 8 control NP samples, GPL17810) and GSE176205 (6 IDD and 3 control samples, GPL20301) datasets were obtained from the GEO database and used as training and validation sets, respectively.^1,13 These datasets were selected due to their relatively large sample sizes, clearly defined case-control groups, and comprehensive annotation. Both were generated on standardized platforms with reliable normalization, ensuring high data quality and reproducibility for subsequent analysis. In addition, 2030 MRGs were obtained from MitoCarta 3.0 database and gene set enrichment analysis (GSEA) database,¹⁴ while 35 MPRGs were retrieved from the molecular signatures database (MSigDB).¹⁵

Identification of Differentially Expressed Genes (DEGs) and Key Module Genes

The “limma” R package (v 3.56.0)¹⁶ was applied to analyze differential gene expression between IDD and normal (CK) samples (p < 0.05 and |log₂ Fold Change (FC)| > 0.5). The top 10 up-regulated and down-regulated genes, ranked by |log₂FC|, were visualized in volcano and heatmaps generated with the “ggplot2” (v 3.4.4)¹⁷ and “ComplexHeatmap” (v 2.14.0)¹⁸ R packages. The weighted gene co-expression network analysis (WGCNA) can investigate the associations between gene modules and traits, mine key genes, and facilitate deciphering biological mechanisms, disease pathologies, and identifying potential therapeutic targets. The Wilcoxon test was first used to assess the differences in MPRG expression between IDD and CK groups in the training set, and the significantly altered genes were designated as MPRGs1 for subsequent WGCNA. The “WGCNA” R package (v 1.71)¹⁹ was then employed to identify modules highly correlated with MPRGs1. Specifically, MPRGs1 scores were calculated across all training set samples using the GSVA algorithm, and differences between groups were tested by the Wilcoxon test with visualization in box plots via “ggpubr” (v 0.6.0).²⁰ Hierarchical clustering based on Euclidean distance of expression profiles was performed to detect and exclude outliers. A soft threshold (β) was determined by setting R² = 0.85, and modules were constructed with a minimum size of 100 genes and mergeCutHeight = 0.25, visualized through the “plotDendroAndColors” function. Finally, correlations between MPRGs1 scores and gene modules were computed (|correlation| > 0.3, p < 0.05), with results displayed via the “labeledHeatmap” function, and module genes significantly positively or negatively correlated with MPRGs1 scores were retained for further analysis.

Acquisition and Validation of Key Genes

The “ggvenn” R package (v 0.1.9)²¹ was used to extract the intersection of DEGs, module genes, and MRGs, which were designated as candidate genes. Enrichment analysis of these genes was conducted using the Kyoto Encyclopedia of Genes and Genomes (KEGG) and Gene Ontology (GO) with the “clusterProfiler” (v 4.7.1.003)²² and “org.Hs.eg.db” (v 3.16.0)²³ R packages (p < 0.05). To further refine candidate genes, the “glmnet” R package (v 4.1–4)²⁴ was employed to perform Least Absolute Shrinkage and Selection Operator (LASSO) regression analysis, and genes with the smallest lambda value and non-zero coefficient were retained. Meanwhile, the “caret” R package (v 6.0–93)²⁵ was used for Support Vector Machine-Recursive Feature Elimination (SVM-RFE), and genes with the highest classification accuracy were identified. The overlapping genes from the two algorithms were then extracted with the “ggvenn” R package and defined as characteristic genes.

Subsequently, expression differences of characteristic genes between IDD and CK groups in both the training and validation sets were compared using the Wilcoxon test (p < 0.05), and results were visualized with the “ggpubr” package. Characteristic genes that exhibited consistent expression trends with significant differences in both GSE70362 and GSE176205 datasets were identified as key genes. To further evaluate their diagnostic performance, receiver operating characteristic (ROC) curves were generated and the area under the curve (AUC) was calculated with the “pROC” R package (v 1.18.0),²⁶ with AUC values greater than 0.7 considered indicative of diagnostic potential.

Experiment Validation

New Zealand white rabbits weighing 1.8–2.2 kg were selected to establish a lumbar intervertebral disc herniation model using a closed puncture method.²⁷ The model group and control group each comprised 8 rabbits. Under anesthesia and disinfection, an 18-gauge puncture needle was inserted into the lumbar intervertebral disc under C-arm guidance in the model group, with a puncture depth of 5 mm. The control group underwent sham surgery using the same procedure without disc puncture. Postoperatively, animals were monitored daily for wound healing, activity, and feeding behavior. The inclusion criteria were: healthy rabbits that successfully underwent surgery and exhibited no abnormal conditions before modeling. Exclusion criteria included: development of postoperative infection, severe neurological deficits, or unexpected complications; perioperative death; poor general health that might interfere with experimental outcomes. Only animals that met the inclusion criteria and remained healthy throughout the study were analyzed.

Six weeks after modeling, intervertebral discs were harvested for gross observation. Successful modeling was confirmed by visible annulus fibrosus rupture and NP protrusion in the model group. NP tissues from both groups were collected for RT-qPCR to validate biomarker expression. Total RNA was extracted from the 16 samples using TRIzol reagent (Yeasen, Shanghai, China), reverse-transcribed into cDNA with the SureScript First-Strand cDNA Synthesis Kit (Yeasen, Shanghai, China), and amplified by quantitative PCR using 2× Universal Blue SYBR Green qPCR Master Mix (Yeasen, Shanghai, China). GAPDH was used as an internal control, and relative gene expression levels were calculated using the 2−ΔΔCt method. The animal experiments were performed at the Animal Experiment Research Center of Zhejiang Chinese Medical University (license number: SYXK-Zhe 2021–0012) and were approved by the Institutional Animal Care and Use Committee (IACUC, Approval No.: IACUC-20250310-18). All procedures followed the Guidelines for the Care and Use of Laboratory Animals (NIH, USA) and the Guidelines for the Ethical Review of Laboratory Animal Welfare of the Ministry of Science and Technology of China (2018 edition).

Construction and Assessment of Nomogram

A nomogram was created via the R package “rms” (v 6.5–0).²⁸ To appraise the performance of the nomogram prediction, the “rms” (v 6.5–0), “ggDCA” (v 1.2)²⁹ and “rmda” (v 1.6)³⁰ R packages were respectively utilized to produce the calibration curve, decision curve analysis (DCA), and clinical impact curve (CIC). Moreover, the Hosmer-Lemeshow test was carried out to assess the calibration capability of the model (p > 0.05).

GSEA

The GSEA was executed for key genes in the training set to probe the biological roles during the IDD growth. Subsequently, the gene set “c2.cp.kegg.v7.5.1.symbols.gmt” was derived from the MSigDB database as the reference gene set. Additionally, in the IDD and CK samples of the GSE70362 dataset, the “stats” (v 4.2.2)³¹ R package was applied to respectively compute the correlation coefficient between the key genes and all the genes, and the key genes were arranged according to correlation coefficients from largest to smallest. Then R package “clusterProfiler” (v 4.7.1.003) was utilized to perform GSEA, and the top 5 GSEA enrichment pathways for each key gene were visualized according to p from smallest to largest by the “enrichplot” R package (v 1.18.0)²³ (p < 0.05).

Immunoinfiltration Analysis

The “ssGSEA” algorithm of the “GSVA” R package (v 1.46.0)³² was applied to assess the 28 immune cells’ infiltration³³ in all samples of the training set, and the “pheatmap” package (v 1.0.12)³⁴ was used to draw a heatmap showing the immune cell score in the IDD and CK samples. The infiltration difference of 28 immune cells between 2 groups was contrasted by the Wilcoxon test, and a box plot was created by the “ggplot2” (v 3.4.4) (p < 0.05). Subsequently, the immune cells exhibiting significant differences were recorded as differential immune cells. The relations among differential immune cells and between differential immune cells and key genes in the GSE70362 dataset were explored through the R package “stats” (v 4.2.2), and correlation heatmaps were plotted using the packages “ggcorrplot” (v 0.1.4)³⁵ and “ggplot2” (v 3.4.4), respectively (|cor| > 0.3, p < 0.05).

Network Construction of Key Genes

The JASPAR database was applied to anticipate the transcription factor (TF) of key genes. Then, the miRNA of the key genes was forecasted through the miRTarBase database. Eventually, the TF-mRNA and miRNA-mRNA network were created by Cytoscape (v 3.9.1).

Drug Prediction and Molecular Docking

To delve into the interrelationships between key genes and potential drugs during the IDD development, the DSigDB was used to investigate the potential drugs associated with key genes, and the network of key genes and drugs was visualized by Cytoscape (v 3.9.1). In this study, the most significantly associated drug was used for molecular docking (p < 0.05). Subsequently, the relevant protein crystal structures were retrieved from the PDB database; the potential drug’s molecular structure was procured from the PubChem database. Then the cb - dock website was employed to carry out molecular docking.

Statistical Analysis

Bioinformatics analyses were performed using R language (v4.2.2). The differences between two groups were examined using the Wilcoxon rank-sum test, with statistical significance set at p < 0.05. Given the limited sample size, a 10-fold cross-validation strategy was applied during the construction of LASSO regression and SVM-RFE models to ensure efficient data utilization, model stability, improved generalizability, and reduced risk of overfitting. Furthermore, to control for type I error in multiple comparisons, the Benjamini-Hochberg false discovery rate (FDR) correction method was applied, and adjusted p-values were reported in the results.

Results

The 484 DEGs and 608 Key Module Genes Were Identified

A total of 484 differentially expressed genes (DEGs), including 199 up-regulated and 285 down-regulated genes, were identified, with the top 10 up-regulated (eg, IGFBP3, COL10A1, IGFBP1) and down-regulated genes (eg, IBSP, SCGB1D2, SCGB2A2) visualized in the volcano plot and heatmap (Figure 2A and B and Supplementary Table 1). To further explore genes related to macrophage polarization, Wilcoxon test analysis revealed that CRIPT, FAM241A, and SOAT1 were significantly differentially expressed between the IDD and CK groups, and the MPRGs1 score also showed a notable difference between the two groups. Notably, after Bonferroni correction, significant differences in MPRGs between the two groups were still observed (p = 0.032). This result indicated that the observed differences were not false-positive outcomes caused by multiple comparisons, further verifying the association between abnormal MPRG expression and the occurrence and progression of IDD (Figure 2C and Supplementary Table 2). Hierarchical clustering prior to WGCNA revealed no outlier samples, allowing all samples to be included in subsequent analysis. With β = 8, the scale-free topology fit index (R²) exceeded 0.85 and produced a co-expression network with good mean connectivity (Figure 2D). Finally, WGCNA identified 12 co-expression modules, among which the brown module (cor = 0.425) was significantly positively correlated and the pink module (cor = −0.452) was significantly negatively correlated with the MPRGs1 score, leading to the identification of 608 key module genes associated with IDD (Figure 2E).

Figure 2 Identification of Differentially Expressed Genes and Weighted Gene Co-expression Network. (A) Volcano plot of differentially expressed genes (DEGs). Red dots represent DEGs with increased expression levels in the IDD group, blue dots represent DEGs with decreased expression levels in the IDD group, and gray dots represent non-significant genes. (B) Heatmap of DEG expression levels. The upper heatmap shows the data distribution, where different colors indicate distribution densities: red areas represent regions with higher density, while blue areas represent regions with lower density. The middle part represents the grouping of samples, with blue indicating the control group and red indicating the disease group. The lower heatmap shows the gene expression levels in different samples, with blue indicating low expression levels and yellow indicating high expression levels. (C) Comparison of GSVA scores between groups. Blue represents the CK group; red represents the IDD group. The inter-group test was performed using the Wilcoxon test. (D) Selection of the soft-thresholding power for WGCNA. The left panel shows the scale-free fit index (y-axis) under different soft-thresholding powers (x-axis), with the red reference line representing R2 = 0.85. The right panel displays the network connectivity under different soft-thresholding powers. (E) Heatmap showing the correlation between gene modules and traits. The color blocks on the far left represent modules, and the color bar on the far right indicates the correlation range. In the heatmap in the middle, a darker color indicates a higher correlation, with red representing a positive correlation and blue representing a negative correlation; the numbers in each cell denote the correlation coefficient and statistical significance.

The Key Genes Were Identified

From the intersection of 484 DEGs, 608 key module genes, and 2030 MRGs, a total of 13 candidate genes were obtained, including CYP1B1, FLVCR1, and PLA2G4A (Figure 3A), which were mainly involved in cytochrome c release from mitochondria, tetrapyrrole binding, tryptophan metabolism, and ovarian steroidogenesis (Figure 3B and C; Supplementary Tables 3 and 4). Subsequently, LASSO regression analysis demonstrated that when log (lambda.min) = −2.6832, the model achieved the optimal fitting effect, yielding seven genes with non-zero coefficients (eg, CYP1B1, FLVCR1, and KYNU) (Figure 3D and E). Meanwhile, the SVM-RFE algorithm identified six genes with the highest model accuracy, including CYP1B1, FLVCR1, and KYNU (Figure 3F). By intersecting the two algorithms, four characteristic genes were screened, namely CYP1B1, FLVCR1, KYNU, and SOX4 (Figure 3G). Further Wilcoxon test analysis revealed that FLVCR1 and SOX4 exhibited consistent expression trends in both the GSE70362 and GSE176205 datasets and showed significant differences between IDD and CK samples (p < 0.05) (Figure 3H and I). Importantly, ROC curve analysis indicated that FLVCR1 (AUC = 0.875) and SOX4 (AUC = 0.844) displayed good diagnostic performance (AUC > 0.7) (Supplementary Figure 1), and were therefore selected as key genes for subsequent analyses.

Figure 3 Functional Enrichment and Gene Screening. (A) Venn diagram showing the overlap among DEGs, WGCNA module genes, and mitochondria-related genes (MRGs). The numbers represent the number of genes, and the percentages indicate the proportion of genes. (B) KEGG enrichment analysis results. The x-axis represents log10 (p-value), and the y-axis denotes different KEGG pathways. The size of the circles indicates the number of genes enriched in the corresponding pathway, with larger circles representing a greater number of enriched genes. (C) GO enrichment analysis results. Different colors represent different GO categories, and the longer the bar, the greater the number of genes enriched in the corresponding GO term. (D) Cross-validation error curve for feature gene selection using the LASSO regression algorithm. The plot illustrates the cross-validation error of the model as a function of the regularization parameter lambda, which is used to determine the optimal lambda value (lambda.min) that minimizes the model bias. (E) LASSO logistic regression coefficient penalty plot. The plot shows the coefficients of different feature genes as a function of log(lambda.min), where log(lambda.min) denotes the log-transformed value of lambda.min when the bias of the LASSO model is minimized. Under this log(lambda.min) value, the model achieves relatively high fitting performance. (F) Feature gene screening using the SVM-RFE algorithm. The x-axis represents the number of genes; the y-axis denotes the accuracy rate after cross-validation. The solid blue dot indicates the point at which the model achieves the highest accuracy rate. (G) Venn diagram showing the intersection of genes identified by LASSO and SVM-RFE. (H) Expression levels of key genes in training (GSE70362) dataset. In the box plots, blue represents the CK group and red represents the IDD group. Asterisks (*) indicate inter-group differences in the expression levels of characteristic genes after the Wilcoxon test. *** p < 0.001 and ** p < 0.01. (I) Expression levels of key genes in validation (GSE176205) dataset. In the box plots, blue represents the CK group and red represents the IDD group. Asterisks (*) indicate inter-group differences in the expression levels of characteristic genes after the Wilcoxon test. * p < 0.05 and ns (non-significant) p > 0.05.

Experiment Validation

The procedure of C-arm-guided closed puncture of rabbit lumbar intervertebral discs was successfully completed (Figure 4A). Six weeks post-operation, dissection revealed significant annulus fibrosus rupture, NP protrusion, and disc degeneration in the model group, while the control group maintained intact annulus fibrosus and NP (Figure 4B and C). RT-qPCR analysis of biomarker expression between groups demonstrated that the expression trends of FLVCR1 and SOX4 in the model group were consistent with our previous analytical results, showing statistically significant differences compared to the control group (Figure 4D and E).

Figure 4 Validation in Intervertebral Disc Herniation Model. (A) Establishment of rabbit lumbar disc herniation model using C-arm fluoroscopy-guided puncture. The black arrow indicates the position of the puncture needle. (B) Nucleus pulposus tissue in control group. (C) Nucleus pulposus tissue in model group. (D) Differences in the relative mRNA expression levels of FLVCR1 gene between the CK group (8 samples) and the IDD group (8 samples). The inter-group statistical test was performed using Student’s t-test, *** p < 0.001. (E) Differences in the relative mRNA expression levels of SOX4 gene between the CK group (8 samples) and the IDD group (8 samples). The inter-group statistical test was performed using Student’s t-test, *** p < 0.001.

The Good Nomogram Was Constructed and Evaluated

A nomogram was constructed, which could clearly show the relationship between the total points of all variables and the risk of IDD development; the higher the total points, the greater the risk of developing IDD. Through the calibration curve, DCA and CIC, the research evaluated the nomogram constructed based on FLVCR1 and SOX4. The evaluation demonstrated that the nomogram could predict the probability of disease occurrence well, and according to the Hosmer-Lemeshow test, it had excellent calibration efficiency (p = 0.617) (Figure 5A and D).

Figure 5 Construction and Evaluation of Nomogram Model. (A) Nomogram based on key genes. Points (Scoring Axis): The top scale represents “scoring values” ranging from 0 to 100, which is used to convert gene expression values into corresponding scores. FLVCR1 and SOX4: The corresponding scales indicate the value ranges of the expression levels of these two genes. The corresponding position can be found on the axis based on the actual expression value of FLVCR1, and the score can be obtained by vertically mapping to the “Points” axis. Total Points (Total Scoring Axis): Add the scores of FLVCR1 and SOX4 together; the resulting total score can be located on this axis (ranging from 0 to 160). Risk (Risk Axis): Vertically map the total score to this axis (ranging from 0.1 to 0.9), and the corresponding intervertebral disc-related risk probability can be obtained. (B) Calibration curve of the nomogram. The values on the X-axis represent the probability of the event (eg, disease) predicted by the model. The values on the Y-axis represent the actual probability of the event occurring. Hosmer-Lemeshow p = 0.617; since the p-value is greater than 0.05, it indicates that the nomogram has a good degree of fit between the predicted probability and the actual proportion of disease occurrence, with no significant difference. Apparent (solid blue line): The calibration curve of the nomogram in the current dataset (without bias correction), which reflects the correspondence between the predicted probability and the actual disease proportion. Bias-corrected (solid black line): The calibration curve after bias correction, which more accurately reflects the predictive performance of the nomogram in external or validation populations. Ideal (dashed black line): The ideal calibration curve, representing the scenario where the predicted probability is completely consistent with the actual disease proportion (ie, the predicted values and actual values completely overlap on the 45° line). B = 1000 repetitions, boot: The Bootstrap method with 1000 repetitions was used for calibration validation. Mean absolute error = 0.055: The mean absolute error is 0.055, which reflects the average deviation between the predicted values and actual values; a smaller value indicates better calibration performance. n = 24: The sample size is 24 cases. (C) Decision curve analysis (DCA) of the nomogram. X-axis (High Risk Threshold): The high-risk threshold, ie, the critical value of disease risk at which intervention is deemed necessary. Y-axis (Net Benefit): Net benefit, which reflects the additional benefit obtained by using the predictive indicator for decision-making (under a specific high-risk threshold) compared to “no prediction at all (None)” or “intervening in all study subjects (All)”. FLVCR1 (Orange line): Net benefit curve when prediction is based solely on the expression level of the FLVCR1 gene. SOX4 (blue line): Net benefit curve when prediction is based solely on the expression level of the SOX4 gene. Nomogram (black line): Net benefit curve when prediction is based on the nomogram (integrating factors such as FLVCR1 and SOX4). All (gray line): Net benefit curve assuming all study subjects are classified as high-risk and receive intervention (usually used as a reference to reflect the scenario of “over-intervention”). None (gray line, close to the X-axis): Net benefit curve assuming no prediction or intervention is performed (reflecting the baseline scenario of “no intervention at all”). (D) Clinical impact curve of the nomogram. Lower X-axis (Cost:Benefit Ratio): Represents the ratio of cost to benefit, ranging from (1:100) to (100:1), which reflects the scenarios under different proportions of cost input and benefit return. Upper X-axis (High Risk Threshold): Denotes the high-risk threshold, ranging from (0.0) to (1.0), which serves as the criterion for defining individuals classified as high-risk. Y-axis (Number high risk, out of 1000): Displays the number of high-risk individuals among every 1000 study subjects. Solid red line (Number high risk): Shows the variation in the total number of individuals classified as high-risk among every 1000 study subjects under different high-risk thresholds and cost-benefit ratios. Dashed blue line (Number high risk with event): Represents the variation in the number of individuals who are both classified as high-risk and have actually experienced the target event (eg, disease onset) among every 1000 study subjects under different high-risk thresholds and cost-benefit ratios. The solid red line and dashed blue line correspond to the two statistical metrics: “number of high-risk individuals” and “number of high-risk individuals with event”, respectively.

The 27 Pathways of FLVCR1 and 55 Pathways of SOX4 Were Found

FLVCR1 and SOX4 were significantly enriched in 27 and 55 pathways, respectively, and the study showed a score map illustrating the top 5 enriched pathways for the 2 key genes, among which “systemic lupus erythematosus” and “pathways in cancer” were the top pathway in the GSEA results of FLVCR1 and SOX4, respectively (p < 0.05). Most importantly, among the top 5 pathways of FLVCR1 and SOX4, “chronic myeloid leukemia” was a co-enrichment pathway (Figure 6A and B and Supplementary Tables 5 and 6). These results facilitate our comprehension of the biological pathways involved in the progression of IDD and lay a foundation for further investigation into the pathogenesis of IDD in the future.

Figure 6 GSEA Analysis of Key genes. (A) GSEA analysis results for FLVCR1. The top section shows the names of the enriched KEGG pathways. The upper part illustrates the calculation process of the enrichment score (ES): an ES value is calculated each time a gene is encountered from left to right, and these values are connected to form a line. A particularly prominent peak on the far left/right represents the ES value of the gene set in terms of phenotype. The middle part: each line represents a single gene in the gene set and its ranked position in the gene list.The lower part is the distribution map of the rank values of all genes; the ordinate is the “ranked list”, which refers to the ranking based on Spearman correlation. (B) GSEA analysis results for SOX4.

The Relationship Between the 6 Differential Immune Cells and the Key Genes Were Identified

In addition to the above findings, immunoinfiltration analysis was performed. The study found that the immune scores of central memory CD4 T cells and plasmacytoid dendritic cells were higher, while those of neutrophils and activated B cells were lower. It was found that 6 immune cells were significantly different between the IDD and CK samples by Wilcoxon test, including activated dendritic cells, Monocyte, etc. (p < 0.05) (Figure 7A and B). Subsequently, it was also displayed that the correlation between type 1 T helper cells and myeloid-derived suppressor cells (MDSC) was the strongest (cor = 0.75 and p < 0.05) (Figure 7C). Moreover, except for type 1 T helper cells, FLVCR1 was negatively correlated with all the other 5 differential immune cells, while SOX4 exhibited a significant positive correlation with 6 differential immune cells (p < 0.05). More importantly, monocytes had the most significant correlation with FLVCR1 (cor = −0.76 and p < 0.05), which provided new directions for understanding the pathogenesis (Figure 7D).

Figure 7 Immune Cell Infiltration Analysis. (A) Immune cell score comparison between healthy controls and disease samples. Blue represents the CK group, and red represents the IDD group. The color of the heatmap reflects the immune score: a higher immune score corresponds to a color closer to red, while a lower immune score corresponds to a color closer to blue. (B) Differential immune cell infiltration between healthy and disease groups. Blue represents the CK group, and red represents the IDD group. The inter-group statistical test was conducted using the Wilcoxon test, ** p<0.01, * p<0.05. (C) Correlation heatmap of differentially infiltrated immune cells. Red represents a positive correlation, and blue represents a negative correlation; a darker color indicates a higher correlation. Asterisks (*) reflect the significance level, while the numbers represent the correlation coefficients. Specifically, *** p<0.001, ** p<0.01, * p<0.05. (D) Correlation between key genes and immune cells. The x-axis represents key genes, and the y-axis represents differential immune cells. Large dots indicate a p-value < 0.05, while small dots indicate a p-value > 0.05. The color of the dots represents the correlation: blue denotes a negative correlation, red denotes a positive correlation, and gray denotes no correlation.

Network of Key Genes Was Conducted

The 5 TFs of FLVCR1, and the 3 TFs and 17 miRNAs of SOX4 were predicted from the JASPAR and miRTarBase databases, respectively, while the miRNAs of FLVCR1 were not predicted. Then, the TF-mRNA and miRNA-mRNA networks were constructed according to the TFs, miRNAs and mRNAs by using Cytoscape. Figure 8A and B showed that SOX4 had regulatory relationships with many TFs and miRNAs, such as YY1- SOX4, hsa-miR-335-5p - SOX4; however, only 3 TFs related to FLVCR1 were obtained in this study. It was worth noting that FOXC1 was the common TF of FLVCR1 and SOX4, which suggested that FOXC1 may influence the occurrence of IDD by regulating FLVCR1 and SOX4.

Figure 8 Regulatory Network of Key Genes. (A) Key gene–transcription factor (TF) regulatory network. Red nodes represent key genes, and yellow nodes represent TFs. (B) Key gene–miRNA regulatory network. Red nodes represent key genes, and yellow nodes represent predicted miRNAs.

The Relationship Between Key Genes and Drugs Was Investigated

The drugs associated with FLVCR1 and SOX4 were predicted from the DSigDB database. The analysis indicated that 23 and 48 drugs were associated with FLVCR1 and SOX4, respectively (Figure 9A). It was found that the binding energy between SOX4 and MCDF was −6.2, while the binding energy between FLVCR1 and Remifentanil was −6.6 kcal/mol. The binding relationship between FLVCR1 and Remifentanil was closer, which may furnish a reference for the IDD management in the future (Figure 9B and C and Supplementary Table 7).

Figure 9 Drug Prediction for Key Genes. (A) Predicted drugs targeting FLVCR1 and SOX4 from the DSigDB database. Red nodes represent key genes, and yellow nodes represent predicted potential drugs. (B) Interaction modes between SOX4 and candidate drugs. The left panel shows the overall view of the binding between the protein (with its overall structural outline outlined by yellow lines) and ligand molecules (colored small molecular structures). The right panel is a locally enlarged view of the binding region, where the yellow helical structures represent the secondary structure of the protein encoded by the key gene, and the small red and cyan molecules correspond to the structures of different drugs. TRP-69 denotes the amino acid residue involved in the interaction between the protein and the drugs. (C) Interaction modes between FLVCR1 and candidate drugs.

Discussion

IDD is driven by mitochondrial dysfunction, macrophage polarization, and chronic inflammation, yet its molecular mechanisms remain poorly understood.³⁶ In this study, transcriptomic analysis identified FLVCR1 and SOX4 as key genes linked to mitochondrial homeostasis and immune regulation in IDD, with subsequent validation confirming their central role in disease progression. Although diagnostic modeling demonstrated relatively high performance (AUC>0.8) and the nomogram incorporating these genes showed good predictive accuracy, these results should be interpreted cautiously due to the limited sample size and indirect immune infiltration estimation by ssGSEA. Overall, FLVCR1 and SOX4 emerge as potential biomarkers bridging mitochondrial stress and immune dysregulation in IDD, while additional studies are needed to validate their clinical applicability.

FLVCR1 encodes a mitochondrial heme exporter essential for iron balance and phospholipid biosynthesis through the Kennedy pathway, playing pleiotropic roles in development and degeneration.³⁷ Its deficiency disrupts choline metabolism and mitochondrial function, causing embryonic lethality, while variants contribute to neurodevelopmental disorders.³⁸ In skeletal biology, FLVCR1 cooperates with TCOF1 to maintain mitochondrial integrity and delay bone aging.³⁹ Downregulation in IDD correlates with mitochondrial dysfunction, ferroptosis, and NP apoptosis, and may regulate macrophage polarization through ROS signaling.⁴⁰ In our study, FLVCR1 expression was significantly reduced in IDD tissues, supporting previous reports that its downregulation may disrupt mitochondrial metabolic homeostasis, exacerbate ferroptosis and apoptosis in NP cells, and indirectly alter ROS-mediated macrophage polarization, collectively driving IDD onset and progression. By contrast, SOX4, an HMGB-family TF, drives degeneration via NF-κB/TGF-β signaling, MMP activation, and chondrocyte senescence.⁴¹ In IDD, it promotes ECM breakdown and apoptosis through miR-335-5p regulation, while suppression of SOX4 by miR-499a-5p or PRP-derived EV-delivered miR-25-3p has been shown to attenuate NP cell degeneration.^41–44 Moreover, SOX4 was found to be upregulated in puncture-induced IDD mouse models, and icariin alleviated low back pain by inhibiting SOX4 via the PI3K/Akt pathway.⁴⁵ Consistently, our analyses and rabbit model validation revealed significant SOX4 upregulation in IDD, suggesting that SOX4 may act as a central regulatory factor by promoting NP cell dysfunction, ECM degradation, and apoptosis, as well as indirectly modulating immune microenvironmental changes such as Th17 cell polarization. Together, these findings highlight FLVCR1 downregulation and SOX4 upregulation as complementary molecular events that synergistically disrupt mitochondrial homeostasis, immune regulation, and tissue integrity, thereby driving the pathological progression of IDD. Furthermore, regulatory network analysis predicted FOXC1 as an upstream TF of both genes, raising the possibility that it acts as a “master switch” in mitochondria–immune crosstalk during IDD, a novel hypothesis requiring further experimental validation.

Immune infiltration analysis revealed six immune cell types with significant differences between IDD and control samples, including activated dendritic cells, CD56dim natural killer cells, central memory CD4 T cells, MDSCs, monocytes, and type 1 T helper cells. Among them, FLVCR1 exhibited the strongest negative correlation with monocytes (cor = −0.76), suggesting that it may suppress monocyte infiltration into degenerative disc tissue, thereby reducing the release of pro-inflammatory cytokines such as TNF-α and IL-6 following monocyte activation.⁴⁶ Combined with its known roles in maintaining mitochondrial homeostasis and inhibiting ferroptosis in NP cells,⁴⁰ FLVCR1 may collectively attenuate local inflammatory responses and exert anti-inflammatory protective effects. In contrast, SOX4 showed significant positive correlations with differential immune cells, particularly monocytes, central memory CD4 T cells, and MDSCs, indicating that it may exacerbate disc degeneration by recruiting pro-inflammatory immune cells and activating NF-κB/TGF-β signaling pathways.⁴⁷ This could enhance monocyte differentiation into M1 macrophages and amplify the inflammatory potential of central memory CD4⁺ T cells, thereby promoting ECM degradation and disrupting the immune-inflammatory microenvironment. Monocytes, as key components of innate immunity, include classical, intermediate, and non-classical subsets,⁴⁸ with distinct chemokine responses.⁴⁹ Their infiltration has been reported as both increased and decreased in IDD,^2,50 highlighting contradictory observations that may be explained by the opposing effects of FLVCR1 and SOX4 on monocyte regulation. Elevated FLVCR1 expression appears to inhibit monocyte infiltration and attenuate inflammation, whereas predominant SOX4 expression enhances monocyte infiltration and accelerates degeneration, forming a “balancing regulatory axis” in the IDD inflammatory microenvironment. Dysregulation of this axis may drive the transition from homeostatic maintenance to inflammatory progression. Given the complexity and heterogeneity of immune infiltration, with macrophages broadly polarized into pro-inflammatory M1 and anti-inflammatory M2 phenotypes and monocytes exhibiting similar heterogeneity, single-cell RNA sequencing has provided evidence of increased monocyte/macrophage populations in degenerated discs and revealed their extensive interactions with CD82⁺ progenitor NP cells via MIF and NF-κB signaling.^50–53 These findings underscore the importance of scRNA-seq for elucidating the cell-specific expression of key genes and establishing a precise framework for targeted immunotherapy in IDD.

Drug prediction identified 23 potential compounds targeting FLVCR1 (eg, Remifentanil) and 48 compounds targeting SOX4 (eg, MCDF). Molecular docking demonstrated strong binding affinity between FLVCR1 and Remifentanil, and previous studies have confirmed that abnormal FLVCR1 expression leads to mitochondrial dysfunction in disc cells and increased release of inflammatory cytokines such as IL-1β.⁴⁰ This suggests that Remifentanil may exert protective effects by locally binding to FLVCR1, thereby improving mitochondrial function within the disc and suppressing inflammatory pathways, while its role as an opioid analgesic may also influence systemic immune responses, given that opioids have been reported to modulate peripheral immune cells, which could further impact the local disc immune microenvironment.⁵⁴ Meanwhile, MCDF, an aryl hydrocarbon receptor (AhR) agonist, may attenuate NF-κB-driven inflammation by activating antioxidant signaling pathways, consistent with prior evidence linking opioid receptors to mitochondrial regulation and AhR activation to oxidative stress mitigation, thereby providing novel therapeutic avenues for the management of IDD.^2,55

This study highlights FLVCR1 and SOX4 as central regulators of mitochondrial-immune crosstalk in IDD and identifies Remifentanil and MCDF as potential therapeutic candidates. Nonetheless, several limitations should be acknowledged. First, immune infiltration was inferred using ssGSEA rather than direct experimental validation, which may not fully capture cell-level features or actual infiltration patterns. Future studies will employ flow cytometry to sort and quantify specific immune cell subsets, and use immunohistochemistry/immunofluorescence to examine the localization and abundance of immune cells within disc tissues, thereby elucidating the regulatory mechanisms of the IDD immune microenvironment. Second, the relatively small dataset size may have limited statistical power and increased susceptibility to sample heterogeneity. To overcome this, multicenter collaborations with larger cohorts and clinical follow-up data will be pursued to validate gene expression and strengthen the translational value of these findings.