In SCI experimental models, urinary dysfunction is a common complication in animals. Female rats are typically selected for SCI studies owing to their relatively short and straight urethra, which simplifies postoperative management and minimizes the likelihood of urinary tract infections. For this study, 18 female SPF-grade SD rats, aged 8 weeks and weighing approximately 220 ± 20 grams, were sourced from Zhuhai Baishitong Biotechnology Co., Ltd. The animals were housed under standardized conditions, including controlled temperature, humidity, and a consistent light–dark cycle, with unrestricted access to food and water.

The rats were randomly assigned to three experimental groups: a sham surgery group (n = 6), a 2-day post-SCI group (n = 6), and a 7-day post-SCI group (n = 6). A SCI impact model was established to replicate human SCI. All procedures were performed at the Experimental Animal Research Center of Zhujiang Hospital, Southern Medical University (license number: SYXK (Guangdong) 2023-0215). The study received approval from the Ethics Committee of Zhujiang Hospital, Southern Medical University (approval number: LAEC-2024-120; approval date: June 21, 2024). Establishing the SCI model is a fundamental step in neuroscience research. In this study, a precise spinal cord impact device was employed to induce SCI in rats. Following anesthesia with isoflurane to ensure pain relief and unconsciousness, the rats were positioned prone, and the dorsal hair was shaved to prepare the area for surgery.

The surgical procedure involved a midline incision approximately 1.5 cm in length on the back. Using bony landmarks for guidance, the paraspinal muscles were dissected bluntly to fully expose the T10 vertebral lamina. A laminectomy was subsequently performed at the T10 level to facilitate the creation of the SCI model.

The study employed the 68099II spinal cord precision impactor from Shenzhen Ruiwoerd Life Science and Technology Co., Ltd., which offers precise control over impact speed, depth, and dwell time, ensuring accurate SCI induction. The impactor settings included a size 2 impact head, an impact velocity of 2 m/s, a depth of 1.2 mm, and a dwell time of 0.5 s. This controlled impact generated a significant hematoma at the targeted spinal segment.

Following SCI, the hematoma underwent dynamic changes over time. During the hyperacute phase, bleeding and inflammation exacerbated the hematoma, which gradually resolved in the later stages. This progression was typical of the natural course after SCI. The contusion also induced temporary involuntary spasms in the hind limbs and tail rigidity, serving as indicators of a successfully established and standardized SCI model.

Postoperatively, the muscles and skin were sutured in layers. To support the animals’ physiological function, manual bladder expression was performed three times daily for 1 week following the injury to ensure adequate urine elimination. In the sham group, rats underwent a laminectomy at the T10 level, preserving the spinal cord and maintaining normal BBB (Basso, Beattie, and Bresnahan) locomotor scores postoperatively.

At 2 and 7 days post-SCI, animals were anesthetized with a 40 mg/kg intraperitoneal injection of 1% pentobarbital sodium. Following cardiac perfusion with phosphate-buffered saline, spinal cord samples were harvested for subsequent qPCR analysis.

The gene expression omnibus (GEO) database¹ provided SCI-associated microarray expression profiles, specifically GSE109902 and GSE166009. The GSE109902 dataset, derived from the GPL22396 platform, included 20 spinal cord tissue samples from SCI rats and 10 from control rats. The GSE166009 dataset, based on the GPL27597 platform, contained 12 spinal cord tissue samples, consisting of nine from SCI rats and three from controls. Additionally, 1,136 mitochondria-related human genes were obtained from the MitoCarta 3.0 database² and subsequently mapped to rat homologs, resulting in 1,119 rat homologous genes (MRGs) (Supplementary Table S1). A total of 431 ferroptosis-related human genes retrieved from the literature were similarly converted to their rat homologs, yielding 437 rat genes (FRGs) (Supplementary Table S2) (12).

Expression data of MRGs and FRGs from the GSE109902 dataset were used to calculate single-sample gene set enrichment analysis (ssGSEA) scores via the GSVA package (v 1.42.0) (13). The Wilcoxon test was then applied to compare MRG and FRG scores between SCI and control groups (p < 0.05). A co-expression network was constructed using MRG and FRG scores as traits within the GSE109902 dataset through the weighted gene co-expression network analysis (WGCNA) package (v 1.7.1) (14). Outlier samples were first removed using hierarchical clustering based on Euclidean distance. The scale-free network evaluation coefficient (R²) and mean connectivity were subsequently used to determine the optimal soft-threshold power (β). Genes exhibiting similar expression patterns were grouped into the same module using the hybrid dynamic tree cutting algorithm, with parameters set as minModuleSize = 20, deepSplit = 2, and mergeCutHeight = 0.15. Different colors were assigned to modules for visualization. A heatmap of trait-module correlations was generated, emphasizing modules significantly associated with both MRG and FRG scores (|cor| >0.3 and p < 0.05). Genes within these modules were identified as key module genes, strongly correlated with MRG and FRG scores.

Differentially expressed genes (DEGs) between SCI and control samples in the GSE109902 dataset were identified using the DESeq2 package (1.34.0) (15), with thresholds set at |log2 Fold Change (FC)| >0.5 and p < 0.05. A comprehensive visualization of DEG distribution was achieved by generating a volcano plot and heat map through the ggplot2 package (v 3.4.1) (16) and ComplexHeatmap-package (v 2.14.0) (17). A Venn diagram, constructed with the VennDiagram package (v 1.7.1) (18) was used to identify overlapping genes between key modular genes and DEGs, which were then classified as candidate genes.

These candidate genes were subjected to Gene Ontology (GO) annotation and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis via the clusterProfiler package (v 4.2.2) (19), with pathways exhibiting p < 0.05 considered statistically significant. The co-expression network of candidate genes was constructed by mapping the genes onto GeneMANIA.³

Two machine learning algorithms, least absolute shrinkage and selection operator (LASSO) and Boruta, were applied to identify feature genes associated with SCI from the candidate genes in the GSE109902 dataset. For LASSO analysis, the glmnet package (v 4.1–4) (20) was employed, with a lambda value of 0 selected as optimal for feature gene selection. Boruta analysis was performed using the Boruta package (v 8.0.0) (21) with default parameters. Subsequently, the common feature genes identified by both algorithms were determined through overlapping the results using the VennDiagram package (v 1.7.1).

Expression trends of the common feature genes were further compared between SCI and control samples in the GSE109902 and GSE166009 datasets using the Wilcoxon test. Emphasis was placed on genes exhibiting stable expression, defined as those with consistent expression patterns across both datasets and significant differential expression between groups (p < 0.05), which were selected as biomarkers for this study.

To assess tissue specificity, gene expression levels were analyzed using the online tool genotype-tissue expression (GTEx, https://gtexportal.org/home/).

An artificial neural network (ANN) diagnostic model was developed to predict SCI likelihood based on the expression levels of identified biomarkers. After normalizing data from the GSE109902 dataset, the model was constructed using the neuralnet package (v 1.44.2) (22), employing the min-max method and configuring three hidden layers. To evaluate the model’s predictive accuracy, the receiver operating characteristic (ROC) curve was generated using the pROC package (v 1.18.0) (23). An area under the curve (AUC) exceeding 0.7 was considered indicative of high accuracy.

To elucidate the signaling pathways related to the biomarkers, GSEA was conducted on the GSE109902 dataset. Initially, Spearman correlation analysis was performed between the biomarkers and other genes across all samples using the psych package (v 2.1.6). Genes were then ranked based on their correlation coefficients, from highest to lowest. Concurrently, GSEA was carried out using the clusterProfiler package (v 4.2.2), with the “c2.cp.kegg.v7.4.symbols.gmt” gene set from the Molecular Signatures Database (MSigDB, http://www.gsea-msigdb.org/gsea/msigdb/index.jsp) as the reference. The thresholds for statistical significance were set at adj. p < 0.05 and |Normalized Enrichment Score (NES)| >1.

Additionally, the biomarkers were analyzed using the GeneMANIA database (see text footnote 3) to predict genes associated with the biomarkers’ functions and the biological processes in which they participate.

Immune cell infiltration plays a critical role in characterizing disease defense mechanisms, offering insights into the extent of immune involvement in specific pathologies. In this study, infiltration levels of 22 immune cell types were quantified in samples from the GSE109902 dataset using the CIBERSORT algorithm (24). Differential immune cell infiltration between SCI and control samples was then assessed via the Wilcoxon test (p < 0.05). Additionally, Spearman correlation coefficients were calculated to examine the relationships between the immune cell types, as well as between biomarkers and immune cell infiltration, based on the expression matrix of biomarkers and immune cell levels.

The SCI process is typically associated with the release of inflammatory factors (25). To explore the relationship between biomarkers and inflammatory factors, 200 inflammation-related genes (IRGs) were retrieved from the MSigDB database using the search term “HALLMARK_INFLAMMATORY_RESPONSE.” These 200 IRGs formed the background gene set for subsequent analysis. ssGSEA was then conducted on the GSE109902 dataset via the GSVA package to calculate the IRG score for each sample. Differences in IRG scores between the SCI and control groups were assessed (p < 0.05). Additionally, the relationship between biomarkers and IRGs was examined by calculating Spearman correlation coefficients and their significance, with the thresholds set at |cor| >0.3 and p < 0.05.

To investigate the regulatory factors targeting the biomarkers, upstream transcription factors (TFs) were identified through the miRNet database.⁴ Additionally, microRNAs (miRNAs) regulating the biomarkers were detected using miRecords,⁵ miRTarBase,⁶ and TarBase.⁷ The upstream long non-coding RNAs (lncRNAs) of the identified miRNAs were further explored in the same databases to elucidate the molecular mechanisms governing the biomarkers. Finally, the TF-biomarker and lncRNA-miRNA-mRNA (biomarker) networks were visualized using Cytoscape software (v 3.10.1) (26) based on the predicted results.

Potential drugs targeting biomarkers were identified through an in-depth analysis of the Enrichr database⁸ to support drug prediction and inform therapeutic strategies for SCI. The relationship between biomarkers and drugs was further examined by obtaining biomarker protein structures from the Protein Data Bank (PDB) database.⁹ Molecular docking was conducted using the Chemical and Biological Docking (CB-Dock) platform,¹⁰ focusing on interactions with binding energies below 5 kcal/mol, indicative of stronger biomarker-drug affinity. The docking results were subsequently visualized using PyMOL software (v 3.0.3) (27) to facilitate structural interpretation and analysis.

Spinal cord samples from T10 were collected at 2 and 7 days post-SCI, immediately frozen, and stored at −80°C (six samples per group). Total RNA was extracted using TRIzol reagent (Thermo Fisher Scientific, Inc.) and reverse transcribed into complementary DNA (cDNA). Thermal cycling was performed under the following conditions: denaturation at 95°C for 30 s, annealing at 95°C for 3 s, and extension at 60°C for 30 s. GAPDH mRNA served as the endogenous control, and results were quantified using the 2−ΔΔCT method. The primers utilized in this study were listed in Table 1. Since the experiment involved comparing data from postoperative day 2 and postoperative day 7 with the control group, we selected one-way analysis of variance (one-way ANOVA) as the statistical analysis method.

Data processing and statistical analysis were conducted using R software (v 4.2.2). Group differences were assessed using the Wilcoxon test, with statistical significance set at p < 0.05.

The expression matrices of MRGs and FRGs in the GSE109902 dataset revealed significant differences in both the MRGs and FRGs scores between the SCI and control groups, with the SCI group exhibiting a lower MRGs score and a higher FRGs score (p < 0.05) (Figures 1A,B), as determined by the ssGSEA algorithm. WGCNA was subsequently employed to identify gene clusters strongly associated with the MRGs and FRGs scores. Cluster analysis revealed no outlier samples. When the β-value was set to 11, the R² surpassed 0.85 (red line) and the mean connectivity approached zero, indicating optimal network consistency with a scale-free distribution and biological relevance (Figure 1C). Hierarchical clustering analysis identified 39 co-expression modules (excluding gray modules) (Figure 1D). Using a threshold of |cor| >0.3 and p < 0.05, four modules—MEblue, MEred, MEgreenyellow, and MEdarkmagenta—were found to be significantly correlated with the MRGs and FRGs scores (Figure 1E). In total, 3,098 genes within these modules were identified as key module genes associated with both the MRGs and FRGs scores.

Using the criteria of |log₂FC| >0.5 and p < 0.05, 234 DEGs, consisting of 129 upregulated and 105 downregulated genes, were identified between SCI and control samples in GSE109902 (Figures 2A,B). Overlapping 3,098 key modular genes with the DEGs revealed eight mitochondrial and ferroptosis-associated genes, designated as candidate genes: U4, AABR07048397.1, Hcrt, Rps18l1, Cdca2, SNORD115, RGD1560034, and Il1a (Figure 2C).

These candidate genes underwent enrichment analysis, resulting in 301 GO terms, including 280 biological processes (BPs), two cellular components (CCs), and 19 molecular functions (MFs), as well as 11 KEGG pathways (p < 0.05). The top 10 GO terms were predominantly related to “regulation of mitotic nuclear division,” “regulation of nervous system processes,” “temperature homeostasis,” “positive regulation of ion transport,” and “positive regulation of neutrophil extravasation,” among others (Figure 2D). KEGG analysis highlighted the involvement of candidate genes in “inflammatory bowel disease,” “graft-versus-host disease,” and “osteoclast differentiation,” among other pathways (Figure 2E).

Furthermore, a GeneMANIA database search identified 20 functionally similar genes (e.g., Hertrl, Mki67, Hertr2), and a co-expression network was constructed to elucidate their roles in various processes, including “regulation of interleukin 6 production,” “interleukin 6 production,” and “cytokine-mediated signaling pathways” (Figure 2F).

Two machine learning algorithms were employed to identify the most significant feature genes associated with SCI from the eight candidate genes. The Boruta algorithm highlighted five key genes (U4, AABR07048397.1, Hcrt, Rps18l1, Cdca2) based on feature importance (Figure 3A). In LASSO analysis, conducted at a lambda_min value of 0.05360385, eight feature genes were selected: U4, AABR07048397.1, SNORD115, RGD1560034, Hcrt, Rps18l1, Il1a, and Cdca2 (Figure 3B). An overlap analysis between the two algorithms revealed five common feature genes: U4, AABR07048397.1, Hcrt, Rps18l1, and Cdca2 (Figure 3C).

The expression level analysis revealed significant differences in the five genes between the SCI and control groups in the GSE109902 dataset (p < 0.05) (Figure 3D). Similarly, Hcrt and Cdca2 exhibited significant differential expression between groups in the GSE166009 dataset (p < 0.05) (Figure 3E). qPCR validation was performed on rat spinal cord tissue samples collected at days 2 and 7 post-SCI surgery, showing a marked increase in the expression of both Hcrt and Cdca2 following the procedure. These results align with the expression patterns observed in the GSE166009 validation dataset, confirming the reliability of the screening process (Figure 3F). Notably, both Hcrt and Cdca2 demonstrated elevated expression in SCI samples compared to control samples across the two datasets. However, no expression of U4, AABR07048397.1, or Rps18l1 was detected in the GSE166009 dataset. Consequently, Hcrt and Cdca2 were selected as biomarkers for subsequent analyses.

The GTEx database was utilized to assess the expression levels of the two biomarkers across various tissues under physiological conditions. Hcrt exhibited higher expression in the brain, particularly in regions such as the hypothalamus, nucleus accumbens, frontal cortex, amygdala, and spinal cord (Figure 4A). In contrast, Cdca2 was predominantly expressed in tissues like the testis, EBV-transformed lymphocytes, cultured fibroblasts, and esophageal mucosa (Figure 4B). Notably, Cdca2 was also detectable in the brain, with the highest expression observed in the spinal cord.

An ANN model was developed using the GSE109902 dataset to improve differentiation between SCI and control samples by integrating the expression of two biomarkers (Figure 4C). The ROC curve yielded an AUC value of 0.775, suggesting that the ANN model exhibited robust accuracy in predicting SCI (Figure 4D).

The potential biological mechanisms of the biomarkers were analyzed using GSEA, identifying 19 enriched pathways for Hcrt and 16 for Cdca2 (adj. p < 0.05 and |NES| >1). Shared enriched pathways included “ubiquitin-mediated proteolysis,” “endocytosis,” and the “neurotrophin signaling pathway” (Supplementary Tables S3, S4). Cdca2 was also associated with pathways such as “oxidative phosphorylation,” “cytokine-cytokine receptor interaction,” “inositol phosphate metabolism,” and “sphingolipid metabolism” (Figure 5A). Hcrt expression showed significant correlations with the “ErbB signaling pathway,” “insulin signaling pathway,” “mTOR signaling pathway,” and “primary immunodeficiency,” among others. The top five pathways for each biomarker were visualized according to pathway significance (Figure 5B).

Additionally, co-expressed gene networks for the biomarkers were constructed using the GeneMANIA platform, comprising 39.02% co-expression, 25.61% physical interactions, 21.99% predicted interactions, 8.15% pathways, 2.86% co-localization, and 2.37% shared protein domains (Figure 5C). Both biomarkers were linked to 20 functionally related genes (e.g., Mki67, Hcrtr1, Hcrtr2), participating in biological processes such as “hormone activity,” “chromosome segregation,” and “mitotic spindle” organization.

The infiltration levels of 22 immune cell types in each sample from the GSE109902 dataset were calculated using the CIBERSORT algorithm, as illustrated in Figure 6A. The Wilcoxon test revealed significant differences in T follicular helper cell levels between SCI and control samples, with lower levels observed in SCI samples (p < 0.05) (Figure 6B). Correlation analysis indicated that T follicular helper cells were negatively associated with naive CD4⁺ T cells (cor = −0.51, p = 0.004), activated mast cells (cor = −0.53, p = 0.003), and resting memory CD4⁺ T cells (cor = −0.63, p < 0.001) (Figure 6C). Biomarker-immune cell correlation analysis further demonstrated a positive association between T follicular helper cells and Hcrt, alongside a negative association with Cdca2 (Figure 6D).

The inflammatory response plays a key role in the pathophysiology of SCI, with immune factors critically influencing its development and progression. The relationship between two specific biomarkers and 200 IRGs was assessed. Analysis using the ssGSEA algorithm revealed a significant difference in IRG scores between SCI and control samples in the GSE109902 dataset (p < 0.05), with higher scores detected in the SCI samples (Figure 6E). Correlation analysis showed that Hcrt exhibited significant positive correlations with 10 IRGs (e.g., KLF6, RAF1, LY6E), while Cdca2 demonstrated significant positive correlations with six IRGs (e.g., IL1R1, ATP2B1, PSEN1) and significant negative correlations with two IRGs (BST2 and IL1B) (|cor| >0.3, p < 0.05) (Figure 6F). Notably, both Hcrt and Cdca2 were strongly positively correlated with PSEN1.

Analysis of the miRNet database identified 110 TFs regulating the biomarkers, among which seven TFs (CEBPB, FOXC1, GATA2, HNF4A, REST, and SP1) were recognized as co-regulators of Hcrt and Cdca2 within the TF-biomarker network (Figure 7A). Additionally, predictions from individual databases identified 54 miRNAs targeting the biomarkers and 533 lncRNAs targeting these miRNAs. Within the lncRNAs-miRNAs-mRNAs network, numerous interaction pairs were observed (Figure 7B). For instance, lncRNAs such as MALAT1, NEAT1, and XIST were predicted to regulate Hcrt expression via hsa-miR-30a-5p, while MALAT1 and XIST were also linked to the regulation of Cdca2 expression through miRNAs like hsa-miR-124-3p and hsa-miR-132-3p. These interactions suggest potential therapeutic targets for SCI.

Drug prediction via the Enrichr database identified five potential candidates associated with Cdca2: pinosylvin, zinc acetate dihydrate, hydroquinone, lucanthone, and dasatinib (Figure 7C). Notably, no drugs were predicted to interact with Hcrt. Molecular docking analysis demonstrated that all four drugs exhibited binding energies to Cdca2 below −5 kcal/mol (Table 2), suggesting stable interactions. The docking patterns of these four drugs with Cdca2 were subsequently visualized (Figures 7D–G).