Introduction

Spinal cord injury (SCI), a devastating disease of the central nervous system (CNS), is characterized by severe functional loss that leads to long-term disability and associated complications1,2. The complex pathophysiology after SCI significantly impedes functional recovery3. Current treatment strategies, including acute surgical decompression, pharmacological interventions, and cell therapy coupled with neurorehabilitation, prevent secondary injury and decelerate disease progression4,5. However, these approaches do not promote neurological regeneration, resulting in a lack of curative therapies for SCI6,7,8.

Recent studies have revealed significant restoration of neural and motor function in neonatal mice after SCI. Notably, there are obvious differences in the recovery of neonatal mice and adult mice after SCI, which may be related to the differences in the expression of genes and related substances between them. The expression of neurotransmitters in glutamatergic interneurons in neonates is excitatory phenotype, but in adults it is inhibitory. Related single-cell sequencing has also found that some cells with ependymal gene expression profiles are only significantly enriched in neonatal injured spinal cord9,10. SCI in neonates exhibits significant upregulation of genes associated with axon growth, cell proliferation, and myelination, this spontaneous functional recovery observed in neonates is reduced in adults, which may be attributed to the expression of genes in most adult neurons that restrict axon growth and limit the energy available for regeneration11,12,13. In addition to neurons, glial cells also play an important role in the success or failure to SCI recovery in neonatal or adult mice, respectively14. Further studies are necessary to elucidate the contribution of reactive glial cells to the different regenerative capacity observed between neonatal and adult SCI.

Microglia, resident macrophages in the CNS, play a pivotal role in SCI recovery, particularly in neonates15,16,17. The results of one study have suggested that microglia in neonatal mice play an important role in coordinating SCI recovery. Transient activation of neonatal microglia promotes scar-free healing, and transplantation of neonatal microglia can significantly improve spinal cord lesions in adults18. In adult mice, although many studies have affirmed the therapeutic effect of microglia on SCI, the recovery of SCI in adults cannot achieve significant results due to the differentiation of microglia into the classically activated microglia (M1) and its pro-inflammatory effect19. Therefore, it may be worth studying the differences in microglial transcriptome genes between neonatal mice and adult mice after SCI. This approach will provide valuable insights into the differential microglial response between neonates and adults, revealing novel therapeutic targets for SCI.

We used two publicly available high-throughput sequencing microglial datasets from the Gene Expression Omnibus (GEO) database (https://www.ncbi.nlm.nih.gov/geo/)20. These datasets included GSE150871 (neonatal mice) and GSE113566 (adult mice). Differentially expressed genes (DEGs) were identified via selection of samples at the corresponding 3 days post-injury (dpi) time point within each dataset. Subsequently, DEGs that exhibited discordant expression patterns between the two groups were isolated. Functional enrichment analysis of the screened DEGs was performed using Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways to elucidate relevant biological functions and enriched pathways. Next, we constructed a protein-protein interaction (PPI) network using Cytoscape software to identify hub genes. Furthermore, we constructed a transcription factor (TF)-hub gene-microRNA network to explore the regulatory mechanisms underlying these hub genes. A drug profile database was used to identify potential therapeutic drugs and compounds that target the identified hub genes. Finally, quantitative real-time polymerase chain reaction (qRT-PCR) was used to validate the expression levels of hub genes in neonatal and adult mouse SCI models. Figure 1 shows the schematics of the experimental procedures and validation strategies.

Fig. 1
figure 1

The schematics of the experimental procedures and validation strategies.

Full size image

Materials and methods

Mice

In this experiment, animals were 7-week-old female C57BL/6 mice and 2-day-old female C57BL/6 mice (Jinan Pengyue Biotechnology Co., Ltd., Jinan, China). Mice were housed under controlled conditions, including a 12-hour light/dark cycle, a temperature of 23 ± 2 °C, a relative humidity of 50 ± 5%, and free access to food and water at all times. All animal procedures, including euthanasia, were performed in accordance with the animal care regulations of Qingdao University. All animal procedures were conducted in accordance with the protocol approved by the Laboratory Animal Welfare Ethics Committee of Qingdao University (NO.202403C572020240620116)and the guidelines for the use of animals in the study according to the Animal Research: In Vivo Experiment Report (ARRIVE). All surgeries are performed under anesthesia and every effort is made to reduce the pain of the animal.

Data acquisition

Microglial transcriptomic data were obtained from the GEO database, including two datasets (GSE150871 and GSE113566) generated using the GPL17021 (Illumina HiSeq 2500, Mus musculus) platform. GSE150871 comprises 5 samples from neonatal mouse microglia (postnatal day 2) at three time points post-injury (0, 3, and 5 dpi). We selected samples from 0 (control) and 3 dpi for further analysis. Dataset GSE113566 investigated transcriptional changes in microglia of adult mice (6–8 weeks old) after SCI at 3, 7, and 14 dpi. Samples from 0 (control) and 3 dpi were also selected for subsequent analyses.

DEG analysis

Microglial gene expression profiles between 3 dpi and control samples in GSE150871 were compared using the Seurat (version 5.0.1, https://github.com/satijalab/seurat) and harmony packages within R statistical software (version 4.3.2)21. This analysis yielded an expression matrix of microglia-related genes, followed by DEG analysis. DEGs were defined as genes with a | log2(fold change) | > 1 and an adjusted p-value < 0.05. For GSE113566, the DESeq2 package in R was used to identify DEGs between 3 dpi and control samples. DEGs were defined using similar criteria (| log2(fold change) |> 1 and adjusted p-value < 0.05). Volcano plots for DEGs were generated using the online Chiplot platform (https://www.chiplot.online/Volcano_plot.html)22. Heatmaps depicting DEG expression patterns were constructed using the ComplexHeatmap package (version 2.13.1)23.

Comparison of DEGs between GSE150871 and GSE113566

To examine changes in gene expression profiles between the two datasets, we performed Venn diagram analysis using a bioinformatics platform designed for data analysis and visualization (https://www.bioinformatics.com.cn/en?p=5)24. We calculated the intersection of upregulated and downregulated DEGs in both datasets. Subsequently, these screened genes (named Gene set 1–4) were subjected to GO and KEGG pathway enrichment analyses.

Functional enrichment analysis of identified gene sets

GO enrichment analysis categorized genes based on their biological processes (BPs), cellular components (CCs), and molecular functions (MFs). KEGG enrichment analysis focused on identifying enriched pathways within the dataset. The Metascape database (https://metascape.org/) was used for GO and KEGG enrichment analyses25. Stringent criteria were applied, with a minimum overlap of 3 genes and a minimum enrichment score of 1.5. Pathways with p-values < 0.01 were considered significantly enriched. Metascape automatically clustered the enriched terms, and the term with the highest ˗ log(p-value) within each category was selected as the representative. The top 20 enriched GO terms and KEGG pathways were visualized using a combination of bubble diagrams, string diagrams, and Sankey diagrams generated with a bioinformatics platform (https://www.bioinformatics.com.cn/en).

PPI network

To identify interacting genes and potential hub genes, we constructed a PPI network using the Search Tool for the Retrieval of Interacting Genes (STRING, https://cn.string-db.org/) database26. Interactions with a combined score > 0.4 were considered significant. The resulting PPI network data were retrieved and visualized using Cytoscape software (version 3.8.0)27. Three centrality measures, including betweenness, closeness, and degree, were used for hub gene selection. Hub genes were selected based on the cumulative score obtained from all three centrality measures. The top 10–20 genes with the highest total scores were identified as hub genes through the Cytohubba plugin (http://apps.cytoscape.org/apps/cytohubba) in Cytoscape software.

Construction of TFs-hub Gene-miRNAs network and identification of drugs and compounds

We used the miRwalk database (http://miRwalk.umm.uni-heidelberg.de/) to identify potential miRNAs that regulate hub genes28. The analysis was restricted to miRNAs predicted to target the 3’ untranslated region of hub genes with a miRwalk score of ≥ 0.95. To ensure accuracy, miRNA interactions were validated using the miRTarBase database. The NetworkAnalyst 3.0 platform (https://www.networkanalyst.ca/) was used to predict TFs regulating hub genes29. TF-gene interactions were extracted from ENCODE ChIP-seq data using the BETA Minus algorithm. Only TF-gene interactions with a peak intensity signal < 500 and predicted regulatory potential score < 1 were used30,31. Finally, we constructed and visualized the TF-hub gene-miRNA regulatory network using Cytoscape software.

The Enrichment platform (https://maayanlab.cloud/Enrichr), linked to the Drug Signatures database, was used to predict drugs and compounds targeting the hub genes32,33. The top 10 results were ranked by their p-values and visualized as a bar chart.

SCI model

Female C57BL/6 mice aged 7 weeks and female neonatal C57BL/6 mice aged 2 days were used to establish SCI models17,34. After an intraperitoneal injection of 80 mg/kg pentobarbital for anesthesia, the surgical area was shaved and disinfected with 75% ethanol; a midline incision was made in the thoracic vertebra. Laminectomy was performed at the T10 vertebral level to expose the spinal cord. Tips of 0.1-mm forceps were carefully inserted on either side of the spinal cord to include the full width of the spinal cord. Then, the spinal cord was fully crushed for 2 s, while the dura matter was kept intact. The muscles and skin layers were sutured, and the mice were placed on a heating pad until full recovery. Manual bladder emptying was performed twice daily.

RNA extraction and qRT-PCR analysis

Spinal cord tissues were harvested from mice at 3 days post-SCI for RNA analysis. Mice were anesthetized with an overdose of pentobarbital intraperitoneally and then perfused with pre-chilled PBS through the heart until the mice’s liver and effluent blood become transparent, mice died painlessly. Quickly dissect a 1 cm piece of spinal cord (centered on the T10 laminectomy site) and place in liquid nitrogen for later use. Total RNA was extracted using TRIzol® reagent (Invitrogen, Carlsbad, CA, USA). Reverse transcription of RNA to complementary DNA (cDNA) was performed using the All-in-One First-Strand Synthesis Master Mix Kit (share-Bio, Shanghai, China). Next, qRT-PCR was performed using ChamQ Universal SYBR qPCR Master Mix (Vazyme, Nanjing, China). Glyceraldehyde 3-phosphate dehydrogenase (GAPDH) was used as an internal control for the normalization of expression levels, and the 2-ΔΔCT method was used to determine the relative expression levels of target genes. Primer sequences used for qRT-PCR are presented in Table 1.

Table 1 Primers of qRT-PCR used in this study.
Full size table

Statistical analysis

GraphPad Prism software (version 8.0.2; GraphPad Software, San Diego, CA, USA) was used for statistical analysis of qRT-PCR data. Student’s t-test was used to compare gene expression levels between groups. Data are presented as means ± standard errors of the mean. P-values < 0.05 were considered statistically significant.

Results

Identification of DEGs after SCI

Two publicly available transcriptomic datasets (GSE113566 and GSE150871) were retrieved from the GEO database to identify DEGs after SCI in neonatal and adult mice. For GSE113566, transcriptome data from mice (6–8 weeks old) at 3 days post-SCI were analyzed. In total, 1,281 DEGs, including 856 upregulated and 425 downregulated genes, were identified with an adjusted p-value < 0.05 and | log2(fold change) | > 1 (Fig. 2A). Similarly, for GSE150871, transcriptome data from neonatal mice (postnatal day 2) at 3 days post-SCI were analyzed using R software packages. In total, 1,458 DEGs, including 1,010 upregulated and 448 downregulated genes, were identified with adjusted p-value < 0.05 and | log2(fold change) | > 1 (Fig. 2B). Heatmaps were generated to visualize the expression patterns of DEGs in each dataset (Fig. 2C and D).

Fig. 2
figure 2

Analysis of DEGs after SCI in GSE113566 and GSE150871 datasets. (A) Volcano plot depicting the distribution of 1,281 DEGs in the GSE113566 dataset. Pink points represent upregulated genes, whereas green points indicate downregulated genes. (B) Volcano plot showing the distribution of 1,458 DEGs in the GSE150871 dataset. Yellow points represent upregulated genes, whereas blue points indicate downregulated genes. (C) Heatmap showing the distribution of DEGs in the GSE113566 dataset. The color gradient from pink to green indicates the relative expression level, with pink and green indicating upregulated and downregulated expression, respectively. (D) Heatmap illustrating the distribution of DEGs in the GSE150871 dataset. The color gradient from yellow to blue represents the relative expression levels, with yellow and blue indicating upregulated and downregulated expression, respectively.

Full size image

Comparative analysis of DEGs after SCI between neonatal and adult mice

To examine the DEG changes after SCI between neonatal and adult mice, we used a Venn diagram to visualize the overlap between DEGs identified in datasets GSE113566 and GSE150871. We compared upregulated DEGs in neonates with upregulated and downregulated DEGs in adults. This analysis revealed 196 DEGs upregulated in both groups, 11 DEGs upregulated in neonates but downregulated in adults (Gene set 1), and 803 DEGs upregulated exclusively in neonates but not differentially expressed in adults (Gene set 3) (Fig. 3A, C, and E). Similarly, we compared downregulated DEGs in neonates with those in adults. We identified 48 DEGs downregulated in both groups, 38 DEGs downregulated in neonates but upregulated in adults (Gene set 2), and 362 DEGs downregulated exclusively in neonates but not differentially expressed in adults (Gene set 4) (Fig. 3B, D, and F). Consequently, the above four Gene sets we considered worthy of further analysis.

Fig. 3
figure 3

Venn diagram depicting the overlap of DEGs identified in neonatal and adult mice. Gene set 1 (11 DEGs) is upregulated in neonates but downregulated in adult mice. Gene set 2 (38 DEGs) is downregulated in neonates but upregulated in adult mice. Gene set 3 (803 DEGs) is upregulated in neonates but not differentially expressed in adult mice. Gene set 4 (362 DEGs) is downregulated in neonates but not differentially expressed in adult mice.

Full size image

Functional enrichment analysis of gene sets

To elucidate the potential biological functions and underlying mechanisms associated with the four identified Gene sets, we used GO and KEGG enrichment analyses. Due to the limited number of DEGs in Gene set 1, GO or KEGG analyses were not applicable. Therefore, we summarized the specific information and reported functions of these genes based on a comprehensive literature review (Table 2). For Gene set 2, GO enrichment analysis revealed significant enrichment of BPs related to DNA metabolism, replication, mitotic cell cycle progression, regulation of cell cycle phase transitions, homologous recombination, and base-excision/postreplication repair (Fig. 4A–B; Table 3). Gene set 3, including 803 DEGs, exhibited significant enrichment of BPs associated with immune regulation and response (cell activation, leukocyte migration), interleukin (IL)-1 and tumor necrosis factor (TNF) superfamily cytokine production, mitogen-activated protein kinase cascade signaling, apoptosis, neuroinflammation, cell-cell adhesion, and endocytosis (Fig. 4C–D; Table 3). Gene set 4 (362 DEGs) exhibited enrichment of BPs involved in DNA metabolism, double-strand break repair via homologous recombination, cell cycle regulation, cell morphogenesis, and cell adhesion (Figs. 4E–F; Table 3). Detailed enrichment results for molecular function and cellular component categories are provided in Supplementary file. KEGG pathway analysis of Gene set 2 (38 DEGs) identified significant enrichment of pathways associated with DNA replication, homologous recombination, cell cycle, base excision repair, and p53 signaling (Fig. 5A; Table 4). Furthermore, Gene set 3 exhibited enrichment of key immune signaling pathways, including TNF signaling, B cell receptor signaling, phosphoinositide 3-kinase/protein kinase B (PI3K-AKT) signaling, and peroxisome proliferator-activated receptors (PPAR) signaling pathways (Fig. 5B; Table 4). The main enrichment pathways in Gene set 4 included those involved in cell cycle regulation (Cell cycle pathway), DNA repair (Fanconi anemia pathway), DNA replication, microRNA regulation in cancer, adipocyte lipolysis regulation, and axon guidance (Fig. 5C; Table 4).

Table 2 Differentially expressed genes (DEGs) upregulated in neonates but downregulated in adult mice (Gene set 1).
Full size table
Fig. 4
figure 4

GO enrichment analysis. (A,C,E) Bubble diagrams depict GO term enrichment analysis. Horizontal and vertical axes represent fold enrichment and specific GO pathways, respectively. Bubble size corresponds to the number of enriched genes. The color gradient transitions from red to blue, denoting enrichment significance (high to low). (B,D,F) String diagrams depict the distribution of screened DEGs across various enriched GO functions. The left side of each plot indicates the DEGs, whereas the right side represents the different GO terms. Connecting bands visualize genes assigned to their corresponding GO terms. (A,B) Major enriched BPs associated with 38 DEGs in Gene set 2. (C,D) Major enriched BPs associated with 803 DEGs in Gene set 3. (E,F) Major enriched BPs associated with 362 DEGs in Gene set 4.

Full size image
Table 3 Gene ontology (GO) enrichment analysis results based on biological processes (Gene sets 2–4).
Full size table
Fig. 5
figure 5

KEGG pathway enrichment analysis. Bubble diagrams depict enriched KEGG pathways, KEGG pathway Image usage copyright was obtained through Kanehisa laboratories102,103,104,105. Horizontal and vertical axes represent the gene ratio and specific KEGG pathway names, respectively. Bubble size corresponds to the number of enriched genes. The color gradient transitions from red to blue indicate enrichment significance (high to low). Sankey diagrams illustrate the distribution of screened DEGs across various enriched KEGG pathways. The left side of the plot represents the DEGs, whereas the right side indicates the different KEGG terms. Connecting bands visualize genes assigned to their corresponding KEGG terms. (A) KEGG pathways associated with 38 DEGs in Gene set 2 (total of 5 pathways). (B) KEGG pathways associated with 803 DEGs in Gene set 3 (total of 20 pathways). (C) KEGG pathways associated with 362 DEGs in Gene set 4 (total of 6 pathways).

Full size image
Table 4 Kyoto encyclopedia of genes and genomes (KEGG) enrichment analysis results (Gene sets 2–4).
Full size table

PPI network analysis and hub gene identification

To elucidate the complex regulatory mechanisms underlying SCI, we used PPI network analysis to reveal functional relationships, characteristics, and biological pathways between DEGs. We constructed PPI networks in the STRING database using 38, 803, and 362 DEGs from Gene sets 2, 3, and 4, respectively. The PPI network for Gene set 2 (38 nodes, 326 edges) exhibited a high average local clustering coefficient (0.727) and significant PPI enrichment (p-value < 1.0e-16). Similarly, the PPI networks for Gene sets 3 (730 nodes, 3014 edges) and 4 (327 nodes, 397 edges) displayed high average local clustering coefficients (0.412 and 0.325, respectively) and significant PPI enrichment (p-value < 1.0e-16). These networks were visualized in Cytoscape (Fig. 6A, C, and E). Node size and color intensity were correlated with degree (number of interacting partners); larger, darker nodes represented higher-degree proteins. Subsequently, Cytoscape plugins using degree, closeness, and betweenness algorithms were used to identify hub genes within each network (Fig. 6B, D, and F). The top 10–20 genes with the highest total scores from these algorithms were designated as hub genes (Table 5). Finally, we compiled a list of 61 differentially expressed hub genes, including 50 from Gene sets 2–4 and 11 from Gene set 1 (Table 6). These 61 genes are considered key microglia-dependent genes with the most important roles in SCI repair.

Fig. 6
figure 6

PPI network construction and differentially expressed hub gene identification. (A,C,E) PPI networks constructed using the cytoNCA plugin. The size and color depth of each node correspond to its “degree” value (number of connections). (B,D,F) PPI networks constructed using the cytoHubba plugin, with three algorithms of betweenness, closeness, and degree. The color gradient from red to yellow indicates the scores assigned by each algorithm (high to low).

Full size image
Table 5 The top 10–20 genes with the highest total scores from the three algorithms (betweenness, closeness, and degree) identified as hub genes.
Full size table
Table 6 The main hub genes identified within the microglia network.
Full size table

Construction of TF-hub gene-miRNA network and drug prediction

To further investigate the potential regulatory mechanisms associated with SCI, we analyzed the interactions between TFs and miRNAs of the 61 differentially expressed hub genes. We uploaded these genes to the NetworkAnalyst platform to predict and establish a TF-hub gene network. In total, 43 hub genes interacted with TFs, forming a network with 76 nodes and 220 interactions (Fig. 7A). Next, we used miRwalk to construct a hub gene-miRNA network encompassing 53 nodes and 83 interactions (Fig. 7B). Subsequently, both networks were integrated into Cytoscape to generate a comprehensive TF-hub gene-miRNA network diagram (Fig. 7C). The construction of TF-hub-miRNA network diagram reveals the complex interaction between transcription factors and microRNAs in gene expression regulation, which is helpful for understanding the molecular mechanism of gene expression regulation in specific disease states. Transcription factors regulate the transcription rate of genes by binding to specific DNA sequences, while miRNAs regulate gene expression at the post-transcriptional level by degrading target gene mRNA or inhibiting its translation, and they are involved in the regulation of driving disease processes and may therefore be potential drug targets.

Fig. 7
figure 7

TF-hub gene-miRNA network, drug prediction and qRT-PCR validation of hub genes. (A) TF-hub gene network (red V-shapes: hub genes, purple ellipses: TFs). (B) Hub gene- miRNA network (red ellipses: hub genes, green ellipses: miRNAs). (C) Integrated TF-hub gene-miRNA network (red V-shapes: hub genes, purple ellipses: TFs, green triangles: miRNAs). (D) The top 10 candidate drugs identified by 61 hub genes. The color gradient from brighter to darker indicates the significance of each drug (high to low). (E) Genes from Gene sets 1–4 subjected to further validation using qRT-PCR. (F) qRT-PCR validation of 11 genes (from panel (E)) in neonatal mice after SCI. (G) qRT-PCR validation of 11 genes (from panel (E)) in adult mice after SCI. Data are expressed as means ± standard errors of the mean (SEMs) (n = 3). *P < 0.05, **P < 0.01, and ***P < 0.001.

Full size image

Furthermore, we investigated potential therapeutic interventions by identifying drugs and compounds that target these microglia-dependent DEGs. The top 10 candidate drugs included troglitazone, resveratrol, bisphenol A, indole-3-carbinol, PD98059, genistein, dasatinib, dexamethasone, phorbol-12-myristate-13-acetate, and methaneselenic acid (Fig. 7D), then we listed the details of the genes targeted by each drug(Table 7).

Table 7 Top 10 candidate drugs and their target microglia-dependent DEGs.
Full size table

qRT-PCR validation of hub genes

To validate the identified hub genes and assess their biological relevance, we performed qRT-PCR on a selection of 11 genes from the 61 microglia-dependent DEGs. These genes (Chek1, Rrm2, Lyve1, Mboat1, Clec4a3, Ccnd1, Cdk6, Zeb1, Igf1, Pparg, and Cd163) were selected by reviewing the comprehensive literature and considering their potential therapeutic role in SCI, with 2–3 genes from each Gene set (Fig. 7E)35,36,37,38,39,40,41,42. These genes have been implicated in promoting SCI recovery in mouse models. qRT-PCR analysis largely confirmed the expression patterns observed in the sequencing data. For example, Chek1 and RRM2 exhibited downregulation in neonatal mice but upregulation in adult mice after SCI. Clec4a3 expression was upregulated in neonates but downregulated in adults. Also consistent with the analysis results, Ccnd1 and Zeb1 were downregulated in neonatal mice and not differential expressed in adults after SCI (Fig. 7F-G).

Discussion

Despite advancements in modern medicine, a definitive cure for SCI remains elusive. This situation presents a significant challenge to global healthcare systems, primarily due to the currently limited understanding of cellular responses to SCI43. Neonatal mice exhibit remarkable regenerative capacity after SCI, and transplantation of microglia isolated from neonates promotes regeneration in adult mice. However, the mechanisms underlying this repair-promoting effect of neonatal microglia remain unclear18,44. Identification of transcriptomic differences between microglia from neonates and adults after SCI could reveal novel therapeutic targets for microglia-dependent treatment strategies. We utilized bioinformatics analysis to compare microglial transcriptomic changes in neonatal and adult mice at 3 days post-SCI. This analysis revealed DEGs, which were subsequently validated through in vivo experiments.

We compared DEGs in neonatal mice with DEGs in adult mice after SCI. Four Gene sets were selected for further analysis, including genes upregulated in neonates but downregulated in adults, genes downregulated in neonates but upregulated in adults, genes upregulated in neonates but not differentially expressed in adults, and genes downregulated in neonates but not differentially expressed in adults. GO and KEGG enrichment analyses revealed enrichment terms in these Gene sets primarily associated with DNA transcription and replication, cell differentiation and regeneration, post-injury immune response, and neuroinflammation.

Acute SCI elicits a robust immune response and neuroinflammation, key contributors to secondary injury. These processes begin within hours of the initial insult and persist for months within the lesion site45. Ischemia, oxidative stress, edema, and excitotoxic glutamate release initiate the cascade within a few hours after injury. Subsequently, a diverse array of inflammatory cells infiltrates the peri-lesional area, triggering substantial upregulation of pro-inflammatory cytokines, including IL-1β, TNF-α, and IL-646,47. At the same time, limited spontaneous repair processes initiate in neonatal and adult mice after primary or secondary injury-induced tissue and cell necrosis48. Consistent with these observations, our enrichment results revealed significant enrichment of DEGs in our identified Gene sets. These DEGs were functionally associated with DNA replication, homologous recombination, regulation of mitotic cell cycle phase transition, cell morphogenesis, base excision repair, cell proliferation and apoptosis, immune response regulation, leukocyte migration and activation, regulation of IL-1 and TNF production, and modulation of neuroinflammatory responses. Notably, several signaling pathways, such as p53, PI3K-AKT, PPAR, and AMP-activated protein kinase (AMPK) pathways, were highly associated with SCI repair through their roles in the regulation of immune inflammation, cell proliferation, apoptosis, and regeneration.

The p53 protein regulates diverse biological processes depending on the nature of cellular stress signals. Recent studies have characterized its involvement in numerous biological pathways, including autophagy, cell metabolism, ferroptosis, and reactive oxygen species production49. Moreover, p53 regulation can promote SCI recovery by regulating antioxidant defenses, apoptosis, and ferroptosis50,51. Ferroptosis, an iron-dependent form of programmed cell death, contributes to organ dysfunction and degeneration52,53. Increased free iron levels exacerbate SCI by inducing glutamate excitotoxicity in spinal motor neurons. Conversely, iron chelation demonstrates therapeutic potential in alleviating SCI54. Recent studies have implicated p53 in ferroptosis regulation. Notably, inhibition of dihydroorotate dehydrogenase (DHODH) can protect spinal cord neurons from ferroptosis by modulating p53 activity, thereby suppressing arachidonic acid 15-lipoxygenase (ALOX15), and promoting neuronal survival and SCI recovery52. Additionally, two other signaling pathways, PI3K/AKT and PPAR-γ, are significantly upregulated in neonatal mice and contribute to various cellular functions, including cell proliferation, survival, and inflammation55,56. Activation of the PI3K/AKT/mTOR pathway with resveratrol has demonstrated significant improvement in neurological function through reductions of cerebral infarction volume and neuronal damage57. PPAR-γ, a ligand-activated nuclear receptor, regulates glucose and lipid metabolism, endothelial function, and inflammation58,59,60. PPAR agonists have demonstrated significant effectiveness in reducing neuronal loss and improving motor function recovery after SCI61,62,63,64. The AMPK pathway, a key regulator of cellular metabolism and inflammation, plays a crucial role in SCI. AMPK regulates NLRP3 inflammasome assembly and activation, influencing inflammation and oxidative stress in various diseases65,66. Modulation of the AMPK signaling pathway can reduce oxidative stress, regulate spinal cord and neuronal glucose metabolism, and promote functional recovery post-SCI67. Therefore, genes associated with these pathways may represent promising therapeutic targets for SCI. For example, CCND1 and Cdk6, genes involved in p53 and PI3K/AKT pathways, warrant further studies. CCND1, a critical cyclin protein, exhibits increased expression 24 h post-injury. It activates CDK4/6, triggering a cascade of events that promote neuronal apoptosis and glial cell reactive proliferation, thereby contributing to the secondary injury process in SCI6,68. Notably, the CDK inhibitor flavopiridol demonstrates potential in restoring cell cycle control. Delayed administration of this drug significantly reduces the lesion volume, preserves white matter integrity, improves motor function recovery after SCI, and mitigates secondary tissue damage69,70.

TFs and miRNAs can regulate gene expression at the transcriptional level. Our TF analysis identified Myc as a potential key player in SCI. c-Myc, a member of the Myc oncogene family, critically regulates cell proliferation and differentiation. Several studies of traumatic SCI have demonstrated that c-Myc knockdown reduces oligodendrocyte precursor cell proliferation and promotes their differentiation and myelination capacity in vitro71. Furthermore, our network analysis revealed several miRNAs with important regulatory roles in SCI. For example, miR-329-3p targets Pik3r1, a gene downregulated in neonatal mice after SCI. Several studies have revealed that miR-329-3p inhibition can reverse lipopolysaccharide-induced neuronal apoptosis, inflammation, and oxidative stress. Conversely, miR-329-3p overexpression attenuates the neuroprotective effects of mesenchymal stem cell-derived exosomes in injured spinal cords72. The upregulated DEG Kit in neonatal mice interacts with miR-222-3p. Notably, miR-222-3p upregulation can inhibit neuronal apoptosis and promote motor function recovery after SCI73. Furthermore, miR-124-3p, the most enriched miRNA in neuron-derived exosomes, interacts with three hub genes: Anxa5, Cdk6, and Chek1. In vivo and in vitro studies have demonstrated that miR-124-3p suppresses M1 microglia and A1 astrocyte activation by inhibiting myosin heavy chain 9 (MYH9) activity, significantly promoting functional recovery after SCI74. Additionally, miR-124-3p promotes M2 macrophage polarization, reducing tissue injury in spinal cord ischemia-reperfusion injury and alleviating neurogenic pain in SCI mice75,76. These findings represent only a fraction of the miRNAs and TFs identified in the SCI repair network. Further studies are needed to explore novel therapeutic targets in SCI.

Resveratrol, a natural polyphenol with antioxidant, anti-inflammatory, antibacterial, and anti-neurodegenerative properties, has emerged as a potential therapeutic candidate; it can promote functional recovery after SCI by modulating gut microbiota and inhibiting microglial activation77,78,79,80,81. Extracellular signal-regulated kinase (ERK) inhibition plays a crucial role in the neuroprotective effects of various drugs82. PD-98,059, a predicted ERK inhibitor, could potentially alleviate neuropathic pain by enhancing ERK inhibition83. Genistein, a natural soy-derived phytoestrogen with anti-inflammatory, antioxidant, and neuroprotective effects, may reduce lipid peroxidation, oxidative damage, and inflammatory response after SCI84,85,86,87. Dexamethasone, a commonly used glucocorticoid, has a well-established role in SCI treatment. However, the efficacies of troglitazone, bisphenol A, indole-3-carbinol, dasatinib, and methane-selenic acid require further investigation.

We identified 11 genes potentially influencing SCI recovery from a pool of 61 DEGs for further verification using qRT-PCR. These genes were Chek1, RRM2, Lyve1, Mboat1, Clec4a3, Ccnd1, Cdk6, Zeb1, Igf1, Pparg, and Cd163. The expression profiles of Chek1, RRM2, Ccnd1, Mboat1, Clec4a3, Igf1, and Cd163 were consistent with the RNA sequencing results in neonatal mice after SCI. Notably, three genes (Chek1, RRM2, and Clec4a3) exhibited expression trends similar to those observed in adult mice. Chek1 is a well-established regulator of cell proliferation and apoptosis, the upregulation of Chek1 may indicate cellular senescence and coordinate DNA damage88,89. However, its specific contribution to SCI remains unclear. RRM2 plays a multifaceted role in regulating purine metabolism, the cell cycle, and DNA repair. The results of recent studies suggest its involvement in protein regulation and modification, as well as angiogenesis90. Notably, RRM2 exhibits a specific elevation in liver cancer and functions as an endogenous ferroptosis suppressor by stimulating glutathione synthesis via glutathione synthase91. This indicates the potential of serum RRM2 as a biomarker for ferroptosis diagnosis and treatment. Considering the emerging role of ferroptosis in various diseases, including SCI, RRM2 warrants further investigation as a potential therapeutic target for SCI. Clec4a3 is a C-type lectin gene within the Aplec cluster. One study has demonstrated that Aplec, which is widely expressed in microglia, upregulated under inflammatory stimuli, is involved in microglial signaling after SCI, however, its expression has no effect on glial cell activation. Moreover, it regulates motor neuron survival and lymphocyte infiltration, indicating its crucial role in modulating immune responses and neuronal cell survival92. Furthermore, Ccnd1, Mboat1, Igf1, and Cd163 emerged as potential therapeutic targets, particularly in neonatal mice. Cyclin D1 (Ccnd1), is upregulated in neurons and microglia after traumatic brain injury (TBI), and its downregulation has important implications for the improvement of neurodegeneration, the recovery of cognitive function, the reduction of lesion volume, and the activation of cortical microglia after TBI93. Mboat1, a phospholipid-modifying enzyme, is particularly intriguing. Ferroptosis, a form of cell death mediated by iron-dependent phospholipid peroxidation and plasma membrane rupture, has been implicated in SCI. Mboat1/2 function as ferroptosis inhibitors by remodeling phospholipids94,95. This suggests a potential neuroprotective role for Mboat1 in SCI, highlighting the need for further studies. Igf1 is a protein hormone with neurotrophic properties. It plays a critical role in CNS development, promoting neuronal differentiation and survival during normal development and aiding recovery after injury or disease96. Numerous studies have revealed its pleiotropic benefits in the CNS, including the promotion of SCI repair by inducing a shift in microglial polarization towards neuroprotective and anti-inflammatory phenotypes97,98. In the chronic stages of SCI, Igf1 therapy may reverse permanent hippocampal alterations and improve cognitive deficits99. It has also been shown to be significantly upregulated in microglia subclusters that are upregulated in ischemic brain tissue, and its high expression can transform microglia into a neuroprotective phenotype100. Similarly, the difference in the expression of this gene in microglia between neonatal mice and adult mice after SCI may be one of the reasons why neonatal mice recover better after injury than adult mice.

Our study has several limitations. First, considering the relatively small animal sample size, hub gene expression should be confirmed in a larger cohort. Second, we solely focused on data from the third day post-SCI. Examination of DEGs across multiple time points would provide a more comprehensive picture. Third, although we focused on genes that exhibit significant changes in neonates compared with adults, genes with less pronounced differential expression in neonates may also be relevant. Finally, the underlying mechanisms of the identified hub genes remain uncertain, emphasizing the need for future studies.

Conclusions

We compared microglial transcriptomic alterations between neonatal and adult mice at 3 days post-SCI. Bioinformatics analysis identified 61 microglia-associated genes implicated in SCI. Furthermore, we investigated transcription factors, miRNAs, and potential drug targets associated with these genes. Subsequently, we selected 11 hub genes, including Chek1, RRM2, Lyve1, Mboat1, Clec4a3, Ccnd1, Cdk6, Zeb1, Igf1, Pparg and Cd163, from the 61 genes for further validation. These findings provide novel insights into neonatal microglial therapy for SCI and will contribute to the identification of novel therapeutic targets for this condition.