This article is mentioned in:

Introduction

Spinal cord injury (SCI) entails severe physical and social consequences for patients and their families, and may be of a traumatic or non-traumatic etiology (1,2). Traumatic SCI, which is primarily caused by traffic accidents, occurs when an external impact acutely damages the spinal cord, and may be temporally divided into four stages: i) Acute at <48 h following injury; ii) subacute at 48 h-14 days following injury; iii) intermediate at 14 days-6 months following injury; and iv) chronic at >6 months following injury (1). Although numerous studies have focused on neuroprotective and neuroregenerative therapies, no major breakthroughs have been achieved (3–5). Therefore, it is necessary to elucidate the underlying mechanism of SCI.

Pathophysiologically, the initial mechanical injury damages neurons and initiates a complex secondary injury cascade, leading to progressive cell death, ischemia and inflammation (6). It has been demonstrated that the transcriptome may reflect the pathophysiological state of the cell (7). In a number of recent bioinformatics studies, transcriptome analysis at different time points post-SCI was performed and various molecular events were characterized (8–12). The immune response, inflammatory-associated functions, vasculature development and neurological functions were demonstrated to serve roles in the development of SCI (6–12). However, certain molecular alterations that occur in a temporal and spatial manner remain to be elucidated, particularly those that occur during the secondary injury process.

In the present study, transcriptome data under accession no. GSE5296 was used to identify SCI-specific molecular programs. Temporally, three different time points were evaluated, including 24 and 72 h, and 7 days post-injury for pathway and functional enrichment analysis. Furthermore, a protein-protein interaction (PPI) network was constructed. Spatially, differentially expressed genes (DEGs) in different locations, including the trauma site (M), and immediately adjacent rostral (R) and caudal (C) regions were determined at the aforementioned time points. The present study revealed molecular mechanisms that may be associated with SCI and provided an insight into potential therapeutic targets for treatment of SCI.

Materials and methods

Transcriptome data

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).

Data preprocessing

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 log2 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.

Identification and analysis of DEGs

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 |log2FC| 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.

Pathway and functional enrichment analysis

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.

PPI network construction

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).

Trend charts of neuronal function- and inflammatory response-associated genes

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 Ca2+ 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.

Results

Data preprocessing and DEG screening

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.

Figure 1. Data normalization of differentially expressed genes. Box plots of gene expression in the post-SCI and sham groups from the site of trauma (A) prior to and (B) following normalization. Box plots of gene expression in the post-SCI and sham groups at the immediately adjacent R and C regions (C) prior to and (D) following normalization. R, rostral; C, caudal; M, trauma site; SCI, spinal cord injury.

Figure 2. Differentially expressed genes between post-spinal cord injury and sham samples at each time-point following spinal cord injury.

KEGG pathway and GO enrichment analysis of up- and downregulated DEGs

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−6), ‘focal adhesion’ (P=3.33×10−6) and ‘regulation of the actin cytoskeleton’ (P=7.98×10−6). At 72 h and 7 days, upregulated DEGs were primarily involved in inflammation- and immunity-associated pathways. ‘Leukocyte transendothelial migration’ (P=1.05×10−9) and ‘natural killer cell-mediated cytotoxicity’ (P=1.02×10−5) were enriched at 72 h. ‘Lysosome’ (P=5.57×10−10), ‘Toll-like receptor signaling pathway’ (P=5.93×10−8) and ‘leukocyte transendothelial migration’ (P=2.48×10−7) 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−5) at 24 h, ‘synaptic vesicle cycle’ (P=2.32×10−7) and ‘calcium signaling pathway’ (P=3.45×10−7) at 72 h, and ‘neuroactive ligand-receptor interaction’ (P=4.95×10−5) and ‘glutamatergic synapse’ (P=2.23×10−4) at 7 days following injury. Other downregulated DEGs were most significantly enriched in ‘steroid biosynthesis’ (P=1.08×10−5).

Figure 3. Kyoto Encyclopedia of Genes and Genomes pathways of (A) upregulated and (B) downregulated differentially expressed genes at three time-points (24 and 72 h, and 7 days) following spinal cord injury.

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−27) at 24 h. Additionally, ‘immune system process’ was the most enriched function at 72 h and 7 days (P=3.59×10−27 and P=7.77×10−68, respectively). By contrast, downregulated DEGs were most enriched in ‘single-organism cellular process’ (P=9.30×10−8) 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−30) and ‘chemical synaptic transmission’ (P=7.68×10−25) were most enriched. At 7 days following injury, ‘synaptic signaling’ (P=5.83×10−9) and ‘anterograde trans-synaptic signaling’ (P=2.51×10−8) were most enriched. ‘Sterol biosynthetic process’ (P=4.96×10−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.

Table I. GO terms enriched by differentially expressed genes at three time-zpoints following spinal cord injury.

Table I.

GO terms enriched by differentially expressed genes at three time-zpoints following spinal cord injury.

A, Injury vs. sham (24 h)
Category	Term	Biological process	No. genes	P-value
Upregulated
	GO:0006950	Response to stress	144	2.30×10−27
	GO:0070887	Cellular response to chemical stimulus	107	9.95×10−22
	GO:0050896	Response to stimulus	233	3.47×10−21
	GO:0010033	Response to organic substance	110	9.68×10−21
	GO:0006952	Defense response	73	3.99×10−20
Downregulated
	GO:0044763	Single-organism cellular process	329	9.30×10−8
	GO:0044699	Single-organism process	363	1.70×10−7
	GO:0048512	Circadian behavior	9	3.21×10−7
	GO:0007622	Rhythmic behavior	9	6.19×10−7
	GO:0007275	Multicellular organism development	154	1.11×10−6
B, Injury vs. sham (72 h)
Category	Term	Biological process	No. genes	P-value
Upregulated
	GO:0002376	Immune system process	148	3.59×10−27
	GO:0006952	Defense response	90	4.21×10−18
	GO:0006955	Immune response	87	5.70×10−18
	GO:0006954	Inflammatory response	52	9.35×10−18
	GO:0002684	Positive regulation of immune system process	64	6.12×10−16
Downregulated
	GO:0099536	Synaptic signaling	56	1.21×10−30
	GO:0098916	Anterograde trans-synaptic signaling	54	3.26×10−29
	GO:0099537	Trans-synaptic signaling	54	3.65×10−29
	GO:0007268	Chemical synaptic transmission	57	7.68×10−25
	GO:0007267	Cell-cell signaling	74	2.96×10−20
C, Injury vs. sham (day 7)
Category	Term	Biological process	No. genes	P-value
Upregulated
	GO:0002376	Immune system process	286	7.77×10−68
	GO:0006952	Defense response	188	7.30×10−52
	GO:0006950	Response to stress	319	6.25×10−45
	GO:0006955	Immune response	166	2.74×10−41
	GO:0002682	Regulation of immune system process	164	1.22×10−39
Downregulated
	GO:0016126	Sterol biosynthetic process	9	4.96×10−10
	GO:0099536	Synaptic signaling	19	5.83×10−9
	GO:1902653	Secondary alcohol biosynthetic process	7	2.09×10−8
	GO:0098916	Anterograde trans-synaptic signaling	18	2.51×10−8
	GO:0099537	Trans-synaptic signaling	18	2.59×10−8

[i] GO, Gene Ontology.

PPI network construction and functional module analysis

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).

Figure 4. Networks of co-expressed genes and PPIs. Venn diagrams of DEGs at 24 and 72 h and 7 days post-SCI, including (A) upregulated and downregulated DEGs following SCI. PPI networks based on (B) co-upregulated and (C) co-downregulated DEGs. The gradient colors from green to red indicate increasing P-values. Circular nodes indicate genes/proteins. Rectangular nodes indicate a biological process or a KEGG pathway. DEGs, differentially expressed genes; PPI, protein-protein interaction; SCI, spinal cord injury; KEGG, Kyoto Encyclopedia of Genes and Genomes.

Trends in the expression levels of neuronal function- and inflammatory response-associated genes

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.

Figure 5. Trend charts of gene expression in R, M and C regions. The trend charts of expression of genes associated with (A) neuronal function, including Chrna7, Syn2, Kcnc1, Atp2b2, Unc13a, Rims1, Camk2g, Camk2a, Trhr and Grm1 and inflammatory response, including (B) IL1b, IL6, IL7, IL4 and IL10, and (C) CD44, Icam1, Cybb, Cyba, Hck, Casp1, Tgfb1, Rac2, Itgb2 and Cxcr4. Significant differences between expression levels at different time points were identified (*P<0.05 and **P<0.01, 24 h vs. 72 h; #P<0.05 and ##P<0.01, 72 h and 7 days). R, rostral; C, caudal; M, trauma site. IL, interleukin; Cd44, CD44 molecule (Indian blood group); Cybb, cytochrome b-245 β-chain; Icam1, intercellular adhesion molecule 1; Cyba, cytochrome b-245 α-chain; Hck, HCK proto-oncogene, Src family tyrosine kinase; Casp1, caspase 1; Tgfb1, transforming growth factor β1; Rac2, Rac family small GTPase 2; Itgb2, integrin subunit β2; Cxcr4, C-X-C motif chemokine receptor 4; Chrna7, cholinergic receptor nicotinic α7 subunit; Syn2, synapsin II; Kcnc1, potassium voltage-gated channel subfamily C member 1; Atp2b2, ATPase plasma membrane Ca2+ transporting 2; Unc13a, unc-13 homolog A; Rims1, regulating synaptic membrane exocytosis 1; Camk2g, calcium/calmodulin dependent protein kinase II γ; Camk2a, calcium/calmodulin dependent protein kinase IIα; Trhr, thyrotropin releasing hormone receptor; Grm1, glutamate metabotropic receptor 1.

Discussion

Experimental modeling of SCI in animals has been widely used to investigate the complex secondary injury cascade (1,19). The GSE5296 database consists of transcriptome data obtained from spinal cord sections from injury sites, and immediately adjacent R and C regions at 0.5, 4, 24 and 72 h, and 7 and 28 days post-injury. In the present study, bioinformatics analysis was used to determine molecular events and pathological states in SCI. Based on analysis of gene expression at different time points, DEGs were determined, and KEGG pathway and GO enrichment analyses were performed. Additionally, a PPI network and gene expression trend charts were constructed to further analyze the molecular processes underlying SCI. The analyses performed in the present study may contribute to a better understanding of SCI.

In the present study, at the M site, there were 374, 707, 1,322, 1,475, 1,724 and 1,342 DEGs identified at 0.5, 4, 24 and 72 h, and 7 and 28 days following injury, respectively. The results of the present study indicated that the majority of alterations in molecular events occurred at 7 days following injury and decreased over time. In addition, Wang et al (8) demonstrated that the number of DEGs decreased in a time-dependent manner (1,942, 396, 188 and 193 DEGs were identified at 3, 7, 14 days and 1 month, respectively). Although the R and C regions exhibited decreased fold changes in gene expression compared with the M region, further studies of these two areas may improve the understanding of the overall process of SCI.

Numerous molecular events were detected by analyzing up- and downregulated DEGs. As demonstrated in a previous study, downregulated DEGs were primarily enriched in neuronal functions (9). The spectra of expression changes in neuronal function- and synaptic transmission-associated genes (Chrna7, Syn2, Kcnc1, Atp2b2, Unc13a, Rims1, Camk2g, Camk2a, Trhr and Grm1) from 24 h to 7 days post-SCI reflects the regulation between the degeneration and survival of injured tissues. Insight into time-dependent alterations in structural and functional neuronal biomarkers may be useful for developing protective or regenerative therapeutic interventions.

As a regulator of degeneration and regeneration of the spinal cord, inflammation is a hallmark of the secondary SCI process (20). In the present study, KEGG pathway and GO enrichment analyses, and the PPI network, revealed that the upregulated DEGs were primarily associated with ‘immune system’ process and the inflammatory response, particularly at 72 h and 7 days post-SCI. As previously demonstrated, there is an association between the severity of SCI and the intensity of the acute inflammatory response, which includes proinflammatory cytokines and immune cells (20,21). A significant increase in the expression of proinflammatory (IL-1b, IL-6 and IL-7) and anti-inflammatory cytokines (IL-4 and IL-10) in the present study reflects both the regulation between degeneration and survival of injured tissues. A protective strategy is to target the process of inflammation and limit the infiltration of immune cells into the injury site (22–24). Notably, another group of inflammation-associated genes, including Cd44, Cybb, Icam1, Cyba, Hck, Casp1, Tgfb1, Rac2, Itgb2 and Cxcr4 exhibited a different temporal pattern compared with proinflammatory genes (IL-1β, IL-6 and IL-7), indicating a complex inflammatory immune microenvironment at the damaged site that requires further analysis. These findings the importance of monitoring inflammation over time following SCI (20).

In the present study, downregulated DEGs at day 7 were primarily enriched in the ‘steroid biosynthesis’ process. Steroids may be functionally divided into cholesterol, corticosteroids, sex steroids, neuroactive steroids and vitamin D (25,26). A number of studies have investigated the association between steroid metabolism and SCI (27–34). Estrogen may attenuate inflammation and promote neural survival and regeneration following SCI (27,28). Statins, known as cholesterol-controlling drugs, may significantly enhance neuronal and oligodendrocyte survival, in addition to decreasing the levels of proinflammatory cytokines (29). Previous studies additionally suggested that individuals with SCI are at an increased risk of vitamin D deficiency (30–32). Furthermore, neuroactive steroids are naturally occurring steroids that impact behavior, alter the excitability of neurons and interact with specific neurotransmitter receptors (33,34). Therefore, targeting steroid biosynthesis as a therapeutic approach for neuroprotective and neuroregenerative purposes merits further investigation.

However, one limitation of the present study was that the raw data did not include gene expression data from samples at 3 or 6 months following injury, limiting the information regarding molecular processes during the progression of SCI. KEGG pathway and GO enrichment analyses of data for the very acute phase (0.5 and 4 h post-SCI) or at 1 month post-SCI were not included in the present study. Additionally, future comprehensive analysis of transcriptome data from the adjacent R and C regions at each time point post-SCI may reveal the complex alterations that occur during the pathophysiological process.

In conclusion, the present study revealed that inflammatory response, immune processes, neuronal-associated functions and ‘steroid biosynthesis’ serve roles in the progression of SCI. Furthermore, the M region exhibited increased fold-changes in the expression of genes associated with inflammatory responses and neuronal function compared with the R and C regions at different time-points post-SCI. However, in vivo and in vitro studies are required to determine the specific roles of these molecular events in SCI.

Acknowledgements

We would like to thank Dr Yu-qing Jiang (Department of Orthopedics, the Affiliated Hospital of Nanjing Medical University, Changzhou No. 2 People's Hospital) for assistance with the present study.

Funding

The present study was supported by grants from the National Natural Science Foundation of China (grant nos. 8177235, 81520108018 and 81472080), the Jiangsu Committee of Science and Technology-Social Development Plan (grant no. BE2017755) and the Nanjing Committee of Science and Technology (grant no. 201505005).

Availability of data and materials

All data generated or analysed during this study are included in this published article.

Authors' contributions

GYY designed the present study. SJZ, WZ, JC and YJL performed the data analysis and statistical analysis. SJZ and WZ wrote and revised the manuscript. GYY supervised the present study. All authors read and approved the final manuscript.

Ethics approval and consent to participate

Not applicable.

Consent for publication

Not applicable.

Competing interests

The authors declare that they have no competing interests.

References