Bioinformatics analyses of differentially expressed genes associated with spinal cord injury: A microarray-based analysis in a mouse model.

Gene spectrum analysis has shown that gene expression and signaling pathways change dramatically after spinal cord injury, which may affect the microenvironment of the damaged site. Microarray analysis provides a new opportunity for investigating diagnosis, treatment, and prognosis of spinal cord injury. However, differentially expressed genes are not consistent among studies, and many key genes and signaling pathways have not yet been accurately studied. GSE5296 was retrieved from the Gene Expression Omnibus DataSet. Differentially expressed genes were obtained using R/Bioconductor software (expression changed at least two-fold; P < 0.05). Database for Annotation, Visualization and Integrated Discovery was used for functional annotation of differentially expressed genes and Animal Transcription Factor Database for predicting potential transcription factors. The resulting transcription regulatory protein interaction network was mapped to screen representative genes and investigate their diagnostic and therapeutic value for disease. In total, this study identified 109 genes that were upregulated and 30 that were downregulated at 0.5, 4, and 24 hours, and 3, 7, and 28 days after spinal cord injury. The number of downregulated genes was smaller than the number of upregulated genes at each time point. Database for Annotation, Visualization and Integrated Discovery analysis found that many inflammation-related pathways were upregulated in injured spinal cord. Additionally, expression levels of these inflammation-related genes were maintained for at least 28 days. Moreover, 399 regulation modes and 77 nodes were shown in the protein-protein interaction network of upregulated differentially expressed genes. Among the 10 upregulated differentially expressed genes with the highest degrees of distribution, six genes were transcription factors. Among these transcription factors, ATF3 showed the greatest change. ATF3 was upregulated within 30 minutes, and its expression levels remained high at 28 days after spinal cord injury. These key genes screened by bioinformatics tools can be used as biological markers to diagnose diseases and provide a reference for identifying therapeutic targets.

Chinese Library Classification No. R447; R363; R741

Introduction

As part of the central nervous system, the spinal cord is crucial for conveying afferent and efferent impulses between the brain and somatic/visceral receptors, as well as executing reflexes. However, its vulnerability and limited capacity for regeneration and self-renewal makes recovery from mechanical trauma difficult. Indeed, severe spinal cord injury (SCI) often results in permanent functional impairment, such as motor/sensory dysfunction or bladder and rectal disturbances (Bastien et al., 2015; Jain et al., 2015; He and Jin, 2016). Thus, SCI can diminish a patient’s quality of life and cause a heavy burden for families (Qiu, 2009). For decades, doctors and patients have been searching for effective interventions and therapies for SCI that do not exhibit serious side effects. However, because of limitations in therapeutic applications, there are currently no available therapies for an effective recovery (Wyndaele and Wyndaele, 2006; Courtine et al., 2011; Yang et al., 2016). Schwab (2002) suggested four principal strategies for SCI repair: promoting regrowth of interrupted nerve fiber tracts, bridging spinal cord lesions, repairing damaged myelin, and restoring nerve-fiber impulse conductivity. Some promising therapeutic interventions, such as cell transplantation and metabolic interventions, have shown effectiveness in animal models (Zhang et al., 2009; Guerrero et al., 2012; Tsukahara et al., 2017; Nordestgaard et al., 2018). Nevertheless, these potential therapies must be further tested in animal models and validated in human clinical trials.

Clarifying the pathological and molecular changes after SCI is crucial because the endogenous mechanisms for repair and intervention may provide potent insight for therapy exploration. Injury can elicit inflammatory stimuli, which then influence the production of proinflammatory cytokines and other mediators (Peifer et al., 2006). Localized immune/inflammatory responses are important contributions to secondary tissue damage and functional deficits after SCI (Ghasemlou et al., 2010), and are also essential for cleaning tissue debris and remodeling and repair after injury. Previous studies (Carlson et al., 1998; Schnell et al., 1999; Hashimoto et al., 2003, 2005) have shown that the lesion phase can be divided into three stages. First, neutrophil infiltration and cell death dramatically increase (Liu et al., 1997). Second, macrophages/microglia accumulate and proliferate, which can result in harmful effects on surrounding tissues. Third, glial scars form, in which astrocytes play a harmful role in tract regeneration by surrounding the injured tissue and producing scar-associated compounds. The characteristics of inflammation and extent of glial scar formation represent the most distinct differences between the acute, sub-acute, and chronic phases of the SCI microenvironment. Hence, microenvironmental changes after SCI may affect prognosis (Nishimura et al., 2013), and must be considered for functional recovery at different stages.

Recently, gene profiling studies have suggested that gene expression and signaling pathways are substantially changed after SCI (Di Giovanni et al., 2003; Byrnes et al., 2011; Lee-Liu et al., 2014), and may influence the microenvironment of the injury site. Microarray analyses have provided novel perspectives on diagnosis, therapy, and prognosis prediction for SCI (Nesic et al., 2002; Xiao et al., 2005). However, differentially expressed genes (DEGs) are rarely consistent across different studies (Nesic et al., 2002; Di Giovanni et al., 2003; Xiao et al., 2005; Byrnes et al., 2011; Lee-Liu et al., 2014; Duran et al., 2017), and many important genes and pathways have not been thoroughly investigated. As DEGs may provide therapeutic targets, exact changes in the gene transcriptome and the underlying cellular and molecular mechanisms still need to be clarified, compared, and critically analyzed for different time points, especially the acute, sub-acute and chronic phases.

In the present study, we focused on alterations in gene expression in the thoracic spinal segment (T8) of the mouse spinal cord at six different time points (0.5, 4, and 24 hours, 3, 7, and 28 days) after SCI, because alterations at these time points are likely related to pathological changes associated with SCI. In this study, we compared gene expression patterns at different time points and analyzed them using bioinformatic methods, particularly from a protein-protein interaction perspective. The aim of our study was to identify potential genes/pathways and clarify the mechanisms underlying SCI at the molecular level.

Materials and Methods

Retrieval of microarray dataset

The microarray dataset, GSE5296, was downloaded from the Gene Expression Omnibus DataSet from the National Center for Biotechnology Information Database (http://www.ncbi.nlm.nih.gov/gds/). GPL1261 (Affymetrix Mouse Genome 430 2.0 Array) (Affymetrix Inc., Santa Clara Valley, CA, USA) was used as the platform. C57BL/6 mice were used as subjects. All animals were deeply anesthetized during surgery under isoflurane anesthesia. Moderate injury was delivered to each SCI mouse at T8. Mice were then sacrificed at the following time points: 0.5, 4, and 24 hours, and 3, 7, and 28 days (n = 3/per subgroup). Tissue samples from the impact site (0.4 cm in length) were collected after sacrifice. Each individual was pooled (4 mice in total), and 12 mice were prepared at each time point. A control group (n = 2/per subgroup) that underwent sham injury was included.

DEG profiling in tissue from the SCI impact site

Raw data were downloaded, and R/Bioconductor software (v3.4.1, R Foundation, Vienna, Austria) (Gregory Alvord et al., 2007; Barrett et al., 2013) was used to check the data quality. The software indicates whether the quality of each selected CEL was sufficient (). gcRobust Multi-array Average (gcRMA) analysis was used as the normalization algorithm for gene expression profiling, and the limma algorithm package was used for DEG identification (Gregory Alvord et al., 2007; Barrett et al., 2013). Spinal cord tissue from sham injury mice was used as the control. DEGs were defined as genes with at least 2-fold changed expression, with P-values < 0.05 between the control group and any subgroup at six time points. Hierarchical clustering analyses were also performed for DEGs.

Quality control of the CEL files in GSE5296 dataset based on R/Bioconductor program.

(A) Summarization of Quality control. (B) Plot of gcRobust Multi-array Average (gcRMA) algorithm. (C) RNA degradation plot. (D) Cluster Dengrogram. (E) Principal components plot.

Click here for additional data file.

Functional enrichment analysis of DEGs

The bioinformatics resource, Database for Annotation, Visualization and Integrated Discovery (DAVID) (Huang da et al., 2009) (https://david.ncifcrf.gov/) was used for functional annotation. The Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis and Gene Ontology (GO) analysis were used to identify genes of interest. Enriched terms for KEGG pathway and GO analysis were collected, including Gene Ontology Biological Process, Cellular Components, and Molecular Function (abbreviated as GO_BP, GO_CC, and GO_MF). Those with P-values < 0.05 were considered significantly enriched terms.

Prediction and analysis of transcription factors

Based on the Animal Transcription Factor Database (TFDB) (http://bioinfo.life.hust.edu.cn/AnimalTFDB/) (Zhang et al., 2015), DEGs that were transcription factors were identified. Interaction networks of DEGs were also constructed using Search Tool for the Retrieval of Interacting Genes/Proteins (STRING) (Franceschini et al., 2013) (v10.5) and Cytoscape software (Smoot et al., 2011) (v3.6.0; National Institute of General Medical Science, Bethesda, MD, USA) to identify protein-protein interactions. The score threshold was set to 0.04.

Statistical analysis

All statistical analyses were performed using GraphPad Prism software (v6; GraphPad Software, San Diego, CA, USA) and SPSS software (v20.0; IBM SPSS Inc., Chicago, IL, USA). P-values < 0.05 were considered statistically significant.

Results

Identification of DEGs from SCI site tissue

The gene expression dataset, GSE5296, was used to identify DEGs from a SCI murine model using a microarray representing approximately 34,000 genes. By comparing the control group and SCI subgroups at each time point, several DEGs were identified, with the numbers shown in . SCI led to upregulation of 2,460 genes in injured spinal cord at 7 days after injury, 109 of which were upregulated at each time point. Moreover, many of these upregulated genes were associated with inflammation. In total, 2010 genes were downregulated at the injury site at 7 days following SCI, 30 of which were downregulated at each time point. Interestingly, the number of downregulated genes was always smaller than the number of upregulated genes at each time point. Overall, we identified 109 upregulated genes and 30 downregulated genes with consistent directions of variation at each time point (). contains the complete list of DEGs.

DEGs at the injury site in mice.

The GSE5296 dataset contained a control group (sham injury) (n = 2/time point) and SCI group with six subgroups at different time points (including 0.5, 4, and 24 h, and 3, 7, and 28 d; n = 3/time point). DEGs between each SCI subgroup and control group were identified using the limma algorithm package of R/Bioconductor software. DEGS were defined as genes with at least 2-fold changed expression levels and P-values < 0.05. DEGs: Differentially expressed genes; SCI: spinal cord injury; h: hours; d: days.

Bi-cluster analysis of DEGs with consistent directions of variation at each time point following SCI compared with the control group.

Each row represents one DEG, and each column represents a tissue sample of T8 spinal cord from injured and sham-injured spinal cord. Control tissue was from sham-injured spinal cord. SCI tissue was from injured T8 spinal cord tissue. “SCI0.5h”, “SCI4h”, “SCI24h”, “SCI3d”, “SCI7d”, and “SCI28d” tissues were collected at 0.5, 4, and 24 hours, and 3, 7, and 28 days following SCI, respectively. The numbers 1, 2, and 3 represent biological repetitions. The heat map was constructed based on DEG expression values in each tissue sample. DEGs: Differentially expressed genes; SCI: spinal cord injury; h: hours; d: days.

Functional annotation of DEGs in SCI site tissue

Using the DAVID tool, enriched KEGG and GO terms were identified from the DEGs. shows the number of KEGG and GO terms that were enriched in each SCI subgroup. The results indicate a peak on day 7, which is similar to the tendency shown by the DEGs. For each SCI subgroup, the three KEGG pathways showing highest changes in enrichment scores are listed in . Similarly, those for the GO terms are shown in . Notably, KEGG analysis revealed that genes related to metabolic pathways were downregulated on day 7, with metabolic pathways also showing the largest number of enriched genes. This was also the finding on day 28. Some genes related to inflammatory processes were upregulated on either day 7 or 28 after SCI, with the majority of these genes upregulated at both time points.

Overview of KEGG and GO terms enriched at each time point of spinal cord injury compared with the control group.

All enrichment analyses were based on upregulated and downregulated genes with P-values < 0.05. KEGG: Kyoto Encyclopedia of Genes and Genomes; GO: Gene Ontology; h: hours; d: days.

DAVID analysis was performed for upregulated DEGs with the same direction of variation at each time point after SCI (see Tables and for details). Accordingly, many inflammation-related pathways were upregulated in injured spinal cord. One pathway showing the greatest change in gene expression profiling was the tumor necrosis factor signaling pathway. Inflammatory response was also top-listed when GO_BP terms were sorted by the number of enriched genes. For these inflammation-related genes, higher expression levels were maintained in the SCI animal model for at least 28 days. Thus, we next focused on these genes. A striking finding was that some C-C motif chemokine-ligand genes were markedly upregulated at each time point. Among these, CCL3 showed obvious upregulation within 30 minutes after SCI, with at least a 6-fold increase in expression at each time point until day 28. Additionally, the Timp3 gene, a tissue inhibitor of metalloproteinases, was found to be upregulated at each time point.

Protein-protein interactions

A protein-protein interaction network was constructed based on expression changes of DEGs with the same direction of variation at each time point after SCI (). The network contained 84 nodes and 410 pairs of connections. The 10 DEGs showing the highest degrees of change were: JUN, FOS, TRP53, PTGS2, CCL2, MYC, TLR2, STAT3, ATF3, and ICAM1, all of which primarily participate in the cytokine response.

Interaction networks of DEGs with the same variation tendency at each time point of SCI.

(A) Interaction network of DEGs with the same variation tendency at each time point of SCI compared with the control group. Circles represent DEGs (red for upregulated DEGs and green for downregulated DEGs) and lines represent interactions. (B) Interaction network of upregulated DEGs with the same variation tendency at each time point of SCI compared with the control group. Circles represent DEGs (red for upregulated DEGs and yellow for transcription factors). DEGs: Differentially expressed genes; SCI: spinal cord injury.

To investigate transcriptome changes after SCI, AnimalTFDB 2.0 was used for transcription factor prediction. To investigate transcriptional activity of transcription factors, we constructed an interaction network from DEGs that were upregulated at various time points. Several transcription factors were found among these DEGs. The network contained 399 regulation modes and 77 nodes involving 15 transcription factors (). Because DEGs with higher-degree distributions may play more important roles in SCI pathogenesis, we reviewed gene expression of transcription factors with the highest degrees. Overall, the 10 upregulated DEGs with the highest degrees were: JUN, FOS, TRP53, PTGS2, CCL2, MYC, TLR2, STAT3, ATF3 and ICAM1, six of which are transcription factors. Among these transcription factors, ATF3 exhibited the greatest change (almost 64-fold increase). Indeed, it was upregulated within 30 minutes after SCI, and its expression levels remained high even on day 28.

Discussion

Spinal cord injury can result in neurological dysfunction involving locomotor, sensory, and autonomic systems. Therefore, understanding its pathogenesis is crucial for further strategies aimed at prevention and therapies. In this study, we used gene expression profiling to identify and characterize a number of DEGs in the T8 segment of the spinal cord at different time points following injury. At each time point, the majority of these DEGs were upregulated, suggesting a greater amount of transcriptional activation.

Notably, KEGG analysis revealed that genes related to metabolic pathways were downregulated on day 7, when the largest number of enriched genes was observed, which was the same on day 28. Moreover, our results show that genes related to inflammatory processes were upregulated within 30 minutes, and continued to be upregulated on even day 28 after SCI, with a peak on day 7. Therefore, inflammation-related pathways were retrieved for Mus musculus (mouse) from the KEGG website (http://www.kegg.jp/). DAVID bioinformatics analysis revealed that almost half of these pathways were enriched for DEGs related to inflammation. Among these pathways, the tumor necrosis factor signaling pathway, rheumatoid arthritis, inflammatory bowel disease, toll-like receptor signaling pathway, phagosome, nuclear factor-kappa B signaling pathway, chemokine signaling pathway, cell adhesion molecules, MAPK signaling pathway, Staphylococcus aureus infection, and HTLV-I infection were identified.

Injury can elicit inflammatory stimuli, which then induce pathological accumulation and proliferation of inflammatory cells, a process that involves neutrophils, macrophages/microglia, and astrocytes (Liu et al., 1997; Carlson et al., 1998; Schnell et al., 1999; Hashimoto et al., 2003, 2005; Peifer et al., 2006; Ghasemlou et al., 2010; Nishimura et al., 2013). Naturally, various mediators are released into the microenvironment, including proinflammatory chemokines/cytokines and proteases (Peifer et al., 2006; Rice et al., 2007). Localized immune/inflammatory responses exacerbate the initial damage and contribute to secondary tissue damage after SCI, which can then result in functional deficits (Ghasemlou et al., 2010). However, these inflammatory responses are also essential for cleaning tissue debris and remodeling and repair after injury. Specifically, neutrophil invasion appears at approximately 6 hours after SCI, as evidenced by upregulation of inflammatory cytokines (Liu et al., 2002), such as interleukin-1α, interleukin-1β, interleukin-6, and tumor necrosis factor α. At 1 day after SCI, neutrophil number reaches peak levels at the injury site (Carlson et al., 1998; Schnell et al., 1999). Meanwhile, macrophages start to emerge at 3 days, reaching a peak at 7 days (Zhu et al., 2017). Macrophages may produce mediators that contribute to subsequent processes. Macrophage infiltration into the injury site at 7 days is best described as foam cells (Zhu et al., 2017), whose function largely depends on the surrounding microenvironment (Amit et al., 2016). Moreover, Guerrero et al. (2012) indicated that blockade of the interleukin-6 pathway in macrophages may promote regeneration of the spinal cord in mice. After the inflammatory reaction, the residual area is surrounded by glial scarring. Astrocytes begin to increase after SCI by day 2 and high numbers are maintained for at least 14 days (Baldwin et al., 1998), a situation that can produce potent inhibitors of neurite outgrowth (Wanner et al., 2013). Moreover, inducible nitric oxide synthetase is synthesized in glial scars, while regeneration of spinal cord tracts is inhibited by glial scar formation. Hence, the lesion phase can be divided into three stages (Hashimoto et al., 2005). First, neutrophils infiltrate and cell death dramatically increases (Liu et al., 1997). Subsequently, macrophages/microglia accumulate and proliferate, which may have dual effects on surrounding tissue. Glial scar formation comprises the final stage, and astrocytes hinder tract regeneration by surrounding the injured tissue and producing scar-associated compounds.

These findings led us to investigate the function of persistent dysregulated genes in the spinal cord after injury. In this study, we identified 109 upregulated genes and 30 downregulated genes that maintain the same variation tendency at each time point. Subsequently, we performed DAVID analysis separately for up- and downregulated DEGs. Not surprisingly, our results show that many inflammation-related pathways are upregulated in the injured spinal cord. The tumor necrosis factor signaling pathway exhibited the greatest change, with high expression levels of inflammation-related genes of this pathway lasting for at least 28 days in our SCI model mice. Strikingly, some C-C motif chemokine ligand genes were also markedly upregulated. Among these, the CCL3 gene showed the most obvious upregulation, beginning within 30 minutes after SCI and lasting for as long as 28 days (more than a 6-fold increase at each time point). This gene participates in the toll-like receptor signaling pathway, cytokine-cytokine receptor interaction, and the chemokine signaling pathway. CCL2 and CCL12 are involved in the tumor necrosis factor signaling pathway, nuclear factor-kappa B signaling pathway, and chemokine signaling pathway. Moreover, it has been reported that matrix metalloproteinases (MMPs) regulate inflammation after SCI (Zhang et al., 2011). Appropriately, our results show that MMP2 and MMP14 are upregulated at every time point. Additionally, TIMP3 (a tissue inhibitor of metalloproteinases) was also upregulated. Because inflammation-related genes and inflammatory responses participate in injury processes, our results provide potential targets for molecular therapies aimed at regulating inflammation.

To investigate transcriptome changes after SCI, transcription factors were predicted based on AnimalTFDB 2.0. An interaction network of upregulated DEGs at each time point was then constructed, with transcription factors identified among these DEGs. The network contained 399 regulation nodes, which was far more than expected. The network also showed that 15 transcription factors were involved in 77 nodes. Moreover, 6 of the 10 upregulated DEGs showing the highest degrees of distribution were transcription factors. Consequently, we focused on transcription factors with higher-degree distributions because these DEGs may play important roles in SCI pathogenesis. Among these six transcription factors, ATF3 showed the greatest change (almost 64-fold increase), which began within 30 minutes after SCI and was still high at 28 days. Thus, according to our network analysis, we predict that ATF3 is the key transcription factor, showing long-term increased expression and interacting with many proteins.

As an activating transcription factor, ATF3 may play important roles in inflammation and carcinogenesis in various diseases (Gilchrist et al., 2006; Lin et al., 2014; Huang et al., 2015; Rao et al., 2015; Bar et al., 2016; Iezaki et al., 2016; Mallano et al., 2016; Udayakumar et al., 2016; Wang et al., 2016; Kaitu’u-Lino et al., 2017; Kim et al., 2017). ATF3 can be induced by the toll-like receptor signaling pathway, and in turn, appears to negatively regulate toll-like receptor-stimulated inflammatory responses (Gilchrist et al., 2006). ATF3 can translocate into the nucleus and recruit histone deacetylase-1 to the promoter of interleukin-6 and tumor necrosis factor, ultimately repressing their expression (Gilchrist et al., 2006). Lai et al. (2013) have already shown that ATF3 regulates the release of some inflammatory molecules to reduce lung injury and decrease mortality rate of mice challenged by lipopolysaccharide. ATF3 can also bind DNA sites with Jun to form homo- or hetero-dimers (Tsujino et al., 2000; Koh et al., 2010). The Jun family protein members may participate in the MAPK signaling pathway, which can be induced in response to neuronal injury (Raivich et al., 2004). These genes were also upregulated at every time point and showed higher degrees in our network analysis, which is consistent with previous studies. Additionally, ATF3 can regulate two aspects of the neutrophil response: inhibition of neutrophil chemokine production and promotion of neutrophil chemotaxis (Boespflug et al., 2014). Furthermore, ATF3 can regulate the canonical transforming growth factor-β signaling pathway in systemic sclerosis and might be a potential target for therapy (Mallano et al., 2016). However, its function in the pathogenesis of SCI has not been fully investigated (Seijffers et al., 2007). Tsujino et al. (2000) suggested that ATF3 connects with SOCS3, which is induced in neurons after SCI. ATF3 has also been proposed to be a potent marker of nerve injury and a novel marker for regeneration (Lindå et al., 2011). STAT3 is an important transcription factor during SCI and also a potential target for therapy (Herrmann et al., 2008).

More or less inevitably, patients with traumatic SCI suffer from locomotor deficits throughout their lives. Currently, no pharmacological or biological therapies have proven to be clinically effective. Moreover, collecting spinal cord samples from SCI patients or administering experimental treatments is not possible. To date, scientific investigation of SCI pathology and therapies have depended mostly on animal models of SCI (Kwon et al., 2010). Therefore, animal models provide promising hints regarding effective therapies and the mechanisms underlying SCI. In this study, our results were obtained at multiple time points, with biological repetitions at each time point. Potential genes related to SCI were screened using bioinformatic methods. Therefore, changes in gene expression after SCI were analyzed from multiple levels. Nonetheless, limitations of microarray analyses should be considered. First, we focused on tissue at the injury site, not one specific cell type. Therefore, different cells types will have contributed to RNA expression profiling as a whole, yet may play distinct functions. Targeting specific cells of interest for thorough analyses may yield more valuable results (Greenhalgh and David, 2014; Sofroniew, 2015). Nonetheless, clarification and comprehensive understanding of these molecular events are indispensable at a systemic level. Additionally, the dataset was retrieved from the existent Gene Expression Omnibus Database, and only moderate damage was introduced in our experimental approach. If more degrees of damage and injury sites were analyzed, again the results and conclusions may be more meaningful. Thus, further relevant studies are needed for validation and enhancement of practical significance at various levels, including animal models and cellular experiments. Despite these limitations, our results may be helpful for clarifying the pathological reaction and exploring new therapeutic approaches for SCI.

Taken together, gene expression profiling was noticeably altered at different time points/stages following SCI. The most remarkably upregulated DEGs were found to be associated with inflammation, including transcription factors such as ATF3. These essential genes can be considered as candidate targets for treatment of SCI. Further experiments, including functional studies, are necessary to integrate various types of data and reveal the exact underlying mechanisms in animal models before any clinical use.

Additional files:

Quality control of the CEL files in GSE5296 dataset based on R/Bioconductor program.

Open peer review reports

Click here for additional data file.

Click here for additional data file.

Click here for additional data file.

Additional Table 1

Differentially expressed genes with consistent directions of variation at each time point after T8 SCI

	Symbols
Up-regulated genes	1810011010Rik
	3930401B19Rik
	A130040M12Rik
	Ada
	Adamts1
	Akr1b8
	Arid5a
	Atf3
	Batf3
	Birc3
	C230007H23Rik
	C5ar1
	Casp4
	Ccl12
	Ccl2
	Ccl3
	Ccl4
	Ccl9
	Cd14
	Cebpb
	Cebpd
	Ch25h
	Cp
	Ctla2a
	Ctla2b
	Ctse
	Cxcl1
	Cxcl10
	Cxcl16
	Cyr61
	Cytip
	D17H6S56E-5
	Dusp3
	Dusp6
	Egr2
	Elavil
	Fnl
	Fos
	Fosl2
	Gadd45g
	Galnt2
	Galnt7
	Gch1
	Glipr1
	Gm20186
	Gpr65
	Gpr84
	H2-Aa
	H2-Ab1
	Hmgb2
	Icam1
	Id1
	Ier2
	Irf2bp2
	Itga5
	Itgam
	Jun
	Junb
	Klf6
	Lasp1
	Lcn2
	Lmna
	Lrgl
	Maff
	Map3k8
	Mcll
	Mpzl2
	Ms4a6d
	Myc
	Nes
	Nfkbiz
	Nras
	Nrros
	Ogfr
	Peak1
	Plala
	Plaur
	Plek
	Plin4
	Ptbp3
	Ptgs2
	Rab20
	Rap2b
	Rgs1 Rhoj
	Rrad
	S100a10
	Sbno2
	Sdc4
	Sec61a1
	Serpine1
	Sgk1
	Snx18
	Socs3
	Srgn
	Stat3
	Stom
	Tgifl
	Thbsl
	Timp3
	Tlrl3
	Tlr2
	Tnfrsf12a
	Tribl
	Trp53
	Tubb6
	Wwtrl
	Xdh
	Zfp36
Down-regulated genes	1700056N10Rik
	4833408G04Rik
	6720422M22Rik
	Agbl5
	AI480526
	B830017H08Rik
	BC030500
	CdhlO
	Cilp
	D630023F18Rik
	Dgkg
	Dnajc21
	Dtxl
	Erdrl
	Faml63a
	Gprl2
	Grin2b
	Hemgn
	Midi1
	Pdp2
	Pou4f1
	Prr36
	Rhox4b
	Smtnl2
	Tfdp2
	Trhde`
	Ttc3
	Upp2
	Wwox
	Xist

Additional Table 2

The 3 KEGG pathway terms with the highest changes in enrichment scores at each SCI sub-group

Sub-groups	DEGs	KEGG terms	P	Fold Enrichment	Genes
0.5 h	Upregulated Genes	TNF signaling pathway	0.0000	8.2443	CXCL1, ICAM1, IL6, CEBPB, CCL2, TNF, PTGS2, SOCS3, MAP2K3, CXCL2, EDN1, NFKBIA, …
		Legionellosis	0.0000	7.8827	CXCL1, IL6, TNF, CXCL2, TLR2, IL1B, NFKBIA, HSPA1A, HSPA1B, ITGAM, CD14
		Leishmaniasis	0.0000	7.6587	FOS, TNF, PTGS2, NCF1, JUN, TLR2, IL1B, NFKBIA, H2-AA, H2-AB1, ITGAM, IL1A
	Downregulated Gene	ECM-receptor interaction	0.0223	6.4983	TNC, COL2A1, COL1A1, SV2C
		Amoebiasis	0.0083	6.1095	RAB7B, ACTN1, GNAS, COL2A1, COL1A1
		Proteoglycans in cancer	0.0123	4.2255	LUM, IGF1, COL1A1, FZD2, MMP2, KDR
4 h	Upregulated Genes	Glycosaminoglycan degradation	0.0003	4.8513	GNS, ARSB, NAGLU, HPSE, GUSB, HEXA, HEXB, GALNS, GLB1
		Leishmaniasis	0.0000	4.7755	PTGS2, C3, TLR2, NFKBIA, TLR4, ITGB2, ITGB1, TGFB1, ITGAM, FOS, MYD88, IL1B, IFNGR1, …
		Staphylococcus aureus infection	0.0000	4.7543	ICAM1, C3AR1, C5AR1, C3, FCGR4, ITGB2, H2-AB1, FCGR1, C1QC, ITGAM, FCGR3, C1QA, …
	Downregulated Genes	Steroid biosynthesis	0.0000	18.8585	TM7SF2, CYP51, SC5D, MSMO1, SQLE, DHCR7, LSS, HSD17B7, DHCR24, NSDHL, FDFT1
		Terpenoid backbone biosynthesis	0.0000	9.9138	FNTB, MVD, HMGCR, FDPS, HMGCS1, MVK, IDI1
		Nicotine addiction	0.0000	8.9578	SLC32A1, GABRA1, GRIN2B, GABRA3, GABRB2, GRIN2C, GRIN2D, GRIN1, CHRNA4, CHRNA7, …
24 h	Upregulated Genes	Malaria	0.0000	6.9027	CSF3, SELP, ICAM1, IL6, CCL2, MET, TLR2, ACKR1, TLR4, ITGB2, TGFB1, VCAM1, CCL12, …
		TNF signaling pathway	0.0000	5.9275	CXCL1, CSF2, CCL2, PTGS2, CSF1, CXCL3, CXCL2, EDN1, NFKBIA, NFKB1, CCL5, CXCL10, …
		Legionellosis	0.0000	5.5222	CXCL1, IL6, RELA, CXCL3, HBS1L, CXCL2, TLR2, NFKBIA, NFKB1, TLR4, ITGB2, HSPA1A, …
	Downregulated Gene	Hedgehog signaling pathway	0.0245	6.2764	PTCH1, GAS1, HHIP, GLI2
		Basal cell carcinoma	0.0001	6.1623	BMP4, WNT5A, WNT4, PTCH1, FZD2, HHIP, GLI2, AXIN2, FZD6
		Fatty acid elongation	0.0302	5.7936	HACD4, ACOT2, ELOVL7, ELOVL6
3 d	Upregulated Genes	DNA replication	0.0000	4.7563	POLE, POLA1, MCM2, MCM3, MCM4, MCM5, MCM6, RFC5, POLD3, DNA2, POLD4, RFC3, …
		Staphylococcus aureus infection	0.0000	3.9536	SELP, ICAM1, C3AR1, C5AR1, FCGR4, FPR1, ITGB2, H2-AB1, FPR2, FCGR1, C1QC, ITGAM, …
		Leishmaniasis	0.0000	3.9016	PTPN6, PTGS2, NCF2, NCF1, NCF4, TLR2, FCGR4, NFKBIA, TLR4, ITGB2, H2-AB1, STAT1, …
	Downregulated Gene	Nicotine addiction	0.0000	7.8355	GABRG1, GABRG2, GABRA2, GABRA1, GABRA4, GABRB2, GABRB1, GRIN1, GABRA5, …
		Cocaine addiction	0.0000	5.0497	CDK5R1, DRD2, MAOA, ADCY5, MAOB, GRIN1, GRIN3B, GRIN3A, GRM3, GRIA2, GRIN2B, …
		Retrograde endocannabinoid signa	0.0000	4.8046	ADCY3, GABRB2, GABRB1, ABHD6, ADCY5, GNG13, KCNJ3, RIMS1, GNG8, CNR1, MGLL, …
7 d	Upregulated Genes	DNA replication	0.0000	4.1957	LIG1, POLE, POLA1, MCM2, POLA2, MCM3, MCM4, MCM5, MCM6, PRIM1, RFC5, POLD3, …
		Other glycan degradation	0.0020	3.4964	AGA, MAN2B2, HEXA, HEXB, MAN2B1, FUCA1, MANBA, GBA, GLB1
		Glycosaminoglycan degradation	0.0014	3.3299	GNS, ARSB, NAGLU, HYAL2, HPSE, GUSB, HEXA, HEXB, GALNS, GLB1
	Downregulated Genes	Nicotine addiction	0.0000	7.7072	GABRG2, GABRG3, GABRA1, GABRB3, GABRB2, GABRA3, GABRB1, GRIN1, GABRA5, …
		Steroid biosynthesis	0.0000	6.7607	CYP51, SC5D, MSMO1, SQLE, DHCR7, LSS, HSD17B7, NSDHL, DHCR24, FDFT1
		Retrograde endocannabinoid signa	0.0000	4.8637	ADCY1, ADCY2, GABRB3, GABRB2, GNAI1, GABRB1, ADCY5, GNG13, RIMS1, KCNJ3, …
28 d	Upregulated Genes	Staphylococcus aureus infection	0.0000	4.3360	ITGAL, ICAM1, C3AR1, C5AR1, C3, FCGR4, ITGB2, H2-AB1, FCGR1, C1QC, ITGAM, FCGR3, …
		Leishmaniasis	0.0000	4.2712	PTGS2, C3, TLR2, TLR4, ITGB2, ITGB1, TGFB1, ITGAM, IRAK4, FOS, MYD88, IFNGR2, …
		Other glycan degradation	0.0016	4.1894	MAN2B2, HEXA, HEXB, MAN2B1, FUCA1, MANBA, GBA, GLB1
	Downregulated Genes	Steroid biosynthesis	0.0000	15.2542	TM7SF2, CYP51, SC5D, MSMO1, SQLE, DHCR7, LSS, HSD17B7, NSDHL, DHCR24, FDFT1
		Nicotine addiction	0.0000	9.8805	GABRG3, GABRA1, GABRA3, GABRB2, GABRA5, GRIN1, GRIA4, GRIN3A, SLC32A1, …
		Terpenoid backbone biosynthesis	0.0014	6.8734	MVD, HMGCR, FDPS, HMGCS1, MVK, IDI1

Additional Table 3

The 3 GO terms with the highest changes in enrichment scores at each SCI sub-group

Sub-groups	DEGs	Category	GO terms	P	Fold Enrichment	Genes
0.5 h	Upregulated Genes	GO_BP	positive regulation of TRAIL-activated apoptotic signaling pathway	0.0012	49.0027	ATF3, TIMP3, PTEN
			cellular response to cycloheximide	0.0012	49.0027	KLF2, MYC, KLF4
			positive regulation of cell aging	0.0403	49.0027	TRP53, LMNA
		GO_CC	neuronal cell body membrane	0.0023	8.8052	ATP2B1, KCNA2, KCNB1, HPCA, CX3CR1
			microtubule associated complex	0.0153	7.5659	NDEL1, MAP1B, PAFAH1B1, RANBP2
			nuclear euchromatin	0.0203	6.8094	JUN, MYC, KLF4, SMARCA4
		GO_MF	C-4 methylsterol oxidase activity	0.0424	46.5227	MSMO1, CH25H
			CCR2 chemokine receptor binding	0.0027	34.8920	CCL12, CCL2, CCL7
			MRF binding	0.0065	23.2613	TSC22D3, CREBBP, BHLHE40
	Downregulated Genes	GO_BP	intramembranous ossification	0.0005	86.9327	COL1A1, AXIN2, MMP2
			endothelium development	0.0282	69.5462	SLC40A1, KDR
			peripheral nervous system axon regeneration	0.0282	69.5462	TNC, MMP2
		GO_CC	cytoplasmic microtubule	0.0360	9.9504	BIRC5, MID1, AXIN2
			nuclear envelope	0.0096	5.9828	TRPC2, AGPAT5, CHIL3, RRM2, TFDP2
			extracellular matrix	0.0009	5.1444	LPL, LUM, TNC, CILP, COL2A1, COL1A1, MFAP4, …
		GO_MF	tRNA (guanine-N2-)-methyltransferase activity	0.0171	115.1551	TRMT1L, TRMT11
			glycine binding	0.0033	34.5465	GLRA1, GRIN2B, SRR
			integrin binding	0.0198	6.9093	FAP, IGF1, ACTN1, KDR
4 h	Upregulated Genes	GO_BP	antigen processing and presentation of exogenous peptide antigen via MHC class I	0.0001	13.8031	H2-K1, IFI30, FCER1G, FCGR1, FCGR3
			antigen processing and presentation of peptide antigen	0.0014	13.8031	SLC11A1, H2-AA, H2-AB1, CTSS
			toll-like receptor 7 signaling pathway	0.0014	13.8031	HAVCR2, UNC93B1, PIK3AP1, TLR7
		GO_CC	NLRP3 inflammasome complex	0.0032	11.3000	CASP4, GSDMD, PYCARD, CASP1
			AIM2 inflammasome complex	0.0273	10.5938	CASP4, PYCARD, CASP1
			CD95 death-inducing signaling complex	0.0273	10.5938	CFLAR, CASP8, FAS
		GO_MF	IgG receptor activity	0.0153	13.6510	FCGRT, FCGR1, FCGR3
			peptidase activator activity involved in apoptotic process	0.0153	13.6510	PYCARD, CTSC, CTSH
			Fc-gamma receptor I complex binding	0.0153	13.6510	LGALS3, FGR, FLNA
	Downregulated Genes	GO_BP	regulation of synaptic transmission, dopaminergic	0.0059	23.2616	CNTNAP4, PNKD, CHRNA7
			positive regulation of long term synaptic depression	0.0097	18.6093	KCNB1, STAU2, IQSEC2
			regulation of cAMP-mediated signaling	0.0097	18.6093	PDE4A, PEX5L, RASD2
		GO_CC	cone cell pedicle	0.0058	23.5944	SLC32A1, RAPGEF4, PCLO
			synaptic cleft	0.0016	15.7296	CDH8, DNM3, GRIN2B, GRIN1
			spectrin-associated cytoskeleton	0.0139	15.7296	ANK1, ANK3, SPTB
		GO_MF	voltage-gated cation channel activity	0.0099	18.3642	GRIK1, GRIN2D, GRIN1
			NMDA glutamate receptor activity	0.0017	15.3035	GRIN2B, GRIN2C, GRIN2D, GRIN1
			oxidoreductase activity, acting on the CH-CH group of donors, NAD or NADP as acceptor	0.0146	15.3035	TM7SF2, DHCR7, DHCR24
24 h	Upregulated Genes	GO_BP	positive regulation of platelet activation	0.0000	20.5477	PTPRJ, SELP, PLEK, PDPN, TLR4
			wound healing involved in inflammatory response	0.0069	20.5477	CD44, HMOX1, CCR2
			positive regulation of monocyte aggregation	0.0069	20.5477	CD44, HAS2, NR4A3
		GO_CC	CHOP-C/EBP complex	0.0062	21.5356	CEBPA, CEBPB, DDIT3
			I-kappaB/NF-kappaB complex	0.0009	17.2285	RELA, RELB, NFKB1, NFKB2
			lamin filament	0.0196	12.9214	EIF6, LMNB1, LMNA
		GO_MF	CCR2 chemokine receptor binding	0.0141	14.9537	CCL12, CCL2, CCL7
			C-C chemokine binding	0.0022	13.2922	ZFP36, CCR5, CCR1, ACKR1
			interleukin-1 receptor activity	0.0022	13.2922	IL1R2, IL1R1, IL18RAP, IL1RAP
	Downregulated Genes	GO_BP	intramembranous ossification	0.0005	22.8308	CTSK, COL1A1, AXIN2, MMP2
			regulation of endothelial cell proliferation	0.0005	22.8308	BMP4, ALDH1A2, TEK, KDR
			urinary bladder development	0.0080	20.5477	WNT5A, RBP4, TRP63
		GO_CC	synaptic cleft	0.0195	13.4059	CDH8, GRIN2B, C1QL1
			AP-2 adaptor complex	0.0246	11.9164	EPS15, TBC1D5, SGIP1
			microfibril	0.0362	9.7498	LTBP1, FBN2, MFAP4
		GO_MF	GPI-linked ephrin receptor activity	0.0156	14.9537	EPHA5, EPHA7, EPHA3
			platelet-derived growth factor binding	0.0003	14.5383	COL3A1, COL1A2, PDGFRB, COL2A1, COL1A1
			G-protein coupled purinergic nucleotide receptor activity	0.0054	10.7360	P2RY12, P2RY13, GPR34, P2RY14
3 d	Upregulated Genes	GO_BP	antigen processing and presentation of exogenous peptide antigen via MHC class I	0.0002	12.7159	H2-K1, IFI30, FCER1G, FCGR1, FCGR3
			antigen processing and presentation of peptide antigen	0.0018	12.7159	SLC11A1, H2-AA, H2-AB1, CTSS
			positive regulation of neutrophil apoptotic process	0.0175	12.7159	CD44, ANXA1, HCAR2
		GO_CC	CHOP-C/EBP complex	0.0162	13.2493	CEBPA, CEBPB, DDIT3
			interleukin-6 receptor complex	0.0162	13.2493	IL6, IL6ST, IL6RA
			collagen type V trimer	0.0162	13.2493	COL5A3, COL5A2, COL5A1
		GO_MF	Fc-gamma receptor I complex binding	0.0178	12.6146	LGALS3, FGR, FLNA
			peptidase activator activity involved in apoptotic process	0.0178	12.6146	PYCARD, CTSC, CTSH
			CCR5 chemokine receptor binding	0.0001	10.8125	NES, CNIH4, STAT1, CCL5, CCL4, STAT3
	Downregulated Genes	GO_BP	spontaneous neurotransmitter secretion	0.0111	16.0729	RIMS2, RIMS1, STX1B
			development of primary male sexual characteristics	0.0111	16.0729	WNT5A, SFRP1, SFRP2
			positive regulation of protein kinase C activity	0.0213	12.0547	WNT5A, AGT, CEMIP
		GO_CC	GABA receptor complex	0.0112	16.0114	GABRG2, GABRA1, GABRB2
			synaptic cleft	0.0000	12.0086	CDH8, GRIN2B, GRIN1, C1QL1, GRIA3, ADGRB3
			kainate selective glutamate receptor complex	0.0343	9.6068	GRIK1, GRIK5, GRIA4
		GO_MF	microsatellite binding	0.0112	16.0202	HEY1, HEY2, HEYL
			neurotransmitter binding	0.0000	12.4602	GRIN2B, SLC6A11, GRIN2D, SLC6A13, GRIN1, …
			NMDA glutamate receptor activity	0.0000	12.0152	GRIN2B, GRIN2C, GRIN2D, GRIN1, GRIN3B, …
7 d	Upregulated Genes	GO_BP	glycosaminoglycan metabolic process	0.0001	8.2679	GNS, NDST1, HEXA, HEXB, FOXC1, CLN6
			positive regulation of platelet activation	0.0010	8.2679	PTPRJ, SELP, PLEK, PDPN, TLR4
			toll-like receptor 7 signaling pathway	0.0064	8.2679	HAVCR2, UNC93B1, PIK3AP1, TLR7
		GO_CC	death-inducing signaling complex	0.0001	8.4970	CFLAR, CASP3, RIPK1, CASP8, FADD, FAS
			ripoptosome	0.0001	8.4970	CFLAR, RIPK1, TICAM1, CASP8, RIPK3, FADD
			condensin complex	0.0009	8.4970	NCAPH, NCAPG, SMC2, SMC4, NCAPD2
		GO_MF	cysteine-type endopeptidase activity involved in execution phase of apoptosis	0.0398	8.3235	CFLAR, CASP3, CASP7
			Fc-gamma receptor I complex binding	0.0398	8.3235	LGALS3, FGR, FLNA
			nicotinate-nucleotide diphosphorylase (carboxylating) activity	0.0398	8.3235	NAMPT, QPRT, NAPRT
	Downregulated Genes	GO_BP	positive regulation of long term synaptic depression	0.0002	11.6884	PPP1R9A, MAPT, KCNB1, STAU2, IQSEC2
			regulation of cAMP-mediated signaling	0.0002	11.6884	GNAI1, PDE4A, PDE10A, PEX5L, RASD2
			regulation of synaptic transmission, dopaminergic	0.0023	11.6884	CNTNAP4, PNKD, CHRNB2, CHRNA7
		GO_CC	GABA receptor complex	0.0212	11.5455	GABRG2, GABRA1, GABRB2
			juxtaparanode region of axon	0.0000	9.2364	EPB41L3, KCNAB2, KCNAB1, KCNA2, KCNA1, …
			synaptic cleft	0.0002	8.6591	CDH8, DNM3, GRIN2B, GRIN1, GRIA3, ADGRB3
		GO_MF	netrin receptor activity	0.0023	11.6774	DCC, UNC5B, UNC5A, UNC5C
			AMPA glutamate receptor activity	0.0023	11.6774	GRIA2, GRIA1, GRIA3, GRIA4
			voltage-gated calcium channel activity involved in AV node cell action potential	0.0207	11.6774	CACNA1G, CACNB2, CACNA1C
28 d	Upregulated Genes	GO_BP	response to interferon-beta	0.0000	11.6059	PLSCR1, IKBKE, IFNAR2, BST2, IFITM1, IFITM2, …
			positive regulation of apoptotic cell clearance	0.0024	11.6059	CCL2, C3, MFGE8, C2
			positive regulation of neutrophil apoptotic process	0.0210	11.6059	CD44, ANXA1, HCAR2
		GO_CC	CHOP-C/EBP complex	0.0200	11.9019	CEBPA, CEBPB, DDIT3
			ripoptosome	0.0006	9.9183	CFLAR, RIPK1, TICAM1, CASP8, RIPK3
			NLRP3 inflammasome complex	0.0052	9.5215	CASP4, GSDMD, PYCARD, CASP1
		GO_MF	IgG receptor activity	0.0212	11.5231	FCGRT, FCGR1, FCGR3
			epidermal growth factor binding	0.0212	11.5231	EGFR, SEC61B, SHC1
			peptidase activator activity involved in apoptotic process	0.0212	11.5231	PYCARD, CTSC, CTSH
	Downregulated Genes	GO_BP	positive regulation of long term synaptic depression	0.0006	19.9801	PPP1R9A, KCNB1, STAU2, IQSEC2
			JUN phosphorylation	0.0091	18.7314	MAPK8, MAPK10, MAPK8IP1
			regulation of synaptic transmission, dopaminergic	0.0091	18.7314	CNTNAP4, PNKD, CHRNA7
		GO_CC	cone cell pedicle	0.0089	18.9058	SLC32A1, RAPGEF4, PCLO
			potassium channel complex	0.0000	13.7497	KCNH1, KCNAB2, KCNAB1, KCNA2, KCNA6, …
			NMDA selective glutamate receptor complex	0.0001	12.6038	DLGAP3, GRIN2B, GRIN2C, GRIN2D, GRIN1, …
		GO_MF	netrin receptor activity	0.0091	18.7457	DCC, UNC5B, UNC5C
			NMDA glutamate receptor activity	0.0002	15.6214	GRIN2B, GRIN2C, GRIN2D, GRIN1, GRIN3A
			voltage-gated cation channel activity	0.0147	14.9966	GRIK1, GRIN2D, GRIN1

Table 1

Top-20 KEGG terms of upregulated differentially expressed genes with the same variation tendency at each time point after spinal cord injury

KEGG terms	P	Fold enrichment	Gene symbols
TNF signaling pathway	0	14.85	CXCL1, FOS, CCL12, ICAM1, CEBPB, CCL2, PTGS2, SOCS3, JUN, MAP3K8, BIRC3, JUNB, CXCL10
Rheumatoid arthritis	0	13.67	CCL12, FOS, ICAM1, CCL3, CCL2, JUN, TLR2, H2-AA, H2-AB1
Leishmaniasis	0	13.62	FOS, PTGS2, JUN, TLR2, H2-AA, H2-AB1, ITGAM
Malaria	0	12.97	CCL12, ICAM1, CCL2, TLR2, THBS1
Thyroid cancer	0.02	12.88	TRP53, NRAS, MYC
Staphylococcus aureus infection	0	12.45	ICAM1, C5AR1, H2-AA, H2-AB1, ITGAM
Bladder cancer	0	12.15	TRP53, NRAS, THBS1, MYC
Inflammatory bowel disease	0	10.55	JUN, TLR2, H2-AA, H2-AB1, STAT3
Toll-like receptor signaling pathway	0	9.86	FOS, CCL3, JUN, MAP3K8, TLR2, CCL4, CD14, CXCL10
Salmonella infection	0	9.58	CXCL1, FOS, CCL3, JUN, CCL4, CD14
Legionellosis	0.01	8.74	CXCL1, TLR2, ITGAM, CD14
Chagas disease (American trypanosomiasis)	0	8.46	CCL12, FOS, CCL3, CCL2, JUN, SERPINE1, TLR2
Pertussis	0	8.41	FOS, ITGA5, JUN, ITGAM, CD14
Colorectal cancer	0.01	7.78	TRP53, FOS, JUN, MYC
p53 signaling pathway	0.02	7.43	TRP53, SERPINE1, GADD45G, THBS1
Small cell lung cancer	0	7.41	TRP53, PTGS2, BIRC3, MYC, FN1
Hepatitis B	0	6.82	TRP53, NRAS, FOS, EGR2, JUN, TLR2, MYC, STAT3
Prolactin signaling pathway	0.02	6.82	NRAS, FOS, SOCS3, STAT3
mmProteoglycans in cancer	0	6.75	TRP53, NRAS, ITGA5, TLR2, THBS1, SDC4, MYC, TIMP3, STAT3, FN1, PLAUR
Phagosome	0	6.44	ITGA5, TLR2, TUBB6, H2-AA, H2-AB1, THBS1, ITGAM, SEC61A1, CD14

P < 0.05, sorted by fold enrichment scores. KEGG: Kyoto Encyclopedia of Genes and Genomes; TNF: tumor necrosis factor; STAT: signal transducer and activator of transcription.

Table 2

Top-three GO terms with highest fold enrichment scores of differentially expressed genes with the same variation tendency at each time point after spinal cord injury

		GO terms	P	Fold enrichment	Gene symbols
Upregulated genes	BP	Response to vitamin b3	0.01	179.03	TRP53, CCL2
		Positive regulation of cell aging	0.01	179.03	TRP53, LMNA
		Negative regulation of natural killer cell chemotaxis	0.02	119.35	CCL12, CCL2
	CC	Ribonucleoprotein complex	0.02	127.26	ZFP36, ELAVL1
		Fibrinogen complex	0.04	54.54	THBS1, FN1
		Phagocytic vesicle membrane	0.02	12.18	TLR2, RAB20, SEC61A1
	MF	Ccr2 chemokine receptor binding	0.02	87.23	CCL12, CCL2
		Ccr5 chemokine receptor binding	0	74.77	NES, CCL4, STAT3
		Lipoteichoic acid binding	0.03	58.15	TLR2, CD14
Downregulated genes	BP	Suckling behavior	0.02	111.96	GRIN2B, POU4F1
	MF	Ligase activity	0.04	8.5	DTX1, MID1, TTC3
		Zinc ion binding	0	6.68	TRHDE, GRIN2B, DTX1, AGBL5, DNAJC21, MID1, TTC3
		Metal ion binding	0.01	2.75	TRHDE, PDP2, GRIN2B, DTX1, DGKG, AGBL5, MID1, TTC3, CDH10

P < 0.05, sorted by fold enrichment score. GO: Gene Ontology; STAT: signal transducer and activator of transcription.