Our study encompassed the GSE213240, GSE20907, and GSE33886 datasets from the Gene Expression Omnibus (GEO) (https://www.ncbi.nlm.nih.gov/). The GSE213240 dataset contained scRNA-seq information about two normal spinal cord tissue samples and six samples with varying degrees of SCI. The datasets GSE20907 and GSE33886, both related to SCI, were used as the training and external validation sets. GSE20907 includes 22 samples evenly distributed between experimental and control groups, while GSE33886 comprises 14 samples (8 control, 6 experimental), as shown in the Supplementary File ( Supplementary Table S1 ). Single-cell analysis was enabled by the “Seurat” package (9). Genes expressed in over three cells and cells exhibiting expression of at least 201 genes were further analyzed. Our study excluded cells that showed gene expression below 500 or above 6,000, as well as cells where mitochondrial RNA made up more than a quarter of the content. Dimensionality reduction was performed utilizing the top 4,000 highly variable genes. The initial 20 principal components were obtained via principal component analysis (PCA) through the JackStrawPlot function. Clustering was carried out employing the UMAP algorithm based on these components. To enhance data comparability and improve statistical efficiency, the Harmony package v.1.2.3 (parameters: θ=2, λ=1, dims=1–20) was employed to correct for batch effects. This procedure effectively eliminated technical artifacts and batch-related variability. During this process, a Seurat clustering resolution of 0.4 was applied. Cell sorts were annotated via the “SingleR” package (10) with MouseRNAseqData as a reference. Marker genes were found utilizing the FindAllMarkers function. MM cells were extracted for re-normalization, dimensionality reduction, as well as clustering. Annotation was performed as per their expression across experimental groups.

To unravel the joint effects of genes and their involvement in biological pathways, Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analyses were enabled by the “ClusterProfiler” package (11). Moreover, Gene Set Enrichment Analysis (GSEA) demonstrated gene sets possibly displaying no evident differences at the single-gene level but exhibiting strong enrichment trends at the pathway level. This approach offers insights into the overall functional behavior of genes under healthy and diseased conditions (12).

Pseudotime analysis enables inference of the dynamic progression of cells or samples through a biological process based on gene expression profiles. Monocle, among the most widely adopted tools, constructs developmental trajectories by embedding single-cell data into a low-dimensional space (13). This method facilitated our investigation of the dynamic state transitions of MM cells during SCI progression.

To investigate intercellular signaling, CellChatDB.mouse was the ligand-receptor interaction reference, and the “CellChat” package was applied to infer cell-cell communication networks (14). Interaction frequency and total interaction strength between cell types were visualized using circular plots. The netVisual_bubble function was further employed to display all signaling pathways among cell populations, allowing us to highlight pathways with strong interactions and identify potential signal senders and receivers. This analysis provides insights into intercellular communication, functional cooperation, and responses to external stimuli, offering a deeper understanding of the mechanisms underlying disease development.

Unlike conventional scRNA-seq analyses, hdWGCNA facilitates network analysis at the isotype level by leveraging long-read single-cell sequencing information (15). hdWGCNA was executed on MM cells, with the parameter fraction set to 0.05 to retain genes expressed in 5% or more of cells. The TestSoftPowers function presented the ideal soft-thresholding power (Minimum value corresponding to the dimensionless topological fitting index (R²≥0.8)). Identifying hub genes - highly interconnected genes within modules - is critical for downstream validation in co-expression network analysis. Therefore, module membership (kME) values were computed using the ModuleConnectivity function and used to identify hub genes within every module. Module-specific hub gene expression levels were measured on a per-cell basis, and modules highly associated with differences between pre-SCI and post-SCI MM cells were selected for further gene extraction.

Candidate disease-associated genes identified via hdWGCNA were further refined using multiple ML algorithms. Four classical algorithms were employed: Random Forest (RF), Survival-SVM, XGBoost, and the Boruta algorithm. These models were integrated to identify shared feature genes. The diagnostic potential of all feature genes and their related nomograms was analyzed through ROC curves with the area under these curves (AUC) as the evaluation standard. An AUC value of approximately 0.65 or higher was considered indicative of moderate to strong predictive performance (16). The GSE142426 dataset was the training set with validation performed on the GSE20907 dataset.

CIBERSORT, an immune cell infiltration analysis tool developed at Stanford University, was employed to analyze the percentages of 22 immune cell sorts among 24 GSE20907 samples (17). Subsequently, Pearson correlation coefficients revealed the relation of characteristic gene expression to immune cell subset infiltration.

RAW264.7 macrophages and BV-2 microglial cells were cultured in standard conditions (37°C, 5% CO₂). The control group (CON group) consisted of cells grown for 48 hours in high-glucose DMEM comprising 10% fetal bovine serum and 1% penicillin/streptomycin. In contrast, those treated for 48 hours with a medium containing 1 mg/mL lipopolysaccharide (LPS) served as the inflammation-mimicking model group (LPS group). Gene expression levels of GNGT2, EMP3, and SGPL1 in both groups were quantified via real-time qPCR (18). The Trizol reagent was employed for the isolation of total RNA, which was then converted into complementary DNA (cDNA) via the Seven All-in-one First Strand cDNA Synthesis Kit II (SM134). Amplification was enabled by the Seven SYBR Green qPCR MasterMix II (SM143). mRNA was quantified with β-Actin as the reference. Relative gene expression was calculated through the 2^−ΔΔCt approach. Statistical significance was evaluated via independent sample t-tests, with P < 0.05 denoting statistical significance.

The primer sequences were: CD206: forward 5’-GTTCACCTGGAGTGATGGTTCTC-3’, reverse 5’-AGGACATGCCAGGGTCACCTTT-3’ IL-10: forward 5’-CGGGAAGACAATAACTGCACCC-3’, reverse 5’-CGGTTAGCAGTATGTTGTCCAGC-3’ TGF-β: forward 5’-TGATACGCCTGAGTGGCTGTCT-3’, reverse 5’-CACAAGAGCAGTGAGCGCTGAA-3’ β-Actin: forward 5’-ACTGCCGCATCCTCTTCCT-3’, reverse 5’-TCAACGTCACACTTCATGATGGA-3’ GNGT2: forward 5’-GCTGTTGAGGATGGAGGTG-3’, reverse 5’-AAGGGATTCTTGTCTTCTGG-3’ SGPL1: forward 5’-GGAAAGCCTCAGGAGCTGTGTA-3’, reverse 5’-CTGCCTCTAACTTCCGCAATCC-3’ EMP3: forward 5’-GTTCCAACTCTACACCATGCGG-3’, reverse 5’-ATCTCCTCGGTGTGGATGGCAT-3’.

Our data analysis was enabled by R 4.4.1. Experimental data were statistically studied via GraphPad Prism (GraphPad Software, La Jolla, CA). Inter-group difference assessment was completed through independent sample t-tests. Asterisks denote the following: *** P < 0.001, ** P < 0.01, * P < 0.05.

Following quality checks and batch effect correction, 58,674 high-quality cells expressing 19,414 genes were retained. Among these, the SCI group contributed 46,650 cells expressing 18,874 genes, while the normal group provided 12,024 cells expressing 18,158 genes. Unsupervised clustering identified 21 distinct cell clusters ( Figures 2A, B ). Subsequent annotation revealed 16 well-defined cells like granulocytes, oligodendrocytes, macrophages/microglia, fibroblasts, monocytes, B, T, endothelial, NK, dendritic and epithelial cells, neurons, astrocytes, erythrocytes, adipocytes, and hepatocytes ( Figure 2C ). The proportions of each cell type were compared across the control and experimental cohorts ( Figure 2D ). The percentage of MM cells was notably increased following SCI. Differential gene expression analysis was subsequently performed with results visualized via volcano plots ( Figure 2E , Supplementary Table S2 ).

An in-depth MM cell analysis was further carried out. Following repeated UMAP analysis, six cellular subpopulations were identified via manual annotation, with the classification being determined by the differential gene expression observed across the two tissue varieties ( Figure 3A ). Four of the six subpopulations (MSCI1, MSCI2, MSCI3, and MSCI4) exhibited elevated expression levels following SCI. Therefore, subsequent investigations were directed toward these four subpopulations. Subsequent statistical analysis revealed significant differences in the MSCI1 and MSCI3 subclusters ( Figure 3B ), suggesting the crucial roles of these two cell populations in SCI onset and progression. A heatmap further presented the expression patterns of subtype-specific marker genes within MM cells ( Figure 3C ). To unravel the possible molecular mechanisms involving MSCI1 and MSCI3 in SCI, GSEA was performed on significantly altered gene sets. The results indicated substantial activation of pathways related to cytoplasmic translation, inflammatory response, precursor metabolite, energy production, and response to external stimuli ( Figure 3D ). Additionally, GO and KEGG enrichment analyses revealed significant involvement in biological processes such as regulation of the nerve growth factor pathway, adaptive immune response, phagocytic apoptosis signaling, inflammatory response, and in the pathogenesis of diseases, including glioma ( Figures 3E, F ).

Pseudotime trajectory analysis revealed the dynamic changes in MM cells during the advancement of SCI. This analysis clearly illustrated that, with disease progression, the developmental trajectory of these cells can be delineated into seven distinct phases ( Figure 4A ). The inferred pseudotemporal trajectory was constructed based on integrated functional metrics ( Figure 4B ). Upon annotation of the cell subpopulations, it became evident that the MSCI1 and MSCI3 subsets play prominent roles throughout disease progression ( Figure 4C ).

Given these findings, MSCI1 and MSCI3 were selected for further in-depth investigation. To systematically explore the interactions between MSCI1, MSCI3, and other MM cell subpopulations, our study employed the CellChat algorithm, which integrates the expression levels of signaling ligands, receptors, and their cofactors to infer intercellular communication networks. The analysis revealed robust and extensive interactions between MSCI1, MSCI3, and other cellular subtypes ( Figure 4D ). Notably, MSCI1 and MSCI3 exhibited enhanced secretory activity and demonstrated greater outgoing interaction strengths compared to other subgroups ( Figure 4E ). The investigation of ligand-receptor pairing between MSCI1/MSCI3 and different cells uncovered several highly expressed molecular axes, including CD99-CD99, SPP1-CD44, and CCL3-CCR1 ( Figures 4F–I ). Further analysis indicated that signaling pathways associated with immune response, inflammatory processes, and cell migration, particularly those involving CD99, CCL3, and SPP1, were markedly upregulated in communications involving MSCI1 and MSCI3 ( Figures 4J–L ).

To ensure biological relevance in the WGCNA, a soft-thresholding power of 12 was used to approximate a scale-free topology ( Figure 5A ). Upon achieving optimal scale - Independence and mean connectivity, 14 gene modules were noted ( Figure 5B ). Within each module, highly connected genes - termed hub genes - were determined by calculating the link of individual genes to the module eigengene (KME value). Notably, the turquoise module was considerably enriched in MSCI1 and MSCI3 ( Figure 5C ). The expression activity of hub genes within each module was visualized across individual cells ( Figure 5D ). Furthermore, our study examined the inter-subgroup expression patterns of each module ( Figure 5E ) and the inter-module correlation structure ( Figure 5F ). The turquoise module exhibited a positive correlation with the purple and yellow modules, and prior cell interaction analysis demonstrated that MSCI1 and MSCI3 were strongly associated with oligodendrocytes and fibroblasts. Figures 4E, F indicate that Ptn-Ncl and Cd99-Cd99 are highly expressed. The expression of the Ptn and Cd99 genes was subsequently validated ( Supplementary Figure S1 ), thereby corroborating this observation and providing insight into potential novel signaling pathways between these cell types.

In the previous analysis, 200 hub genes were identified within the turquoise module. Building upon these results, four ML algorithms were further employed to mine for SCI-associated central genes using the GSE20907 dataset. In the XGBoost algorithm, parameters were set as max_depth=5 and eta=0.2, and the top 14 genes were identified based on the feature importance evaluation ( Figure 6A ). Subsequent SVM-RFE analysis indicated that prediction efficacy was greatest with nine genes ( Figure 6B ). The RF model displayed the lowest misclassification rate when the number of decision trees reached 53, identifying 12 candidate genes at this point ( Figure 6C ). Using parameters (perc = 100, alpha = 0.01, max_features < 0.3) in Boruta analysis, 19 genes that significantly contributed to the prediction task were effectively screened out ( Figure 6D ). Venn diagrams present the overlapping genes identified through the four algorithms, resulting in the selection of three intersecting genes: EMP3, GNGT2, and SGPL1 ( Figure 6E ).

After comparing spinal cord tissue datasets before and after injury, EMP3, GNGT2, and SGPL1 were markedly upregulated following SCI ( Figures 6F–H ). Integration of spatial transcriptomic data with single-cell gene expression patterns demonstrated the main enrichment of these three genes in the MSCI1 and MSCI3 clusters ( Figures 6I–L ) ( Supplementary Figure S2 ). To evaluate their diagnostic potential, ROC curve analyses were carried out for EMP3, GNGT2, and SGPL1, which yielded AUC values of 0.951, 0.917, and 0.896 ( Figures 6M–O ). Validation using an external dataset (GSE33886) demonstrated AUCs of 0.671, 0.640, and 0.742 ( Figures 6P–R ).

Infiltration of immune cells, which is crucial for inflammatory reactions, significantly and independently influences the pathophysiology of SCI. CIBERSORT, a computational method for estimating the relative abundance of immune cell sorts within a single sample, was employed via the LM22 signature matrix, which profiles 22 immune cell subsets with diverse phenotypes and functional states. Comparative analysis revealed marked differences in macrophage infiltration between normal and SCI tissues, suggesting pronounced macrophage involvement in SCI samples ( Figures 7A, B ). The intercorrelations among all immune cell types were further analyzed ( Figure 7C ). M2 macrophages were found to be positively correlated with CD8 T cells, while exhibiting a negative correlation with plasma cells, neutrophils, and activated NK cells. In addition, activated NK cells, plasma cells, and naïve CD4 T cells demonstrated positive correlations. These observations provide preliminary insights into the interactions among immune cell subsets. Using expression data of the genes identified via ML and the immune cell infiltration matrix, our study examined the associations between SGPL1, GNGT2, EMP3, and immune cells. Notably, all three genes exhibited strong correlations with M2 macrophages ( Figure 7D ), providing theoretical support for targeting M2 macrophages in future studies on therapeutic strategies for SCI recovery.

In the RAW264.7 macrophage in vitro inflammation model, the mRNA levels of CD206 (P = 0.001), IL-10 (P = 0.0033), TGF-β (P = 0.0233), GNGT2 (P = 0.0018), SGPL1 (P = 0.0082), and EMP3 (P = 0.0083) were markedly upregulated in the LPS-treated group. Similarly, in the BV-2 microglial inflammation model, elevated expression trends were also observed for CD206 (P = 0.0001), IL-10 (P = 0.003), TGF-β (P = 0.0078), GNGT2 (P = 0.0163), SGPL1 (P = 0.0101), and EMP3 (P = 0.0426). The corresponding qPCR results are presented in Supplementary Tables S3–S6 , and Figure 8 .