Exploration of Hub Genes and Immune Cell Infiltration Characteristics Associated With Spinal Cord Injury in Mice

Introduction

Spinal cord injury (SCI) is a serious disease caused by direct or indirect external forces on the spinal cord. In addition to the devastating impact of the primary injury, secondary injury is a huge obstacle to the recovery of neurological function.¹ Immune cell infiltration and inflammatory cytokine release driven by multiple causes after SCI amplifies the secondary damage and impedes neurological recovery.^2,3 Therefore, an in-depth understanding of the molecular mechanisms and immune cell infiltration associated with the development of secondary injury after SCI is critical for identifying novel therapeutic targets for neurological recovery.

The developmental process of SCI and its subsequent recovery of neurological function is associated with various biological processes mediated by immune cells. In the absence of the role of microglia and astrocytes inherent in the spinal cord, the infiltrated neutrophils often cause and enhance the inflammatory response and aggravate the degree of injury, which affects the synaptic regeneration of neurons and hinders the recovery of neurological function in the spinal cord.^2–4 Interestingly, a specific number of neutrophils can stimulate neuronal axon regeneration and promote neurological recovery.⁵ Macrophages of monocyte and microglia origin infiltrate the center of the injury and are activated upon exposure to multiple inflammatory factors at the site of injury, activated macrophages can result in axonal regression and death.^6,7 In addition, increased M1-type pro-inflammatory macrophage transformation under inflammatory conditions can exacerbate the local inflammatory response and inhibit nerve function recovery.⁸ The inhibition of M1 macrophage differentiation has been shown to reduce the local inflammatory response to spinal cord injury and promote neurological recovery.⁹ In adjunction, infiltrating T lymphocytes play an indelible role in the recovery of injury and neurological function; it has not only direct toxic effects on neurons and glial cells but also exerts pro-inflammatory effects through the release of pro-inflammatory cytokines.¹⁰ Another study demonstrated that CD4 T lymphocyte infiltration was significantly reduced and that mice recovered better in terms of neurological function.¹¹ These studies together demonstrated that immune cells are essential for spinal cord injury. Although the development of SCI and the subsequent recovery of neurological function has been studied widely, there are still many potential therapeutic targets that need to be explored further.

In this study, we screened the differential genes. Moreover, we conducted the weighted correlation network analysis (WGCNA) to screen out the candidate key genes, and obtained the most relevant genes to spinal cord injury according to least absolute shrinkage and selection operator (LASSO) and Receiver Operating Characteristic Curve (ROC) algorithms. Furthermore, we analyzed the infiltration of relevant immune cells and the association of each key gene (2010004M13Rik, Cdkn1c, Nox4, and Gpr101) with immune cells. The results of the bioinformatic analysis were validated in animal models in order to identify new hub genes for spinal cord injury and potential therapeutic targets for neurological recovery (Figure 1). This study not only provide new ideas for the treatment of spinal cord injuries, but may also advance the development of precision medicine and personalized treatment.

Figure 1 Flowchart of the present study.

Methods

Microarray Data Download and Processing

Two gene datasets related to spinal cord injury from mus musculus were downloaded from the public database GEO (Gene Expression Omnibus, http://www.ncbi.nlm.nih.gov/geo/): GSE5296 and GSE47681. Both sequencing platform files are GPL1261 (Affymetrix Mouse Genome 430 2.0 Array). The GSE5296 database consisted of 6 sham-operated control samples and 9 mice with spinal cord injury, while the GSE47681 dataset comprised 9 sham-operated mice and 7 with spinal cord injury. The probe information was converted to the corresponding gene symbols according to the annotation information on the respective platform files of the database. Then, the two datasets were normalized using the “SVA” (version 3.50.0) package in the R program (version 4.1.0) to obtain the training group dataset.

Identification of DEGs

The “limma” package (version 3.58.1) in the R program was used to analyze the DEGs between the spinal cord injury group and the healthy control group in the training data set, and the screening criteria were |logFC| > 1, P-value < 0.05, and adjusted P-value < 0.05. The “pheatmap” (version 1.0.12) and “ggplot2” (version 3.2.2) packages were used to visualize the differential gene results and create DEGs heat maps and volcano maps. The details of DEGs can be seen in the Supplementary Table 1.

DEGs Enrichment Analysis and Protein-Protein Interaction (PPI) Network

Enrichment analysis of differential gene pathways and processes was performed using the Metascape online tool (www.metascape.org; version 3.5.20230501) using the following ontology sources: GO Biological Processes, KEGG pathway, Reactome gene Sets, CORUM, and WikiPathways, while all genes in the training dataset comprised the enrichment background. A P-value < 0.01, a minimum count of 3, and an enrichment factor > 1.5 (the enrichment factor is the ratio between the observed counts and the counts expected by chance) were set as filtering conditions, and the genes were grouped into clusters based on their similarity. We also constructed a PPI network using the Cytoscape online tool (https://cytoscape.org/; version 3.10.1) on DEGs.

Construction of Co-Expression Networks

WGCNA is a method for constructing gene co-expression networks based on gene expression data. It was performed on the training set expression profile data using “WGCNA” (version 1.72–1) in the R software. First, the data set was sorted and filtered to check the completeness of the data. Then, the “pickSoftThreshold” function was used to select a desirable soft threshold (b) to represent the strong correlation between gene networks, truncate the disordered neighborhoods between genes, and finally transform them into a topological overlap matrix (TOM) to measure the network connectivity between genes and analyze the dynamic tree cut by hierarchical clustering, such that the genes with similar differential expression can be included in the same module. Next, the modules were combined with phenotypic data to calculate gene significance (GS) and module correlation (MM), which were used to measure the importance of genes and clinical information as well as analyze the correlation between modules and models.

Identification of Hub Genes

Hub genes are the genes with the highest intermodule connectivity. GS > 0.5 and MM > 0.8 were set as the conditions to filter out the module signature genes. Then, we used the “venn” (version 1.11) package of the R program to intersect the module signature genes with the differential genes and used the LASSO regression algorithm to identify the final hub genes. Subsequently, BOXPLOT was used to analyze the expression levels of hub genes in the SCI and the healthy control groups. ROC curves were plotted using the “pROC” (version 1.18.4) package to assess the accuracy of the Hub genes as disease markers.

Immune Infiltration Analysis

Immune cell infiltration analysis was performed in SCI mice using the CIBERSORT method at a threshold of P-value < 0.05 to select the samples with high reliability and obtain the relative content of 24 immune cell infiltrates in each sample. The correlations between immune cells were analyzed using the “corrplot” (version 0.92) software package. The “vioplot” (version 0.4.0) software package was used to analyze and visualize the differences between the immune cells of the experimental and control samples. The correlations between key genes and immune cells were analyzed using the “reshape2” (version 1.4.4), “ggpubr” (version 0.6.0) and “ggExtra” (version 0.10.0) packages with the standard P-value < 0.05.

Animal Model

BALB/c male mice (20–24 g, 6–8 weeks) provided by the Nantong University Animal Center were used for the experiments. All procedures were conducted in accordance with the “Guiding Principles in the Care and Use of Animals” (China) and were approved by the Laboratory Animal Ethics Committee of Nantong University (S20230315-004). The animals were housed in 12 h light/dark cycles, constant temperature (25°C) and humidity (40%), with free access to food and water. The mice were anesthetized by intraperitoneal injection of Pentobarbital Sodium (40mg/kg), and the spinal cord injury group was exposed by laminectomy at T9-T10. Subsequently, the spinal cord was impacted at T9-T10 with a blunt head using a 10 g keratin needle dropped freely at 5 cm, and the establishment of the spinal cord injury model was confirmed when the mice showed bilateral hind limb twitching and tail wagging. In the control group, only laminectomy was applied. (Figure 2) After successful modeling, the incision was sutured, and benzylpenicillin (20000UI/kg, ip) and buprenorphine (0.1 mg/kg, sc) were given to prevent infection and analgesia. The bladder was emptied artificially (twice a day). All exper

Figure 2 A: Spinal cord injury model; (B) control group with simple laminectomy model.

Tissue RNA Extraction and qRT-PCR

The two groups of mice were executed 3 days after surgery, and the spinal cord tissues (0.15 cm above and below the injury site) were removed from T9-T10, crushed by a tissue grinder, and total tissue RNA was extracted from each sample by adding the corresponding volume of Trizol (Vazyme) according to the weight, and subsequently, qRT-PCR was performed using cDNA reverse transcription reagent kit (1st Strand cDNA Synthesis Kit (by stem-loop), Vazyme) to reverse transcribe the RNA into cDNA, Real-time RT-qPCR was done using gene-specific primers and SYBR^® (2×SYBR Green qPCR Master Mix (None ROX), Servicebio). The relative expression of mRNA was calculated by the 2^−ΔΔCT method against GADPH as an internal reference.

Statistical Analysis

All experiments contained at least three biological replicates, and all results were illustrated as mean ± standard error (SD). A statistical significance was detected using Student’s t-test or one-way analysis of variance (ANOVA) by SPSS software (version 19; IBM Corp., USA). Bar charts were created using the GraphPad Prism software (version 9.0; GraphPad Prism, USA). “*” indicates the P-value < 0.05. “**” indicates that the P-value < 0.01. “***” indicates that the P-value < 0.001.

Results

The Differentially Expressed Genes in Mouse With Spinal Cord Injury

By setting |logFC| > 1, P-value < 0.05, adjusted P-value < 0.05 as the screening conditions, we identified 522 DEGs between the spinal cord injury and the healthy control groups, including 440 up and 82 downregulated genes. (Figure 3)

Figure 3 Differentially expressed genes between mice with SCI and healthy controls. (A) Heatmap of the top 50 upregulated and downregulated genes. (B) Volcano plot for DEGs between SCI tissues and healthy controls.

DEGs Enrichment and PPI Analysis

The biological functions and signaling pathways associated with SCI were obtained by GO and KEGG analysis. The results showed that the DEGs were enriched in the inflammatory response, regulation of cytokine production, neutrophil chemotaxis and degranulation, angiogenesis, and cell death. The KEGG results showed that the differential genes were highly correlated with the TNF signaling pathway and osteoclast differentiation (Figure 4A-B). Also, the PPI network correlations of the DEGs were demonstrated with a PPI enrichment P value <1.0e-16. The network consisted of 239 edges and 367 nodes with tight connections among the nodes (Figure 4C).

Figure 4 Functional enrichment analysis and PPI network of DEGs. (A) GO enrichment analysis. (B) The network of GO enrichment analysis. (C) Protein-protein interaction (PPI) network.

Co-Expression Network Construction and Hub Module Identification

A co-expression network was constructed using 20814 genes. The samples were processed to remove outliers and missing values, and a soft threshold = 3 (R²= 0.86, slope = −1.1) was selected to construct a scale-free network and build a co-expression matrix (Figure 5, Figure 6A). Four gene modules were obtained by dynamic hybrid shearing, and the correlation between modules and SCI clinical traits was shown using heat maps (Figure 6B and C). The turquoise module had the highest correlation (Cor= 0.74, P = 2e-06), in which the GS of the turquoise module was highly correlated with the module membership (MM) (Cor= 0.89, P < 1e-200) (Figure 6D). Therefore, the signature genes in the turquoise module were selected for subsequent analysis.

Figure 5 Determination of the soft threshold power in the weighted gene co-expression network analysis (WGCNA). (A) Left: analysis of the scale-free fit index for various soft threshold powers. Right: analysis of the mean connectivity for various soft threshold powers. (B) Histogram of connectivity distribution and checking the scale-free topology.

Figure 6 Construction of Weight Gene Co-Expression Network and Identification of the Key Module. (A) Clustering dendrogram of genes. In the figure, each branch represents a gene and on the bottom each color represents a co-expression module. (B) Gene significance in the modules. (C) Heatmap of the association between module genes and clinical traits. (D) Module membership and gene significance analyses of the turquoise module.

Identification of Hub Genes

A total of 478 candidate Hub genes were screened out from the turquoise module using the screening criteria of GS > 0.5 and MM > 0.8 and subsequently intersected with DEGs to obtain 404 intersecting genes (Figure 7). Then, the LASSO regression algorithm was used to obtain 4 hub genes: 2010004M13Rik, Cdkn1c, Nox4, and Gpr101.

Figure 7 LASSO Analysis and Validation of Hub Genes. (A) Venn diagram of intersecting genes between DEGs and the turquoise module. (B) Coefficients distribution trend of LASSO regression. (C) Distribution of hub genes in cross validation.

Hub Gene Expression Levels and Diagnostic Value Analysis

The expression levels of 4 hub genes in the SCI and control groups were analyzed using BOXPLOTS. We found that the expression levels and amounts of 2010004M13Rik (P < 0.001) and Nox4 (P < 0.001) were significantly higher in the spinal cord injury group than in the healthy control group, while the expression levels and amounts of Cdkn1c (P < 0.001) and Gpr101 (P < 0.001) were significantly higher in the healthy control group than in the spinal cord injury group (Figure 8A). Next, we established mouse spinal cord injury models, obtained corresponding normal and injured spinal cord tissues, extracted RNA for reverse transcription, and performed qRT-PCR (Figure 9); the results were consistent with those of the previous analysis. In the ROC curve analysis, the area under the curve (AUC) of the hub genes was used to assess the accuracy of spinal cord injury marker genes. The results showed that the AUC values of all four Hub genes were > 0.9 (Figure 8B), indicating that these are key genes with high diagnostic accuracy for spinal cord injury.

Figure 8 (A) Expression levels of hub genes. (B) Validation of hub genes in the diagnostic value (Higher the value of AUC, better the diagnostic value of the gene). “***” indicates that P-value < 0.001.

Figure 9 Expression levels of hub genes in qRT-PCR. (A) Cdkn1c is higher in the control group than in the SCI group. (B) Gpr101 is higher in the control group than in the SCI group. (C) Nox4 is higher in the SCI group than in the control group. “**” indicates that P-value < 0.01. “***” indicates that P-value < 0.001.

Analysis of Immune Cell Infiltration

We used the CIBERSORT algorithm to further investigate the differences in immune cell infiltration between the spinal cord injury group and healthy controls. The analysis of immune cell infiltration showed that eosinophils (P < 0.05) and GammaDelta T cells (P < 0.05) were infiltrated at significantly higher levels in the control group than in the spinal cord injury group, and M0 macrophages (P < 0.05), M2 macrophages (P < 0.05), T Cells CD4 Follicular (P < 0.05), and DC Actived (P < 0.05) infiltrated at significantly higher levels in spinal cord injury than in the control group. At the same time the relationship between individual immune cells is visualised (Figure 10).

Figure 10 Assessment of Immune Cell Infiltration. (A)Heatmap. (B) Relationship between immune cells. (C)Violin diagram of immune cell infiltration.

Relevance of Hub Genes to Immune Cells

Analysis of the correlation between four Hub genes and immune cells, the results show that in terms of immune cells, those significantly positively correlated with 2010004M13Rik were DC Actived, M0 Macrophage, T Cells CD4 Follicular and those significantly negatively correlated were GammaDelta T Cells, Monocyte, T Cells CD4 Memory and other cells (Figure 11A). Positively correlated immune cells with Cdkn1c are GammaDelta T Cells, T Cells CD4 Naive, Monocyte and negatively correlated immune cells with Cdkn1c are M0 Macrophage, M1 Macrophage (Figure 11B). For Gpr101, there was a positive correlation with GammaDelta T cells, T cell CD4 Naïve and M2 macrophages, and a negative correlation with M0 macrophages DC Actived and T cell CD4 Follicular (Figure 11C). In terms of NOX4, the positively correlated cells were M0 Macrophage, M1 Macrophage, DC Actived and the negatively correlated cells were T Cells CD4 Naive, Plasma Cells, GammaDelta T Cells (Figure 11D).

Figure 11 Lollipop diagram of relationship between immune cell infiltration and four hub genes. (A) The association between immune cell infiltration and 2010004M13Rik. (B) The association between immune cell infiltration and Cdkn1c. (C) The association between immune cell infiltration and Gpr101. (D) The association between immune cell infiltration and Nox4.

Discussion

In the present study, we statistically analyzed the DEGs between SCI and normal tissue in the GSE5296 and GSE47681 datasets to obtain DEGs. The exploration of the relevance of differential genetic, biological processes, and signaling pathways revealed that DEGs were mainly enriched in the inflammatory response, regulation of cytokine production, neutrophil chemotaxis and degranulation, angiogenesis, and cell death. Previous studies have shown that a prolonged immune inflammatory response can exacerbate the injury, cause neuronal death, impair axonal growth, and promote spinal cord cavity formation, severely affecting neurological recovery.^10,12,13 Immune cell infiltration is a major condition for the immune-inflammatory response,¹⁴ the chemotaxis and activation of immune cells cannot be distinguished from the mediation of cytokine, neutrophils mediated by C-X-C motif chemokine ligand 1 (CXCL1) are the first inflammatory immune cells to infiltrate the injury site.^2,15 Recent studies have shown that neutrophils are a major source of atypical growth factor oncogenic regulatory protein (Ocm), which promotes the regeneration of neurological axons and the recovery of neurological function.⁵ The activated neutrophils can also release neutrophil extracellular traps (NETs) into the extracellular environment by degranulation, which act synergistically with reactive oxygen species (ROS) to promote an inflammatory response and exacerbate the damage and affect neurological recovery.^16,17 The spinal cord resident immune cells and the first infiltrating neutrophils are stimulated to express inflammatory chemokines (IL-1, TNF, and colony-stimulating factor) and attract peripheral immune cell infiltration. Moreover, these cytokines can also activate immune cells and influence the immune-inflammatory response,^18,19 for example, T helper cells cytokine-1 (Th-1), IFN-γ, and TNF-α stimulate microglia/macrophages to convert to pro-inflammatory M1 type, secrete pro-inflammatory factors to exacerbate inflammatory responses,^20–22 and express major histocompatibility complex II (MHC II), which activates and regulates T lymphocytes to mediate innate and adaptive immune responses.²¹ In addition, the results of the KEGG signaling pathway analysis showed that the DEGs were mainly enriched in the TNF signaling pathway and osteoclast differentiation. TNF signaling pathway is a series of signaling pathways activated by the TNF protein superfamily acting on cells. For example, the function of TNF-a is achieved mainly by binding to two TNF receptors (TNFR1 and TNFR2), which induce NF-kb and MAPK activation via TNFR1-TRADD-TRAF to promote the inflammatory response.^23,24 Osteoclasts, mainly derived from the monocyte/macrophage hematopoietic stem cell lineage,²⁵ are a major cause of bone loss after SCI. While acting as a bone resorber, it also releases inflammatory cytokines, such as interleukin 6 (IL-6), which exacerbate the inflammatory response.²⁶ The DEGs enrichment results and the current studies have demonstrated that these differential genes play a key role in the biological processes related to SCI and recovery of neurological function.

To further explore the Hub genes associated with spinal cord injury and the promotion of neurological recovery, we crossed the highly correlated genes screened by WGCNA with differential genes to obtain crossover genes that are both differential and correlated. Subsequently, 4 hub genes, 2010004M13Rik, Cdkn1c, Nox4, and Gpr101, were obtained utilizing the LASSO algorithm.

Detailed studies on the role of 2010004M13Rik have not yet been reported. Intriguingly, NADPH oxidase 4 (Nox4) is the main isoform of the NADPH oxidase (Nox) family. As one of the primary sources of Reactive oxygen (ROS), Nox4 can directly produce hydrogen peroxide (H₂O₂), which has a direct toxic effect on cells²⁷ and also induces pro-inflammatory cytokine release through the activation of the JAK2-STAT3 pathway and NLRP3 inflammasome to produce the inflammatory response.²⁸ Yu et al promoted the release of inflammatory factors (IL-1b, IL-6, IL-8, IL-17, TNF-a, and MCP-1) by upregulating SAA1 protein expression, elevating NADPH oxidase activity, and promoting inflammation in vascular smooth muscle cells (VSMCs) on the activation of the p38-MAPK/NF-kB pathway, leading to increased inflammation.²⁹ The immune-inflammatory response is a major form of secondary injury following spinal cord injury, and Nox4 expression is significantly upregulated in mice.³⁰ Bermudez et al found that Nox4 upregulation may influence the change from a phase-equilibrium M1/M2 microglia ratio to a high M1 ratio, contributing to inflammation and exacerbating the severity of spinal cord injury.³¹ Based on this phenomenon, Nox4 has been extensively studied as a target for the treatment of inflammatory responses. Some studies have shown that targeting Nox4 by miR-99a reduces the ROS levels and decreases the expression levels of inflammatory markers, such as TNF-α, IL-6, and IL-1β, effectively inhibiting the inflammatory response in the spinal cord injury models and promoting the recovery of their neurological functions.³² Cell cycle protein-dependent kinase inhibitor 1c (Cdkn1c) is a cell cycle inhibitor with a unique domain structure consisting of an N-terminal CDK inhibitory domain (CdK), a proline-alanine repeat (PAPA) structural domain, and a C-terminal proliferating cell nuclear antigen (PCNA) binding domain.^33–35 In recent years, functional mutations in Cdkn1c have been associated with Beckwith–Wiedemann syndrome (BWS), IMAGe syndrome, and Silver–Russell syndrome.^36–38 The Cdkn1c mutation (loss-of-function variant), which causes loss of cell cycle inhibition, is present in 50% of familial cases of BWS and exhibits signs similar to those of BWS patients in mouse models with Cdkn1c knockdown;^36,39 unlike the mutation in BWS, the latter two congenital disorders show gain-of-function mutations in Cdkn1c, resulting in growth retardation.^37,38,40 Cdkn1c inhibits the cell cycle, leading to blocked cell growth and division,⁴¹ and this inhibitory effect on the growth of inflammatory cells in spinal cord injury warrants further study. In addition, Cdkn1c was shown to play a role in promoting cell differentiation. Accumulating evidence suggested that Cdkn1c is consistently expressed at higher levels in differentiated embryonic stem (ES) cells than in these cells when undifferentiated.⁴² Moreover, Mademtzoglou et al demonstrated that Cdkn1c balances the proliferation and differentiation of skeletal muscle stem cells (MuSCs) and promotes muscle regeneration.⁴³ Due to the specific characteristics, neuronal cells cannot complete division to proliferate, making neurological function recovery challenging after spinal cord injury. A large number of neuronal cells are derived from neural stem cells, and the differentiation of skeletal muscle stem cells by Cdkn1c is valuable for studying the role of stem cells differentiating into neurons to promote neural regeneration. G protein-coupled receptor 101 (Gpr101) is an orphan G protein-coupled receptor that is detected on neutrophils, monocytes, and macrophages in both humans and mice.⁴⁴ Gpr101 was first described in the study of X-LAG,⁴⁵ and the majority of studies on Gpr101 are focused on disorders related to growth hormone and lactogen secretion.^46,47 However, novel ideas have emerged in recent years. Magdalena B. Flak et al discovered that Gpr101 inhibits the development of inflammation by binding to N-3 docosapentaenoic acid-derived solving in D5 (RvD5n-3 DPA), enhancing phagocytosis by neutrophils and macrophages while limiting the infiltration of various types of immune cells into inflammatory sites.⁴⁸ Nonetheless, previous studies have conducted a comparative analysis of macrophages collected from Mac^gpr101KO mice and Gpr101^fl/fl mice and found that Gpr101 deficiency leads to an upregulation of the macrophage phenotypic markers IL-10 receptor (IL-10r), MCH II, and IL-23r, as well as pro-inflammatory eicosanoids, such as PGE2 and PGD2, confirming the phenotypic shift of Gpr101-deficient macrophages to the M1 pro-inflammatory phenotype infiltration in spinal cord injury vs normal spinal cord tissue. The results showed increased infiltration.⁴⁹ Together, our findings provide strong support for revealing the complexity of spinal cord injury repair mechanisms, but also provides more precise biomarkers (Cdkn1c, Nox4, Gpr101, etc). for clinical diagnosis and treatment, which promotes the development of precision medicine.

Since immune cell infiltration is crucial to the development of SCI and the recovery of neurological function, we used the CIBERSORT algorithm to calculate the difference in immune cell i DC Actived, M0 macrophages, M2 macrophages, and T Cells CD4 Follicular in the injured tissue.

Dendritic cells (DCs) are a class of antigen-presenting cells that play a critical role in linking non-specific and specific immunity.⁵⁰ Pan et al demonstrated that the maturation of DCs is impaired after spinal cord injury, and the cells are unable to function as antigen-presenting cells compared to the DCs phenotype of patients with spinal cord injury to that of healthy controls.⁵¹ The activation and maturation of DCs to promote neurological recovery after spinal cord injury has been well-documented. Wang et al found that injection of mature DCs stimulated the release of neurotrophic factors to promote neuronal differentiation of neural stem cells.⁵² The main role of T Cells CD4 Follicular is to promote the differentiation of B cells into antibody-secreting plasma cells and memory B cells. Furthermore, B cells are activated by infiltrating T Cells CD4 Follicular after spinal cord injury, and the resulting antibodies inhibit neurological recovery after spinal cord injury.⁵³ As described above, the role of macrophages has both advantages and disadvantages for the development of spinal cord injury and subsequent neurological recovery, which need to be further investigated with respect to the ratio of subpopulations, their duration of presence, and interactions. Together, our findings

In the present study, we identified 4 hub genes in SCI and verified the accuracy of the results by qRT-PCR, clarified the infiltration of immune cells, provided potential immunotherapeutic targets for inflammation control, neuroprotection and recovery after spinal cord injury, and offered a new therapeutic approach to alleviate poor outcomes after spinal cord injury. However, the roles of the characteristic genes and immune cells in spinal cord injury need to be explored by further experiments.