This article is mentioned in:

Introduction

Spinal cord injury (SCI) is a specific type of damage to the central nervous system, resulting in temporary or permanent changes in the function of the spinal cord. SCI is usually associated with traffic accidents (38%) (1), falls (31%) (2) and sports-associated injuries (10–17%) (3). SCI presents an increasing social and economic burden through its treatment and rehabilitation costs (4). Although there are a number of preclinical and clinical studies investigating this disease, its underlying molecular mechanism remains unclear.

Characterization of the molecular mechanisms of neuroplasticity following traumatic SCI may provide insight for the identification of effective treatments for the disease. Scarisbrick et al (5) indicated that kallikrein-6, a member of the kallikrein protease family, affects neural repair and regeneration in traumatic SCI. Tissue plasminogen activator promotes endogenous type 4 disintegrin and metalloproteinase with thrombospondin motifs-induced chondroitin sulfate proteoglycans degradation, advancing neuroplasticity subsequent to SCI (6). The selective inhibition of acid-sensing ion channel 1a provides morphological and functional neuroprotection following traumatic SCI (7). In addition, several studies have indicated that the abnormal expression of microRNAs (miRNAs) may be associated with SCI progression, and may be potential targets for the treatment of this disease (8,9). Despite these previous studies, the molecular mechanisms of neuroplasticity following traumatic SCI are also unclear.

Using the GSE52763 microarray dataset, Shin et al (10) identified that a number of inflammation-associated genes are upregulated in lumbar spinal cord following traumatic SCI, and treadmill locomotor training may partly improve locomotor function. Yang et al (11) described that transforming growth factor-β-induced factor homeobox 1, Ras-related C3 botulinum toxin substrate 2, TYRO protein tyrosine kinase binding protein (TYROBP), and progesterone receptor (PGR) are associated with traumatic SCI. Liu et al (12) suggested that ATPase sarcoplasmic/endoplasmic reticulum Ca2+ transporting 1, Fos Proto-Oncogene, AP-1 transcription factor subunit and glycogen synthase kinase 3β are involved in treadmill locomotor training in locomotor recovery. However, the GSE52763 microarray dataset has not been analyzed using comprehensive bioinformatics methods to reveal the mechanisms of neuroplasticity following SCI. In the present study, the GSE52763 microarray dataset was examined with multiple bioinformatics analyses, including differentially expressed genes (DEGs) screening, time series analysis, enrichment analysis, protein-protein interaction (PPI) network analysis and regulatory network analysis.

Materials and methods

Data source and data preprocessing

The GSE52763 microarray dataset, which was deposited by Shin et al (10), was downloaded from the Gene Expression Omnibus (http://www.ncbi.nlm.nih.gov/geo/) database (13). From GSE52763, 15lumbar spinal cord samples [3 sham (1 week), 4 injury only (1 week), 4 injury only (3 weeks), 4 injury + treadmill (3 weeks)] were selected. Fragmentation was assessed with an Agilent 2100 bioanalyzer using RNA 6000 Nano Chip (Agilent Technologies, Inc., Santa Clara, CA, USA) and hybridized to the arrays containing >22,500 probe sets. The sequencing data platform was [Rat230_2] Affymetrix Rat Genome 230 2.0 Array, and the data was downloaded in March 2018. The microarray dataset was downloaded from a public database; therefore, ethical approval was not obtained.

Data preprocessing (background correction, normalization and expression calculation) was performed using a robust multi-array average method in the affy R package (version 1.38.0; http://bioconductor.org/packages/release/bioc/html/oligo.html) (14). For single genes mapped to several probes, the mean value of the probes was used to represent the unique expression value of this gene.

DEGs analysis

Samples were grouped according to whether injury or motor rehabilitation was done or not. Empirical Bayes linear regression model in the limma R package (Version 3.32.5; http://bioconductor.org/packages/release/bioc/html/limma.html) (15) was used to identify DEGs in injury only (1 week) vs. sham (1 week) and injury + treadmill (3 weeks) vs. injury only (3 weeks) comparison groups, and the P-values of all genes were obtained. P<0.05 and |log2 fold change (FC)|>1 were used as the thresholds.

Time series analysis

DEGs in injury vs. sham groups and treadmill vs. non-treadmill groups were merged and considered DEGs in the injury/rehabilitation process. Genes that exhibited increased/decreased expression in the injury groups compared with the sham group, and decreased/increased expression in the injury + treadmill group compared with the injury groups were identified as significantly altered genes (candidate gene sets) following treadmill rehabilitation training.

The heatmap for the candidate gene set was drawn using R package pheatmap (16) (Version 1.0.8; http://cran.r-project.org/web/packages/pheatmap). Clustering distance was determined by Pearson correlation, and the clustering method was single-linkage clustering.

Functional and pathway enrichment analysis

Using the R package clusterProfiler (17) (Version 3.4.4; http://bioconductor.org/packages/release/bioc/html/clusterProfiler.html), gene ontology (GO)_‘Biological Process’ (BP) (18,19) and Kyoto Encyclopedia of Genes and Genomes (KEGG) (20) pathway enrichment analyses were performed for the two groups of candidate genes. P<0.05 was selected as the threshold.

PPI network analysis

PPI analysis for candidate DEGs was performed by STRING database (21) (Version 10.0; http://www.string-db.org/). The network was visualized by Cytoscape software (22) (Version 3.4.0; http://www.cytoscape.org/). CytoNCA plug-in (23) (Version 2.1.6; http://apps.cytoscape.org/apps/cytonca) was used to analyze the network topology properties of nodes. The degree of each node was calculated and nodes with the highest degrees were determined as significant nodes (hub proteins) in the PPI network (24).

Construction of regulatory network

Based on WebGestalt (25) (http://www.webgestalt.org/) tool, Over-representation Enrichment Analysis was performed to predict miRNA-target and transcription factor (TF)-target pairs for candidate DEGs. P<0.05 was selected as the threshold. The TF-miRNA-target loops were obtained by integrating the results of miRNA-target and TF-target predictions. The TF-miRNA-target network was constructed using Cytoscape software as aforementioned (21). In addition, key network nodes were obtained by performing network topology property analysis.

Results

DEGs analysis

The microarray dataset included 31,099 probes, and the expression values of 14,404 genes were obtained following gene annotation. Subsequent to pre-processing, 159 (96 upregulated and 63 downregulated) DEGs and 105 (14 upregulated and 91 downregulated) DEGs were obtained in the injury vs. sham groups and treadmill vs. non-treadmill groups, respectively. The union of these gene sets included 242 genes.

Time series analysis

A total of 40 genes in cluster 1 presented upregulated expression in the injury (1 week/3 weeks) groups compared with the sham group, and subsequent downregulation in the injury + treadmill group compared with the injury only groups. Conversely, a total of 52 genes in cluster 2 exhibited opposing expression profiles (downregulation following injury and subsequent upregulation following treadmill rehabilitation). The heatmap and line charts of these genes are presented in Fig. 1.

Figure 1. Heatmap and line graphs for candidate genes. (A) Heatmap for candidate genes: Blue represents genes with decreased expression levels; and orange-red represents genes with increased expression levels. (B) Line graphs of cluster 1, indicating increased expression levels following injury that was reversed following treadmill rehabilitation. (C) Line graph of cluster 2, indicating decreased expression levels following injury that was reversed following treadmill rehabilitation.

Functional and pathway enrichment analysis

The genes in cluster 1 were significantly enriched in certain GO_BP processes, including adaptive immune response, activation of immune response, peripheral nervous system axon regeneration and leukocyte chemotaxis (Table I), and specific pathways including Staphylococcus aureus infection, natural killer cell mediated cytotoxicity, and complement and coagulation cascades (Fig. 2A).

Figure 2. KEGG pathways enriched for (A) cluster 1 and (B) cluster 2. The horizontal axes represent the number of enriched genes, and the vertical axes represent the names of KEGG pathways. The darker colors represent decreased P-values. KEGG, Kyoto Encyclopedia of Genes and Genomes.

Table I. Top 10 significantly enriched GO-‘Biological Process’ terms for clusters 1 and 2 for differentially expressed genes.

Table I.

Top 10 significantly enriched GO-‘Biological Process’ terms for clusters 1 and 2 for differentially expressed genes.

Terms	Description	Gene symbol	Count	P-value
Cluster 1
GO:0002250	Adaptive immune response	ADGRE1, C1QA, CD48, FCGR2B, RSAD2	5	2.05×10−4
GO:0002253	Activation of immune response	C1QA, CFH, CLEC7A, LGALS3, RSAD2	5	2.33×10−4
GO:0014012	Peripheral nervous system axon regeneration	TNC, TSPO	2	2.61×10−4
GO:0071294	Cellular response to zinc ion	MT2A, TSPO	2	2.61×10−4
GO:0030595	Leukocyte chemotaxis	CXCL13, ITGB2, LGALS3, PF4	4	3.52×10−4
GO:0030593	Neutrophil chemotaxis	ITGB2, LGALS3, PF4	3	4.58×10−4
GO:0001818	Negative regulation of cytokine production	CIDEA, FCGR2B, SUZ12, TSPO	4	6.61×10−4
GO:1990266	Neutrophil migration	ITGB2, LGALS3, PF4	3	6.87×10−4
GO:0007229	Integrin-mediated signaling pathway	ITGAL, ITGB2, TYROBP	3	7.09×10−4
GO:0002366	Leukocyte activation involved in immune response	CLEC7A, ITGAL, LGALS3, TYROBP	4	7.70×10−4
Cluster 2
GO:1902803	Regulation of synaptic vesicle transport	RIMS1, STXBP1, SYT11	3	8.54×10−5
GO:0051650	Establishment of vesicle localization	RASGRP1, RIMS1, SH3GL2, STXBP1, SYT11	5	1.44×10−4
GO:0045055	Regulated exocytosis	PI4K2A, RASGRP1, RIMS1, STXBP1, SYT11	5	1.47×10−4
GO:0010765	Positive regulation of sodium ion transport	AHCYL1, ATP1B2, CNTN1	3	1.55×10−4
GO:0006836	Neurotransmitter transport	KCNJ10, RIMS1, SLC6A11, STXBP1, SYT11	5	1.59×10−4
GO:0048167	Regulation of synaptic plasticity	CAMK2N2, KCNJ10, MAP−1B, RIMS1, STXBP1	5	1.59×10−4
GO:0001505	Regulation of neurotransmitter levels	GAD2, KCNJ10, RIMS1, STXBP1, SYT11	5	1.67×10−4
GO:0051648	Vesicle localization	RASGRP1, RIMS1, SH3GL2, STXBP1, SYT11	5	2.31×10−4
GO:0099504	Synaptic vesicle cycle	RIMS1, SH3GL2, STXBP1, SYT11	4	3.93×10−4
GO:0048489	Synaptic vesicle transport	RIMS1, SH3GL2, STXBP1, SYT11	4	4.54×10−4

[i] GO, Gene Ontology.

Genes in cluster 2 were significantly enriched in several GO_BP processes, including regulation of synaptic vesicle transport, establishment of vesicle localization, regulated exocytosis and positive regulation of sodium ion transport (Table I), and certain pathways including GABAergic synapse, glucagon signaling pathway, and endocrine and other factor-regulated calcium reabsorption (Fig. 2B).

These results suggested that the majority of the genes in cluster 1 were associated with immune response, and the genes in cluster 2 were associated with signal transduction.

PPI network analysis

The PPI network for candidate DEGs was built, which included a total of 42 nodes and 44 edges (Fig. 3). Among these network nodes, 19 belonged to cluster 1, and 23 belonged to cluster 2. The 3 nodes with the highest degrees were TYROBP (degree=6), CD68 (degree=6), and integrin subunit beta 2 (degree=6), suggesting that these may be hub proteins in this PPI network.

Figure 3. PPI network. Green nodes represent cluster 1, and purple nodes represent cluster 2. The lines between the nodes represent the PPIs. PPI, Protein-protein interaction.

Construction of regulatory network

A total of 7 miRNAs and 12 TFs were included in the TF-miRNA-target network (Fig. 4). There were 21 miRNA-mRNA and 34 TF-mRNA regulatory pairs in the regulatory network. The results of the topological property analysis demonstrated that microtubule associated protein 1B (MAP-1B; degree=11), phosphofurin acidic cluster sorting protein 2 (PACS-2; degree=6), adenosylhomocysteinase-like 1 (AHCYL1; degree=6) and protein phosphatase 2 scaffold subunit alpha (PPP2R1A; degree=6) were regulated by more miRNAs and TFs in comparison with other genes in the regulatory network. PGR, early growth responsive, androgen receptor and tyrosine aminotransferase regulated several of these 4 genes.

Figure 4. TF-miRNA-target regulatory network. The orange triangles represent TF, and the blue inverted triangles represent miRNAs. The green diamonds denote cluster 1, and the purple ovals denote cluster 2. The arrowheads represent the TF-mRNA regulatory associations, and the T arrows represent the miRNA-mRNA regulatory associations. TF, Transcription factor; miRNA, microRNA.

Discussion

In the present study, 159 (96 upregulated and 63 downregulated) DEGs and 105 (14 upregulated and 91 downregulated) DEGs were obtained in the injury vs. sham groups and treadmill vs. non-treadmill groups, respectively. Additionally, 40 genes in cluster 1 were upregulated in the injury (1 week/3 weeks) groups compared with the sham group, and downregulated in the injury + treadmill group compared with the injury only groups. In cluster 2, 52 genes were downregulated in the injury (1 week/3 weeks) groups compared with the sham group, and subsequently upregulated in the injury + treadmill group compared with the injury only groups. The results from the enrichment analysis suggested that genes in clusters 1 and 2 were enriched in immune response and signal transduction, respectively. In addition, a PPI network was built for the candidate DEGs, which involved 19 genes in cluster 1 and 23 genes in cluster 2. In addition, MAP-1B, PACS-2, and AHCYL1 exhibited higher degrees in the regulatory network, and were regulated by miRNAs including miR-34A, miR-34B, miR-34C and miR-449.

MAP-1B has been demonstrated to serve roles in the progression of the nervous system (26). Ma et al (27) also indicated that the differential regulation of MAP-1B is important for development of the central nervous system. MAP-1B is involved in neuronal migration, neuronal differentiation and axonal regeneration (28), and is required for the development of the dendritic spine and maturation of synapses (29). Therefore, MAP-1B may be involved in the development of the pathogenesis of traumatic SCI by affecting synaptic maturation and dendritic spine development, and additionally affecting neuronal migration, neuronal differentiation and axonal regeneration.

Furthermore, Köttgen et al (30) suggested that the PACS proteins may be associated with ion channel trafficking. Ion channel blockers may have potential roles in SCI progression (31). Kawaai et al (32) hypothesized that AHCYL1, also known as IP(3)Rs binding protein released with IP(3) 2 (IRBIT2), contributed to neuronal function and interacted with synaptic molecules, and demonstrated that mice lacking IRBIT2 exhibited an increased locomotor activity. Therefore, PACS-2 and AHCYL1 may be associated with traumatic SCI.

In the non-proliferative stage, the upregulated expression of the miR-34 family is involved in maintaining mature neurons, and miR-34A serves a significant role in neuronal differentiation through arresting cells in G1 phase (33). Aranha et al (34) indicated that miR-34A regulated neural stem cell differentiation in mice. Rokavec et al (35) suggested that the miR-34 family served an important function in neuronal development. Therefore, miR-34A/B/C may be essential in the progression of traumatic SCI by affecting neuronal development.

Zhu et al (36) revealed that electro-acupuncture promoted neural stem cells proliferation and neuron survival via downregulation of miR-449a in rats with SCI. Furthermore, miR-449a is involved in the regulation of autophagy (37), which serves a role in SCI (38). Administration of rosiglitazone may decrease autophagy and promote recovery in experimental traumatic SCI (39). Therefore, miR-449 may be involved in traumatic SCI development. In the present study, MAP-1B, PACS-2 and AHCYL1 exhibited the highest degrees and were regulated by miRNAs including miR-34A, miR-34B, miR-34C and miR-449 in the regulatory network. In light of the aforementioned data, we hypothesized that miR-34A/B/C and miR-449 served roles in the development of traumatic SCI, partly by targeting MAP-1B, PACS-2 and AHCYL1.

The present study explored the mechanisms of neuroplasticity following SCI using comprehensive bioinformatics methods. However, only the genes of rat samples were analyzed, and therefore the genes and results described do not directly apply to humans. Additionally, the lack of in vivo and in vitro experiments was also a major limitation in the present study. Therefore, additional verification analyses are required to confirm the results obtained.

In conclusion, MAP-1B, PACS-2 and AHCYL1 are key genes for the development of traumatic SCI. Furthermore, MAP-1B, PACS-2 and AHCYL1 were regulated by miR-34A/B/C and miR-449 in the progression of traumatic SCI. These data improve the current understanding of the mechanisms of neuroplasticity following traumatic SCI, and may provide promising therapeutic targets for the disease.

Acknowledgements

Not applicable.

Funding

No funding was received.

Availability of data and materials

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

Authors' contributions

HC and YZ were responsible for the conception and design of the research, and drafting the manuscript. ZC performed the data acquisition. BZ performed the data analysis and interpretation. HW and LA participated in the design of the study and performed the statistical analysis. All authors have read and approved the manuscript.

Ethics approval and consent to participate

Not applicable.

Patient consent for publication

Not applicable.

Competing interests

The authors declare that they have no competing interests.

References