Identification of pivotal genes and pathways for spinal cord injury via bioinformatics analysis.

The present study aimed to identify key genes and pathways associated with spinal cord injury (SCI) and subsequently investigate possible therapeutic targets for the condition. The array data of GSE20907 was downloaded from the Gene Expression Omnibus database and 24 gene chips, including 3‑day, 4‑day, 1‑week, 2‑week and 1‑month post‑SCI together with control propriospinal neurons, were used for the analysis. The raw data was normalized and then the differentially expressed genes (DEGs) in the (A) 2‑week post‑SCI group vs. control group, (B) 1‑month post‑SCI group vs. control group, (C) 1‑month and 2‑week post‑SCI group vs. control group, and (D) all post‑SCI groups vs. all control groups, were analyzed with a limma package. Gene Ontology annotation and Kyoto Encyclopedia of Genes and Genomes pathway enrichment analyses for DEGs were performed. Cluster analysis was performed using ClusterOne plugins. All the DEGs identified were associated with immune and inflammatory responses. Signal transducer and activator of transcription 3 (STAT3), erb‑B2 receptor tyrosine kinase 4 (ERBB4) and cytochrome B‑245, α polypeptide (CYBA) were in the network diagrams of (A), (C) and (D), respectively. The enrichment analysis of DEGs identified in all samples demonstrated that the DEGs were also enriched in the chemokine signaling pathway (enriched in STAT3) and the high‑affinity immunoglobulin E receptor (FcεRI) signaling pathway [enriched in proto‑oncogene, src family tyrosine kinase (LYN)]. Immune and inflammatory responses serve significant roles in SCI. STAT3, ERBB4 and CYBA may be key genes associated with SCI at certain stages. Furthermore, STAT3 and LYN may be involved in the development of SCI via the chemokine and FcεRI signaling pathways, respectively.

Introduction

Spinal cord injury (SCI) is caused by incident trauma or stroke, infection, inflammation and several congenital disorders (1,2). In the UK, ~40,000 people suffer from SCI with an annual incidence of 13 per million (3). There are ~60,000 new cases of SCI reported in China and 12,000 new cases in the USA every year (4,5). SCI may cause significant mortality and morbidity (6). There is no effective therapy for SCI and it is prevalent in a number of different countries (7,8). Therefore, it is necessary to identify an effective treatment for SCI, and identify specific genes or pathways that are associated with SCI and may be able to provide novel ways of treating this condition.

A number of achievements associated with molecular mechanism of SCI have been made in previous years. The pathobiology of SCI comprises a primary mechanical insult to the spinal cord and a delayed secondary injury; the secondary injury may instigate apoptotic and necrotic cell death mechanisms (2). Neuroinflammation is a significant component for secondary injury in SCI (9) and the expression of apoptotic genes is increased, whereas the expression of anti-apoptotic genes is reduced (10). Axonal regeneration and functional recovery following SCI may be promoted by the inhibition of ephrin type-A receptor 4 in wild-type mice (11). A temporary block of the interleukin-6 signaling pathway modifies the inflammatory response following SCI and thus promotes the regeneration of the spinal cord (12). Jin et al (13) suggested that the oxidative phosphorylation and central nervous system disease signaling pathways, in addition to the cell cycle and immune system pathways, served significant roles in the progression of SCI. Lai et al (14) indicated that immune response, olfactory transduction and cell cycle pathways were associated with SCI. However, certain key genes and pathways of SCI at different stages have not been completely studied. Further investigation is necessary to elucidate these key genes and pathways to provide possible therapeutic targets for the treatment of SCI.

Previous studies with the microarray data of GSE20907 have indicated that oligodendrocyte transcription factor 1, activating transcription factor 3, jun proto-oncogene, vascular endothelial growth factor A, cell division cycle 20, CD74 molecule, major histocompatibility complex, class II invariant chain, bromodomain containing 2 and SMAD family member 1 are key genes associated with SCI, and that SCI is associated with cardiovascular disease and cancers (13–15). The present study identified the differentially expressed genes (DEGs) at different stages of SCI based on the microarray data of GSE20907. Subsequently, functional enrichment analyses for DEGs were performed. A protein-protein interaction (PPI) network was constructed. The present study aimed to identify certain other key genes and pathways in SCI at different stages and investigate potential therapeutic targets of the disease.

Materials and methods

Microarray data

The array data of GSE20907 was downloaded from the Gene Expression Omnibus database, which is based on the platform of GPL6247 (Affymetrix Rat Gene 1.0 ST Array; Affymetrix; Thermo Fisher Scientific, Inc., Waltham, MA, USA) and was deposited by Siebert et al (16). A total of 24 gene chips, comprising 3-day, 4-day, 1-week, 2-week and 1-month post-SCI together with control propriospinal neurons, were used for the analysis. The sample numbers of every time point are presented in Table I.

Table I.

Sample numbers of every time point.

Time point	Control	Post-SCI
3 days	0	4
4 days	4	0
1 week	4	0
2 weeks	4	2
1 month	4	2

SCI, spinal cord injury.

Data preprocessing

The raw expression data were normalized using the robust multiarray average algorithm in oligo (17). Following the removal of the internal standard probe, the boxplot was performed for the expression values of every chip prior to and following normalization.

DEGs analysis

The DEGs in (A) 2-week post-SCI group vs. control group, (B) 1-month post-SCI group vs. control group, (C) 1-month and 2-week post-SCI group vs. control group, and (D) all post-SCI group vs. all control group, were analyzed with a limma package (18); |log2FC (fold-change)|≥0.5 and false discovery rate ≤0.01 were set as the threshold criteria. Venn diagrams were plotted for these four comparison groups by a Venn diagrams package (19), and the overlaps of DEGs were obtained. The heatmaps were obtained for the expression values of DEGs by Pheatmap function (20) and if the identified DEGs could be distinguished the corresponding group was verified.

Gene Ontology (GO) and pathway enrichment analyses

The database for annotation, visualization and integrated discovery (DAVID) is an integrated data-mining environment and is used for the analysis of gene lists (21). GO (geneontology.org) is a tool used to annotate genes by collecting a defined, structured and controlled vocabulary (22). The Kyoto Encyclopedia of Genes and Genomes (KEGG) is a database used to put associated gene sets into their respective pathways (23).

The DEGs in four comparison groups were submitted to DAVID. Subsequently, GO annotation and KEGG pathway enrichment analyses were performed. Gene counts ≥5 and a P≤0.01 were set as the threshold values.

PPI network analysis

The DEGs were imported into the BisoGenet plugin in Cytoscape (24), known interactions were searched with the Biomolecular Interaction Network Database (BIND) (25) and the network diagrams were plotted. At the same time, cluster analysis was performed using default parameter of the ClusterOne plugin in Cytoscape (26) with P<0.01.

Results

Data preprocessing and DEG screening

The boxplots of samples prior to and following normalization are presented in Fig. 1. Volcano plots for the groups are presented in Fig. 2. The results demonstrated that the expression values of each sample were close to the same following normalization.

Figure 1.

Boxplots of sample data (A) prior to and (B) following normalization. The lateral axis represents names of samples and the longitudinal axis represents expression levels. The horizontal line in the middle of post represents the expression levels of each sample. If the line of each sample was close to the same following normalization, then that indicated that all data were successfully normalized. Ctrl, control; post, post spinal cord injury.

Figure 2.

Volcano plot of (A) 2-week post-SCI group vs. control group, (B) 1-month post-SCI group vs. control group, (C) 1-month and 2-week post-SCI group vs. control group, and (D) all post-SCI groups vs. all control groups. SCI, spinal cord injury, FC, fold-change.

A total of 148 DEGs (31 downregulated and 117 upregulated) were obtained with |log2FC|≥0.5 in group (A). The volcano plot is presented in Fig. 2A. A heat map by clustering analysis is presented in Fig. 3; it indicated that SP110 nuclear body protein and olfactory receptor 1645 separated the 2-week post-SCI group from the control group. The top five genes with maximum change are presented in Table II.

Figure 3.

Heat maps of the differentially expressed genes in 2-week post-SCI group vs. control group. SCI, spinal cord injury; Ctrl, control; post, post-SCI.

Table II.

Top five genes with maximum change in a 2-week time point, 1-month time point, 2-week and 1-month time point, and all samples respectively.

Probe set ID	Log fold-change	P-value	Symbol
2-week
10729312	−0.94	4.21×10−4	ENSRNOT00000011249
10847189	−0.92	9.93×10−3	Olr697
10910167	−0.85	6.73×10−4	ENSRNOT00000020374
10812399	−0.78	1.00×10−3	Hapln1
10846970	−0.75	7.58×10−3	LOC100365431
10761128	1.41	2.53×10−4	Hspb1
10770710	1.83	4.50×10−5	Atf3
10750524	1.95	9.93×10−3	Mx2
10775624	2.06	2.84×10−5	ENSRNOT00000063555
10753425	2.43	5.09×10−3	Mx1
1-month
10709450	−1.63	2.35×10−3	Olr122
10840577	−1.19	1.60×10−3	Cstl1
10707992	−1.17	2.23×10−3	Klhl25
10777123	−0.93	1.97×10−3	ENSRNOT00000037316
10904593	−0.91	8.74×10−3	RGD1565410
10938893	1.50	9.33×10−4	Mir325
10735369	1.55	3.82×10−3	ENSRNOT00000052640
10814301	1.57	4.79×10−3	ENSRNOT00000053252
10887088	1.86	4.58×10−3	Mir544
10775731	2.47	8.38×10−3	Cxcl13
2-week and 1-month
10733056	−1.72	6.50×10−3	Ifi47
10775624	−1.58	9.09×10−3	ENSRNOT00000063555
10840577	−1.46	4.92×10−3	Cstl1
10909985	−1.32	2.97×10−3	ENSRNOT00000053078
10936645	−1.32	7.15×10−3	RGD1563945
10934666	1.41	7.43×10−3	ENSRNOT00000057449
10769215	1.52	3.12×10−3	ENSRNOT00000059261
10796258	1.70	8.65×10−3	ENSRNOT00000052594
10878815	1.83	5.88×10−4	LOC690079
10814301	2.25	4.23×10−3	ENSRNOT00000053252
All samples
10888931	−0.89	3.55×10−4	Abhd1
10709450	−0.72	5.53×10−3	Olr122
10878210	−0.64	1.25×10−3	Kank4
10799241	−0.62	1.87×10−3	Idi1
10708471	−0.59	8.79×10−4	Olr34
10880731	1.52	4.31×10−6	C1qc
10936482	1.59	5.38×10−4	Timp1
10770710	1.62	1.34×10−7	Atf3
10761128	1.64	1.69×10−6	Hspb1
10775731	2.07	2.12×10−4	Cxcl13

A total of 128 DEGs (24 downregulated and 104 upregulated) were obtained with |log2FC|≥0.5 in group (B). The volcano plot is presented in Fig. 2B. It was identified from the results of clustering analysis that cystatin-like 1 and ENSRNOT00000052640 separated the control group from the post-SCI group (Fig. 4). The top five genes with maximum change were presented in Table II.

Figure 4.

Heat maps of the differentially expressed genes in 1-month post-SCI group vs. control group. SCI, spinal cord injury; Ctrl, control; post, post-SCI.

A total of 133 DEGs (56 downregulated and 77 upregulated) were obtained with |log2FC|≥0.5 in group (C). The volcano plot was presented in Fig. 2C. A heat map by clustering analysis is presented in Fig. 5 and it indicated that cell death-inducing DFFA-like effector B and LOC689961 divided the samples into two groups. The top five genes with maximum change are presented in Table II.

Figure 5.

Heat maps of the differentially expressed genes in 1-month and 2-week post-SCI group vs. control group. SCI, spinal cord injury; Ctrl, control; post, post-SCI.

A total of 238 DEGs (11 downregulated and 227 upregulated) were obtained with |log2FC|≥0.5 in group (D). The volcano plot is presented in Fig. 2D. From the results of clustering analysis, it was identified that abhydrolase domain containing 1 and ENSRNOT00000053282 divided the samples into two groups (Fig. 6). The top five genes with maximum change are presented in Table II.

Figure 6.

Heat maps of the differentially expressed genes in all post-SCI groups vs. all control groups. SCI, spinal cord injury; Ctrl, control; post, post-SCI.

The Venn diagrams of the above four sets of DEGs [(A), (B), (C) and (D)] are presented in Fig. 7.

Figure 7.

Venn diagrams of four sets of differentially expressed genes. Post, post-spinal cord injury.

GO and pathway enrichment analyses

A total of 7 GO terms, 6 GO terms and 251 (9 KEGG and 242 GO) terms were obtained in DEGs of (A), (B) and (D), respectively following functional clustering analysis of the DEGs of the above four sets; there was no significant enrichment results in the DEGs of (C).

The functional enrichment analysis of group (A) demonstrated that the DEGs were primarily enriched in GO terms associated with wound and inflammatory responses, including GO: 0009611-response to wounding and GO: 0002526-acute inflammatory response (Table III). The functional enrichment analysis of group (B) demonstrated that the DEGs were primarily enriched in GO terms associated with the immune system and inflammatory response, including GO: 0006955-immune response and GO: 0006954-inflammatory response (Table III). The enrichment analysis of DEGs identified in all samples of the (D) group demonstrated that the DEGs were primarily enriched in the GO terms and pathways associated with the immune system and inflammatory response, including GO: 0006955-immune response and GO: 0006954-inflammatory response, in addition to chemokine signaling pathway and high-affinity immunoglobulin E receptor (FcεRI) signaling pathway. Signal transducer and activator of transcription 3 (STAT3) and src family tyrosine kinase (LYN) are enriched in these two signaling pathways, respectively (Table III).

Table III.

GO and KEGG pathway enrichment analysis for differentially expressed genes.

Term	Description	Counts	P-value
2-week
GO:0009611	Response to wounding	10	1.72×10−3
GO:0002526	Acute inflammatory response	5	2.82×10−3
GO:0006954	Inflammatory response	7	2.92×10−3
GO:0006952	Defense response	9	4.42×10−3
GO:0001775	Cell activation	7	5.00×10−3
GO:0034097	Response to cytokine stimulus	5	6.15×10−3
GO:0051336	Regulation of hydrolase activity	7	7.39×10−3
1-month
GO:0006955	Immune response	11	2.16×10−6
GO:0002684	Positive regulation of immune system process	8	2.13×10−5
GO:0050778	Positive regulation of immune response	6	1.69×10−4
GO:0002252	Immune effector process	5	1.04×10−3
GO:0048584	Positive regulation of response to stimulus	6	1.63×10−3
GO:0006954	Inflammatory response	5	8.46×10−3
All samples
GO:0006955	Immune response	40	1.84×10−25
GO:0006952	Defense response	37	8.08×10−23
GO:0001775	Cell activation	28	1.51×10−19
GO:0009611	Response to wounding	34	6.34×10−19
GO:0006954	Inflammatory response	25	2.35×10−17
KEGG: rno04650:	Natural killer cell mediated cytotoxicity	13	3.22×10−8
KEGG: rno05322:	Systemic lupus erythematosus	9	6.82×10−5
KEGG: rno04610:	Complement and coagulation cascades	8	9.40×10−5
KEGG: rno04662:	B cell receptor signaling pathway	8	1.46×10−4
KEGG: rno04060:	Cytokine-cytokine receptor interaction	12	2.11×10−4

GO, Gene Ontology; KEGG, Kyoto Encyclopedia of Genes and Genomes.

Following a BIND search, no known interaction among DEGs was identified in group (B). The network diagram of DEGs identified in group (A) is presented in Fig. 8A, and STAT3 and cbl proto-oncogene B, E3 ubiquitin protein ligase are in this diagram. Erb-B2 receptor tyrosine kinase 4 (ERBB4) was in the network diagram of DEGs identified in group (C; Fig. 8B). The network diagram of DEGs identified in group (D) is presented in Fig. 8C and cytochrome B-245, α polypeptide (CYBA), STAT3, protein tyrosine phosphatase, non-receptor type 6, S100 calcium binding protein A4 and Cluster of Differentiation 48 were in this diagram.

Figure 8.

Network diagram of (A) 2-week time points, (B) 2-week time and 1-month time series analysis and (C) the comparison between all control groups and all post-SCI groups.

Discussion

In the present study, with the dataset of GSE20907, 148, 128, 133 and 238 DEGs were obtained in the 2-week post-SCI group vs. control group, 1-month post-SCI group vs. control group, 1-month and 2-week post-SCI group vs. control group, and all samples, respectively. The results demonstrated that the DEGs were mainly associated with the immune system and inflammatory response. STAT3, ERBB4, and CYBA were in the network diagram of 2-week post-SCI group vs. control group, 1-month and 2-week post-SCI group vs. control group, and all samples, respectively. The enrichment analysis of DEGs identified in all samples demonstrated that the DEGs were also enriched in the chemokine and FcεRI signaling pathways, and that STAT3 and LYN are enriched in these two signaling pathways, respectively.

SCI may initiate immune responses by producing cytokines and chemokines, and by infiltration of peripheral leucocytes, and thus SCI-induced inflammation may cause a further deterioration in certain functions (27). Furthermore, it is indicated that inflammation serves a part in neuroprotection, axonal regeneration and functional recovery following SCI (28). Nuclear factor kappa-light-chain-enhancer of activated B cells (NF-κB) is a pivotal regulator of inflammation and the suppression of astroglical NF-κB is able to reduce inflammation following SCI (29). Certain other studies also identify the important parts of the immune system and inflammatory response following SCI (30–32). In the present study, the DEGs were primarily associated with immune system and inflammatory response. The results therefore were in accord with previous studies and suggested that the immune system and inflammatory response served significant roles in the development of SCI.

In the present study, ERBB4 and CYBA were in the network diagram of 1-month and 2-week post-SCI group vs. control group, and all samples, respectively. Lindholm et al (33) demonstrated that in the vicinity of the lesion, a number of scattered cells expressed ErbB mRNA following SCI. Neuregulin1/ErbB signaling may be involved in the development of SCI (34). Nakano et al (35) indicated that vestibular and immune defects may be caused by the mutation of the CYBA gene encoding p22phox. As stated above, the immune system serves an important role in the progression of SCI and thus, CYBA is associated with the progression of SCI via the immune system. Thus, the results of the present study were consistent with previous studies, and suggested that ERBB4 and CYBA may be key genes associated with SCI at a certain stage of the condition.

In the present study, the enrichment analysis of DEGs identified in all samples demonstrated that the DEGs were also enriched in the chemokine signaling pathway and that STAT3 is enriched in this signaling pathway. STAT3 was in the network diagram of the 2-week post-SCI group vs. the control group. Chemokines may attract stem cells to the injured site and promote their survival following SCI (36). Taylor et al (37) suggested that cerebrospinal fluid cytokines and chemokines are dysregulated following acute thoracolumbar SCI in dogs. Chemokine-ligands/receptors are key mediator of secondary injury response following SCI (38) and STAT3, a transcription factor, is activated by tyrosine phosphorylation in response to epidermal growth factor and interleukin-6 (39). STAT3 is able to inhibit the spread of inflammation following SCI by regulating certain aspects of reactive astrogliosis, particularly the upregulation of glial scar formation and astrocyte hypertrophy, and the scar-forming astrocytes (40). Therefore, STAT3 may be a key gene associated with SCI and may be involved in the development of SCI via a chemokine signaling pathway.

LYN, a non-receptor src-like tyrosine kinase, is a significant signaling mediator and is able to regulate apoptosis, differentiation and migration (41). Certain studies have suggested that LYN is involved in the process of oligodendrocyte differentiation (42,43) and serves a vital role in myelination in the brain (44). The number of newly-formed oligodendrocytes and the branches of the myelin membrane are decreased in FYN proto-oncogene Src family tyrosine kinase-deficient mice (45). Certain studies have suggested that the death of oligodendrocytes may cause spinal cord dysfunction and that minocycline treatment may reduce delayed oligodendrocyte death following SCI (46,47). In the present study, the enrichment analysis of DEGs identified in all samples demonstrated that the DEGs were also enriched in the FcεRI signaling pathway and that LYN is enriched in this pathway. Although the roles of the FcεRI signaling pathway have not been fully elucidated in SCI, it is hypothesized that LYN may be involved in the development of SCI via the FcεRI signaling pathway.

In conclusion, the immune system and inflammatory response may serve significant roles in the development of SCI. STAT3, ERBB4 and CYBA may be key genes associated with SCI at a certain stage of the disorder. Furthermore, STAT3 and LYN may be involved in the development of SCI, via the chemokine and FcεRI signaling pathways, respectively. However, as the sample size of the present study was small and no experimental verification was performed, further studies on the mechanisms of SCI development are required.