The transcriptome data under accession number GSE5296 based on the GPL1261 platform (Affymetrix Mouse Genome 430 2.0 Array; Affymetrix; Thermo Fisher Scientific, Inc., Waltham, MA, USA) was obtained from the National Center for Biotechnology Information Gene Expression Omnibus (GEO) database (www.ncbi.nlm.nih.gov/geo). In the original dataset (GSE5296), C57BL6 mice were subjected to moderate contusion injury at the T8 spinal segment. Sections of the spinal cord (4 mm in length) were analyzed from the site of the trauma and at the immediately adjacent R and C regions, at 0.5, 4, 24 and 72 h, and 7 and 28 days following injury. A total of 96 chips were available in this dataset data, including 18 SCI samples for each M, R and C region (n=3/time-point), 12 sham-injury samples for each M, R and C region (n=2/time-point), and two samples in each region obtained from naive mice (data not shown).

The Robust Multichip Average algorithm in the Oligo package, version 1.42.0 (http://www.bioconductor.org) was used to preprocess the raw transcriptome data included in the GSE5296 dataset (13). Data were subjected to background correction, normalization, probe summary and log₂ transformation. If there were several probes annotated to same gene symbol, the average value was used to represent the expression level of this gene. There were 45,037 probes in the raw data and 21,812 genes remained following data processing.

In the present study, data were divided into the following paired groups: i) Post-SCI group vs. sham group in the M region at different time-points (0.5, 4, 24 and 72 h, and 7 and 28 days); ii) post-SCI group vs. post-sham group in the R region at different time-points (24 and 72 h, and 7 days); and iii) post-SCI group vs. sham group in the C region at different time-points (24 and 72 h, and 7 days). Fold change |log₂FC| and P-values from a Student's t-test were used to identify the DEGs. An average fold-change >2.0 and P<0.05 were used as cutoff criteria. Subsequently, Venny 2.1 (bioinfogp.cnb.csic.es/tools/venny/index.html) was used to compare lists of DEGs and to construct Venn diagrams.

To identify pathways and biological processes enriched by DEGs, the Database for Annotation, Visualization and Integrated Discovery (DAVID 6.8; http://david.abcc.ncifcrf.gov) was used to perform Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway and Gene Ontology (GO) functional analysis (14). GO terms were identified under categories of biological process. The threshold was set to P<0.05.

A PPI network was constructed using the Search Tool for the Retrieval of Interacting Genes/Proteins (String 10.5; www.string-db.org). The network was visualized using the Cytoscape software platform (Cytoscape 3.6) based on functional analysis information, including fold-change in gene/protein expression, PPIs and GO/KEGG pathway enrichment (15). A default confidence cutoff of 400 was used in the present study. Experimentally determined interactions were presented as solid lines between genes/proteins, and the dashed lines represent database predicted interactions (16).

In the present study, neuronal function and synaptic transmission-associated genes were defined as cholinergic receptor nicotinic α7 subunit (Chrna7), synapsin II (Syn2), potassium voltage-gated channel subfamily C member 1 (Kcnc1), ATPase plasma membrane Ca²⁺ transporting 2 (Atp2b2), unc-13 homolog A (Unc13a), regulating synaptic membrane exocytosis 1 (Rims1), calcium/calmodulin dependent protein kinase IIγ (Camk2g), calcium/calmodulin dependent protein kinase IIα (Camk2a), thyrotropin releasing hormone receptor (Trhr) and glutamate metabotropic receptor 1 (Grm1) (3,4,9). Furthermore, based on literature review, inflammation-associated genes were identified, including interleukin (IL)1β, IL6, IL7, IL4, IL10, CD44 molecule (Indian blood group) (Cd44), cytochrome b-245 β-chain (Cybb), intercellular adhesion molecule 1 (Icam1), cytochrome b-245 α-chain (Cyba), HCK proto-oncogene, Src family tyrosine kinase (Hck), caspase 1 (Casp1), transforming growth factor β1 (Tgfb1), Rac family small GTPase 2 (Rac2), integrin subunit β2 (Itgb2) and C-X-C motif chemokine receptor 4 (Cxcr4) (3,4,9). The fold-change of expression of each gene (post-SCI data vs. sham data in each region) at three different time-points (24 and 72 h, and 7 days) was determined.

Box plots presenting post-SCI and sham surgery data in the different regions at different time points (M region: 0.5, 4, 24 and 72 h, and 7 and 28 days; R and C regions: 24 and 72 h, and 7 days) prior to and following data normalization are presented in Fig. 1. The results demonstrated that the gene expression values in each sample were similar following normalization. Following data pre-processing, DEGs between the post-SCI groups and sham groups in the M region at six time points were analyzed. There were 226, 495, 566, 816, 1,286 and 1,023 upregulated DEGs at 0.5, 4, 24, and 72 h and 7 and 28 days, respectively. Additionally, a total of 148, 212, 756, 659, 438 and 319 downregulated DEGs were identified at the six time points, respectively (Fig. 2). There was an increased number of upregulated DEGs compared with downregulated DEGs at the five time points (0.5, 4 and 72 h, and 7 and 28 days) and the number of upregulated DEGs reached a peak on day 7. By contrast, the number of downregulated DEGs peaked at 24 h and subsequently decreased over time. The above results indicated that gene expression alterations occurred primarily at 24 and 72 h, and 7 days following injury.

A number of studies have investigated alterations in gene expression between 0.5 and 6 h following SCI in mice (17,18). The present study focused on DEGs at 24 and 72 h, and 7 days following injury. The top five enriched KEGG pathways of up- and downregulated DEGs at each time point (24 and 72 h, and 7 days) are presented in Fig. 3A and B, respectively. At 24 h, upregulated DEGs were enriched in pathways including ‘extracellular matrix-receptor interaction’ (P=2.08×10⁻⁶), ‘focal adhesion’ (P=3.33×10⁻⁶) and ‘regulation of the actin cytoskeleton’ (P=7.98×10⁻⁶). At 72 h and 7 days, upregulated DEGs were primarily involved in inflammation- and immunity-associated pathways. ‘Leukocyte transendothelial migration’ (P=1.05×10⁻⁹) and ‘natural killer cell-mediated cytotoxicity’ (P=1.02×10⁻⁵) were enriched at 72 h. ‘Lysosome’ (P=5.57×10⁻¹⁰), ‘Toll-like receptor signaling pathway’ (P=5.93×10⁻⁸) and ‘leukocyte transendothelial migration’ (P=2.48×10⁻⁷) were enriched at 7 days following injury. Notably, downregulated DEGs were enriched in the neuronal function and synaptic transmission-associated pathways, including ‘neuroactive ligand-receptor interactions’ (P=1.43×10⁻⁵) at 24 h, ‘synaptic vesicle cycle’ (P=2.32×10⁻⁷) and ‘calcium signaling pathway’ (P=3.45×10⁻⁷) at 72 h, and ‘neuroactive ligand-receptor interaction’ (P=4.95×10⁻⁵) and ‘glutamatergic synapse’ (P=2.23×10⁻⁴) at 7 days following injury. Other downregulated DEGs were most significantly enriched in ‘steroid biosynthesis’ (P=1.08×10⁻⁵).

The GO terms (biological process) of up- and downregulated DEGs are summarized in Table I. The results demonstrated that upregulated DEGs were most enriched in ‘response to stress’ (P=2.30×10⁻²⁷) at 24 h. Additionally, ‘immune system process’ was the most enriched function at 72 h and 7 days (P=3.59×10⁻²⁷ and P=7.77×10⁻⁶⁸, respectively). By contrast, downregulated DEGs were most enriched in ‘single-organism cellular process’ (P=9.30×10⁻⁸) at 24 h following injury. Furthermore, downregulated DEGs were enriched in neuronal function- and synaptic transmission-associated biological process terms at 72 h and 7 days following injury. At 72 h ‘synaptic signaling’ (P=1.21×10⁻³⁰) and ‘chemical synaptic transmission’ (P=7.68×10⁻²⁵) were most enriched. At 7 days following injury, ‘synaptic signaling’ (P=5.83×10⁻⁹) and ‘anterograde trans-synaptic signaling’ (P=2.51×10⁻⁸) were most enriched. ‘Sterol biosynthetic process’ (P=4.96×10⁻¹⁰) was the most significantly enriched biological process at 7 days following injury. Collectively, in the present study, upregulated DEGs were predominantly associated with immune and inflammatory functions, while downregulated DEGs were involved in neuronal function, synaptic transmission and steroid biosynthesis.

Venn diagram analysis of DEGs at 24 and 72 h, and 7 days is presented in Fig. 4A. There were 138 common upregulated and 20 overlapping downregulated DEGs at these three time points. Subsequently, a more comprehensive bioinformatics analysis was performed using Cytoscape software, a tool for predicting PPI networks. The results revealed that immune and inflammatory functions were enriched in co-upregulated DEGs (Fig. 4B). Additionally, neuronal functions and ‘steroid hormone biosynthesis’ were enriched in the co-downregulated DEGs (Fig. 4C).

The results of pathway, functional enrichment and PPI network analyses demonstrated that inflammatory- and neuronal-associated functions serve roles in post-SCI mice in the M region at three different time points (24 and 72 h, and 7 days). Alterations in these functions were further investigated via temporal and spatial analysis of the expression of numerous genes. Neuronal function- and synaptic transmission-associated genes, including Chrna7, Atp2b2, Rims1, Camk2g, and Trhr, were downregulated at 72 h compared with the expression at 24 h post-SCI in the M region (Fig. 5A). Furthermore, neuronal function associated genes, Chrna7, Syn2 and Unc13a, were significantly upregulated on day 7 compared with the expression at 72 h post-SCI. Similar alterations were observed in R and C regions at three different time points post-SCI, although the overall fold-changes with time were minimal (Fig. 5A). Genes involved in inflammatory processes exhibited different and complex alterations (Fig. 5B). The expression levels of inflammatory-associated genes (IL1b and IL6) decreased significantly at 72 h following injury in the M region compared with levels at the 24 h time interval. However, similar alterations were not observed in the R and C regions. By contrast, the expression levels of other genes associated with inflammatory processes (Cd44, Cybb, Cyba, Hck, Casp1, Itgb2 and Cxcr4) significantly increased at 7 days post-SCI in the M region compared with expression levels at the 72 h time interval. Similar trends were observed in the C region, although not in the R region (Fig. 5C). Therefore, the results of the present study suggest that alterations in expression of genes associated with inflammatory and neuronal functions were primarily observed in the M region at different time points post-SCI and these two events occurred in a temporally and spatially specific manner respectively, as reflected in the trend charts.