Datasets GSE151371, GSE45006, and GSE162610 were sourced from the GEO (Gene Expression Omnibus) database (https://www.ncbi.nlm.nih.gov/geo/) (Edgar et al., 2002). The GSE151371 dataset includes gene expression profiles of peripheral white blood cells (WBCs) along with clinical data from 38 SCI patients and 10 healthy individuals (Kyritsis et al., 2021). Probes were mapped to their respective gene symbols using the annotation information from the GPL20301 [Illumina HiSeq 4,000 (Homo sapiens)] platform. The GSE45006 microarray expression dataset for SCI rats was generated using the GPL1355 [Rat230_2] Affymetrix Rat Genome 230 2.0 Array platform. In this dataset, SCI was induced in rats using an aneurysm clip impact-compression injury model on the thoracic spinal cord (T7). This dataset includes the Sham group (n = 4) and the SCI group, with samples collected at 1, 3, and 7 days post-injury (n = 4 per time point). Although GSE151371 and GSE45006 differ in species and tissue origin, GSE45006 was incorporated as a supportive external validation dataset to assess the consistency of glycolysis-related gene expression patterns at the transcriptomic and pathway levels, rather than to establish direct tissue or species equivalence. The scRNA-seq dataset GSE162610 includes data from normal and injured spinal cords at 1, 3, and 7 days post-injury of wild-type mice. This dataset was used to characterize the cellular distribution of candidate hub genes within the spinal cord microenvironment, thereby providing complementary cell-type–specific context for the transcriptomic findings derived from bulk RNA-seq analysis.

The GeneCards database (https://www.genecards.org/) provides extensive gene information for humans (Stelzer et al., 2016). By searching with the keyword “glycolysis”, we selected genes classified as ‘Protein Coding’ and with a relevance score exceeding 1.00, thereby identifying glycolysis-related genes (GRGs). This search resulted in the identification of 1,572 GRGs, which are listed in Supplementary Table S1.

To quantify glycolytic transcriptional activity in each sample, a predefined set of glycolysis-related genes was extracted from the transcriptomic data, and principal component analysis (PCA) was performed using the expression matrix of these genes. PCA identifies orthogonal linear combinations of genes that successively maximize variance within the selected gene set. In this study, the first two principal components (PC1 and PC2) captured the dominant axes of coordinated transcriptional variation among glycolysis-related genes. As PC1 and PC2 together explained the majority of the variance within this gene set, a composite glycolysis score was constructed by summing the PC1 and PC2 scores for each sample. This PCA-based score was used as an empirical measure to summarize overall glycolytic transcriptional activity, consistent with previous studies (Sun et al., 2021). This score summarizes coordinated transcriptional variation within GRGs.

In this study, WGCNA was used to identify gene modules strongly associated with glycolysis, utilizing the “WGCNA” package (version 1.73) in R software (Langfelder and Horvath, 2008). Hierarchical clustering was applied to the study samples to identify and exclude outliers. A soft-thresholding power of β = 12 was chosen to build a scale-free network and compute the topological overlap matrix. The correlation between each module and the glycolysis score was evaluated. Modules with a high correlation coefficient were identified as key modules strongly associated with glycolysis. A correlation coefficient exceeding 0.5 and a significance level of P < 0.05 were set as the criteria for module-phenotype correlation. Genes from the key module were compared with glycolysis-related genes and DEGs, leading to the identification of a set of glycolysis-associated DEGs (GR-DEGs).

The raw GEO data was normalized using the ‘NormalizeBetweenArray’ function in R. Differentially expressed genes (DEGs) between the HC and SCI groups in the GSE151371 dataset were identified using the “Limma” package, which was found to be the most efficient tool. A log fold change (LogFC) threshold of >0.5 and an adjusted P-value <0.05 were applied to determine gene differential expression.

In this study, we utilized the Bagged Trees, Bayesian Networks, Support Vector Machine (SVM), Random Forest, XGBoost, Boruta, Learning Vector Quantization (LVQ), and 1,000 iterations of 10-fold cross-validation with Least Absolute Shrinkage and Selection Operator (LASSO) algorithms to pinpoint the key genes linked to glycolysis. Genes that were consistently recognized as feature genes by at least 7 different algorithms, based on their classification efficacy, were classified as key hub genes involved in the interaction between SCI and glycolysis. The interactions between these crucial intersecting genes were subsequently visualized using the “Upset” application from the R package. It should be noted that these machine learning algorithms were used for feature selection and robustness assessment rather than for constructing a single optimized predictive classifier.

To investigate the regulatory pathways and biological functions linked to the two hub genes, single-gene GSEA analysis was conducted. A significance threshold of adjusted P-value <0.05 was applied. The top five enriched pathways for both upregulated and downregulated genes were displayed separately, based on their enrichment scores.

Transcription factors (TFs) are proteins that interact with DNA at specific sequences to control the process of transcription. Through their ability to recognize unique DNA motifs, TFs influence chromatin configuration and transcriptional activity, orchestrating the regulation of gene expression across the genome (Lambert et al., 2018). MicroRNAs (miRNAs) are a class of short, endogenous non-coding RNAs that regulate gene expression by either degrading target mRNAs or preventing their translation (Bartel, 2004). To elucidate the dysregulation of gene expression in diverse physiological and disease contexts, it is essential to explore the intricate transcriptional regulatory networks that govern the interactions between TFs and miRNAs (Bartel, 2004; Zaborowski and Walther, 2020; Iwama et al., 2011). NetworkAnalyst (http://www.networkanalyst.ca) is an online tool designed for advanced meta-analysis of gene expression data (Xia et al., 2015). IRS1 and UBTD1 were analyzed using NetworkAnalyst to construct a TF–miRNA coregulatory network. Regulatory interaction data curated from the literature were gathered from RegNetwork (http://www.regnetworkweb.org/) (Liu et al., 2015). The relevant results were visualized using Cytoscape.

Using the NetworkAnalyst platform, compounds that interact with the hub genes were retrieved from The Comparative Toxicogenomics Database (CTD), and a protein-chemical interaction network was then built based on these interactions. Cytoscape was employed to visualize the relevant results.

The analysis of single-cell RNA sequencing data was mainly carried out using the “Seurat (v 5.1.0)” package (Hao et al., 2021). To maintain data quality, cells with fewer than 200 detected genes or more than 20% mitochondrial gene content were excluded. Batch effects were corrected using the “Harmony” package (v.1.2.1) (Korsunsky et al., 2019). Data normalization was performed using the ‘lognormalization’ function, and highly variable genes were selected with the ‘FindVariableFeatures’ function. We applied the tSNE algorithm to classify cell types, generating a unified map that reveals structural patterns across various scales, making it especially valuable for high-dimensional data. Clustering was performed with the Louvain algorithm, and marker genes specific to each cell subpopulation were identified through the ‘FindAllMarkers’ function. Cell type annotation was initially carried out automatically using the “SingleR” package.

The animals used in this study were sourced from the Laboratory Animal Center at Xi’an Jiaotong University (Animal License No. SCXC [Shan] 2023–002). All experimental protocols were approved by the Animal Ethics Committee of Yan’an University and carried out in strict compliance with both institutional and national regulations governing the ethical treatment and use of animals. A total of 12 female Sprague Dawley (SD) rats (8 weeks old) were randomly assigned to two groups: the Sham group and the SCI group. In the Sham group, a laminectomy was performed without inducing spinal cord injury. The SCI model was created following a modified version of the Allen method, as outlined in the existing literature (Liu et al., 2021). Briefly, the animals were anesthetized by intraperitoneal injection of 1% pentobarbital sodium (40 mg/kg, intraperitoneal injection), and then positioned on a surgical table. The T10 vertebra was exposed by removing the lamina using bone cutting forceps. To induce the injury, a 10 g weight was dropped from a height of 3 cm directly onto the exposed spinal cord. To prevent infection, penicillin (Dingjian, Sichuan, China) was administered subcutaneously for three consecutive days at a daily dose of 2.5 mg/kg. Bladder voiding was manually facilitated twice daily through the application of abdominal pressure until the animals were able to urinate spontaneously. Seventy-two hours following SCI induction, the SD rats were anesthetized through an overdose of 1% pentobarbital sodium administered via intraperitoneal injection. A segment of spinal cord tissue, 1.5 cm in length from the center of the injury site, was carefully excised for further analysis.

Total RNA was isolated from spinal tissues using the TRIzol Total RNA Purification Kit (5003050, Simgen, China), and cDNA synthesis was carried out with the NovoScript® Plus All-in-one 1st Strand cDNA Synthesis SuperMix (E047-01A, Novoprotein, China). Target gene expression was quantified using the NovoStart® SYBR qPCR SuperMix Plus (E096-01A, Novoprotein, China). The relative expression levels were normalized to β-actin and calculated using the 2^−ΔΔCt method. The primer sequences are presented in Table 1.

Proteins were extracted from spinal tissues, and their concentrations were quantified using the BCA method. The extracted proteins were separated by SDS-PAGE and subsequently transferred to a PVDF membrane (Merck, Germany). The membrane was blocked with 5% skimmed milk for 1 h to prevent non-specific binding. Primary antibodies were applied to the membrane overnight, including rabbit anti-UBTD1 (1:1,000 dilution, DF15663, Affinity) and rabbit anti-β-Tubulin (1:4,000 dilution, YM8332, ImmunoWay). Following washing steps, the membrane was incubated with HRP-conjugated goat anti-rabbit IgG (H + L) (1:10,000 dilution, GB111738-100, Servicebio) for 2 h at room temperature. After additional washing, the protein bands were visualized using ECL reagent, and the results were captured and documented using a chemiluminescence imaging system.

All analyses were conducted using R version 4.2.1. Data are presented as means ± standard error of the mean (SEM) from at least three independent experiments and were analyzed using GraphPad Prism 9.3.1. The Wilcoxon rank-sum test and Welch’s t-test were employed to compare the expression differences between unpaired samples from the two groups. P-values were interpreted according to the following scale of significance: asterisks denote increasing levels of significance (*P < 0.05, **P < 0.01, ***P < 0.001), while ‘ns’ denotes non-significance.

PCA was performed based on the expression of GRGs to calculate the glycolysis score for each sample. The glycolysis score was significantly higher in SCI samples compared to HC samples (Figure 2A). Using the GSE151371 dataset from GEO, clustering analysis was carried out, and a soft-thresholding power of β = 12 was chosen, ensuring an R ² value greater than 0.9 and a relatively high mean connectivity (Figure 2B). Under the parameter settings of minModuleSize = 100 and mergeCutHeight = 0.25, six modules were identified (Figures 2C,D). The MEbrown module demonstrated the strongest correlation with the score, exhibiting a correlation coefficient of r = 0.76 and a P-value of 3.42e-10, as shown in Figure 2E. Moreover, the GS and MM values in the brown module were strongly correlated, suggesting that this module exhibited the most significant association with glycolysis (Figure 2F). Based on this, we selected the MEbrown module, which consists of 704 genes, as the key module for further analysis.

Differential expression analysis of the GSE151371 dataset identified a total of 1,138 DEGs, with 707 genes being upregulated and 431 genes downregulated, as depicted in the volcano plot (Figure 3A). By intersecting the key WGCNA module genes, DEGs, and GRGs, a total of 13 GR-DEGs were identified (Figure 3B). Additionally, a heatmap was generated to illustrate the expression patterns of the 13 GR-DEGs. The results showed that ALK and IRS1 were downregulated, while PPARG, SLC1A3, GCKR, KL, UBTD1, GGH, GPR87, ANO5, SCO2, DHCR7, and HMGB3 were upregulated in SCI (Figure 3C). The positions of the 13 GR-DEGs on the chromosomes were mapped using the “RCircos” package for visualization (Figure 3D) (Zhang et al., 2013). The correlation analysis of the 13 GR-DEGs revealed significant interconnections among these genes, indicating a strong mutual correlation. As shown in Figures 3E–G, the correlation heatmap illustrates the overall relationships among the 13 GR-DEGs, while Figures 3F,G display the strongest positive and negative correlations, respectively.

The GSE151371 dataset was used to apply machine learning, which facilitated the identification of 13 GR-DEGs associated with the SCI. The intersection of findings from various machine learning algorithms, including Lasso-Logistic regression with 1,000 iterations (Figure 4A), the LVQ algorithm (Figure 4B), Boruta (Figure 4C), Bagged Trees (Figure 4D), Random Forest (Figure 4E), Bayesian Networks (Figure 4F), SVM (Figure 4G), and XGBoost (Figure 4H), revealed six genes that are significantly associated with the interaction between SCI and glycolysis (Figure 4I). These genes are ALK, GGH, IRS1, PPARG, SLC1A3, and UBTD1. To validate the six hub genes, the GSE45006 dataset was utilized, and the analysis demonstrated that only Irs1 and Ubtd1 exhibited significant differential expression between the SCI and Sham groups (Figure 4J).

The relative abundance of 28 immune cell types was quantitatively evaluated using the ssGSEA method, and comparisons were conducted between the SCI and HC groups based on the peripheral blood transcriptomic data. Significant disparities in the abundance of 20 immune cell types between the two groups were identified through the Wilcoxon rank-sum test. Specifically, the SCI group showed significantly higher infiltration levels of activated dendritic cells, CD56dim natural killer cells, central memory CD4 T cells, gamma delta T cells, macrophages, monocytes, neutrophils, and regulatory T cells. In contrast, the infiltration levels of activated B cells, activated CD8 T cells, CD56bright natural killer cells, effector memory CD8 T cells, immature B cells, memory B cells, natural killer cells, natural killer T cells, T follicular helper cells, type 1 T helper cells, type 17 T helper cells, and type 2 T helper cells were significantly lower (Figure 5A). Further correlation analysis revealed a significant positive correlation between IRS1 and the abundance of activated CD8 T cells, CD56bright natural killer cells, gamma delta T cells, and T follicular helper cells in SCI (Figure 5B). Additionally, UBTD1 showed a significant positive correlation with activated dendritic cells, monocytes, and neutrophils, while demonstrating a negative correlation with activated CD8 T cells, effector memory CD8 T cells, immature B cells, memory B cells, and type 2 T helper cells in SCI (Figure 5B). These findings reveal distinct alterations in immune cell populations within systemic circulation of SCI patients, suggesting a potential link between systemic immune responses and the identified hub genes.

A single-gene GSEA analysis was performed to explore the biological functions and associated pathways of the two candidate hub genes. The upregulation of IRS1 is linked to immune-related pathways, platelet-specific genes, and translation initiation processes (Figure 6A). The downregulation of IRS1 is associated with immune functions, neutrophil activity, and the regulation of myeloid cells (Figure 6B). These findings underscore the diverse biological functions of IRS1 in regulating both immune and cellular activities. The upregulation of UBTD1 is associated with immune response regulation, neutrophil activity, and cancer cell proliferation, particularly in the context of infection and cancer (Figure 6C). The downregulation of UBTD1 is linked to immune response mechanisms, particularly cytokine signaling and interferon pathways, while also being linked to diseases like Sezary syndrome and melanoma (Figure 6D). These results emphasize the broad and complex roles of IRS1 and UBTD1 in regulating immune and cellular activities across different biological contexts.

The TF–miRNA coregulatory network for IRS1 is composed of 48 nodes and 47 edges. The analysis indicates that IRS1 is modulated by 23 TFs and 23 miRNAs (Figure 7A). The TF–miRNA coregulatory network for UBTD1 consists of 31 nodes and 30 edges. The findings demonstrate that UBTD1 is influenced by 5 TFs and 25 miRNAs (Figure 7B).

Protein-chemical interaction networks are crucial for understanding the mechanisms underlying disease development and can aid in the drug discovery process. In this study, protein-chemical interaction networks for the candidate hub genes were constructed, consisting of 81 nodes and 84 edges. The compounds that interact with both IRS1 and UBTD1 include Aflatoxin B1, Atrazine, Testosterone, Aspirin, and Calcitriol (Figure 7C).

To explore the potential involvement of IRS1 and UBTD1 in SCI at the single-cell resolution, we performed an analysis of the scRNA-seq data from the GSE162610 dataset. A total of 28,936 cells derived from C57BL/6 mouse spinal cord samples were processed, followed by dimensionality reduction, clustering, and t-SNE visualization of the dataset. Through analysis, 31 distinct cell clusters were characterized, and marker genes were identified using the “FindAllMarkers” package (Figure 8A). Using these marker genes, a total of 162,610 cells were classified into nine distinct cell types: astrocytes, endothelial cells, epithelial cells, granulocytes, macrophages, microglia, monocytes, neurons, and oligodendrocytes (Figure 8B). Irs1 was predominantly expressed in endothelial and epithelial cells (Figures 8C,D). We found Ubtd1 expression across multiple cell types, with elevated levels found in endothelial cells and microglia (Figures 8E,F).

To experimentally validate the roles of IRS1 and UBTD1 in SCI, we conducted RT-qPCR and Western blot analyses to assess the mRNA and protein levels of IRS1 and UBTD1 in spinal tissue from the SCI group compared to the Sham group. qPCR validation revealed a significant upregulation of Ubtd1 in the SCI group versus the Sham group, consistent with the bioinformatic predictions. In contrast, Irs1 mRNA levels did not exhibit a statistically significant difference (Figure 9A). Furthermore, the Western blot analysis demonstrated a notable increase in UBTD1 expression in the spinal tissue of the SCI group compared to the Sham group (Figure 9B), confirming its robust dysregulation at both transcriptional and translational levels. Consequently, based on the cross-validation results, UBTD1 was identified as the primary hub gene for further discussion.