Spinal cord injury (SCI) results in a devastating loss of motor and sensory functions. In the United States, it has been reported that some 300,000 people are living with SCI and nearly 12,000 new patients are diagnosed annually [1]. In China, the incidence of traumatic SCI is approximately 60,000 per year [2]. Despite advances in surgical techniques for the spinal cord, there are still no effective treatments for this devastating neurological disorder. Therefore, it is very important to further understand the molecular changes in SCI in order to develop a better therapeutic program.

The pathological changes following traumatic SCI can be divided into two processes: primary injury and secondary injury. Primary injury is the immediate effect of mechanical damage on the spinal cord, which directly damages various tissue components. Secondary injury is a cascade of effects triggered by the primary injury affecting multiple biological processes and producing extensive temporal changes in gene expression [3]. Considering the complexity of gene expression and signaling pathways, a global analysis is necessary to understand the molecular mechanisms and develop therapeutic strategies for SCI. During the last decade, cDNA microarray and genechip technologies have provided valuable insights into gene changes after SCI [4–6]. However, these technologies have suffered from limitations in resolution, dynamic range and accuracy [7]. With the advances in genome-wide transcriptome analysis, especially high-throughput RNA sequencing (RNA-Seq) technology, we now have a more useful approach to observe whole transcriptome changes [7, 8].

In this study, temporal genome-wide gene expression profiles in injured spinal cords of adult rats were examined using RNA-Seq, and the results were confirmed by real-time quantitative reverse-transcriptase polymerase chain reaction (qRT-PCR). Based on these data, the functions and pathways involved in the pathologic process of SCI in the acute, subacute, and chronic phases were characterized. The bioinformatics analysis confirmed that some of the expressed genes in our study have been shown to play important roles in SCI, thus supporting the reliability of our data. We also identified many additional expressed genes whose functions in SCI have not been well studied. These additional genes may be involved in the pathological process of SCI and deserve further study and discussion.

In this study, we modeled SCI in rats using the New York University Impactor. According to the BBB scores and our previous experience [9, 15], our model produced a moderate level of impairment. To establish a systematic global analysis of gene expression patterns in the injured spinal cord, we used RNA-Seq to characterize the temporal changes after contusive SCI in rats at acute (1 dpi), subacute (6 dpi), and chronic (28 dpi) phase timepoints. The transcriptomes were systematically characterized with the goal of identifying pathways and genes critical in SCI pathology. Before data analysis, the quality of the cDNA library was examined. Our results showed that more than 77% of clean tags from RNA-Seq data mapped uniquely to the genome and that variance in the expressed transcripts as analyzed by PCA was low. These results demonstrated that our data was of sufficient quality to be used for functional analyses.

Comparing each time point to all other time points, hundreds to thousands of differentially expressed genes were obtained in the acute, subacute, and chronic phases of SCI. Hierarchical cluster analysis and K-means clustering analysis further indicated that more than 3000 genes in the two clusters decreased post-SCI and reached the lowest level at 6 dpi. There were more 4000 genes in five other clusters showing an increasing trend. Most of these genes increased post-SCI and reached the highest level at 6 dpi. There were only 243 genes in cluster 2 which expressed at the highest level at 1 dpi, followed by a decrease at 6 dpi and 28 dpi. These results suggest that our analysis may provide new information related to gene expression profiles for the study of SCI pathological mechanisms.

To further characterize the above DEG profiles, GO enrichment analysis was performed. The analysis results showed that for statistically significant differentially expressed genes, most enrichment occurred in immune response, MHC protein complex, antigen processing and presentation, translation-related genes, structural constituent of ribosome, ion gated channel activity, small GTPase mediated signal transduction, cytokine and/or chemokine activity, among others. KEGG pathway analysis showed that the top enriched pathways included ribosome, antigen processing and presentation, retrograde endocannabinoid signaling, axon guidance, dopaminergic synapses, glutamatergic synapses and GABAergic synapses, as well as TNF, HIF-1, NF-kappa B, Toll-like receptor, NOD-like receptor, cAMP, calcium, oxytocin, Rap1, B cell receptor and chemokine signaling pathway. In an early study using a mouse SCI model [16], gene profiling data at acute (2 days) and subacute phases (7 days after injury) showed that “inflammation response”, “cell death and survival”, “nervous system development” and “neurological disease” were the top enriched functional categories. “Acute Phase Response Signaling”, “LXR/RXR Activation”, “Role of Pattern Recognition Receptors in Recognition of Bacteria and Viruses” and “Atherosclerosis Signaling” were the most enriched pathways [16]. Comparing our findings to the mouse gene profiling data, we find that although animal species and data analysis methods were different, there were still many similar findings between the rat and mouse models. For example, “inflammation response” in the mouse model is similar to “immune response” in the rat model. “Role of Pattern Recognition Receptors in Recognition of Bacteria and Viruses” in mice is related to “Toll-like receptor” and “NOD-like receptor” in rats [17, 18]. These results demonstrate that there are indeed some similarities between the mechanisms of SCI in rats and mice. However, we also found many more differences in gene expression and signaling pathways between these two species. The use of high throughput RNA-Seq, coupled with different species and different laboratory results, makes it is difficult to ascertain the causes of these differences. In fact, the pathology of SCI in rats and mice is distinct, particularly with regard to cavity formation and scar/inflammatory responses, and the rat model is more similar to human SCI [19, 20]. Thus, more thorough mechanistic research using the rat model may be a better approach for understanding human SCI.

In addition to differing species, the choice detection methods may also produce important differences. There are at least four other groups that have studied gene expression during acute, subacute, and chronic SCI using methods different from RNA-Seq [21–24]. In the following sections, we will discuss the important similarities and differences in gene expression patterns in the injured spinal cord between our RNA-Seq results and data available from the other labs, as well as the advantages of RNA-Seq and the novel pathways have been discovered.

The first point we would like to make is that our results implicate many important genes previously known to be involved in SCI, such as Socs3 [25], Il17 [26], Tnf [27], Il6 [26], Il1b [28], CD44 [29], Fgf2 [30], Annexin [31], Mmp9 [32] and Bax [33]. In addition to these well-studied genes, many other genes not previously identified were found to be involved in the response to SCI. For example, in the acute phase, among the top ten significantly upregulated genes, we found that colony stimulating factor 2 receptor beta common subunit (Csf2rb), perilipin 2 (Plin2), growth arrest and DNA-damage-inducible protein GADD45 gamma (Gadd45g), strawberry notch homolog 2 (Sbno2), filamin-C (Flnc), B cell lymphoma 3 (Bcl3), glioma pathogenesis-related protein 2 (Glipr2) and tribbles pseudokinase 1 (Trib1) had not been reported in SCI.

Csf2rb is a common subunit of the GM-CSF receptor, IL-3 receptor, and IL-5 receptor [34]. GM-CSF has been shown to inhibit glial scar formation, enhance the integrity of axonal structure and produce a long-term protective effect after SCI [35]. The roles of IL-3 and IL-5 in SCI are not very clear. However, both of them and GM-CSF often co-express in TH2 cells, a subset of CD4⁺ T cells that are characterized by the production of IL-4, IL-5, IL-10 and IL-13 [36]. Therefore, in the acute phase of SCI, upregulation of their co-receptor subunit expression may be beneficial by exerting anti-inflammatory effects.

Trib1 has been shown to interact with a number of proteins, such as MEK-1, MKK4, COP1 and C/EBPalpha [37, 38]. Satoh et al. demonstrated that Trib1 plays a critical role in the differentiation of tissue-resident M2-like macrophages [39]. In our previous study, we demonstrated that macrophages can be polarized into M2 subtypes in the early stage of SCI [40]. This is consistent with the high expression of Trib1, which may play important roles in the differentiation of macrophages after SCI.

Gadd45g is involved in DNA repair, cell cycle control, cellular stress response and apoptosis, among other functions [41]. Sbno2, a novel inflammatory response factor, is predominantly expressed in astrocytes, as well as in the choroid plexus and in some microglia, endothelial cells, and neurons in the CNS [42]. Plin2, also known as adipose differentiation-related protein, is associated with adipocyte differentiation [43]. Flnc, also known as filamin-2 or actin-binding-like protein, is mainly expressed in cardiac and skeletal muscles and functions by crosslinking actin filaments and anchoring membrane proteins to the actin cytoskeleton [44]. Bcl3 is a transcriptional co-activator that is activated through interaction with p50 NF-kappaB homodimers [45]. Knockout of Bcl3 inhibited fiber atrophy and abolished NF-kappa B reporter activity [46]. Glipr2, also known as Golgi-associated plant pathogenesis-related protein 1, is a member of plant pathogenesis related proteins group 1 superfamily. It strongly associates with lipid rafts in the cytosolic leaflet of Golgi membranes [47]. However, little is known about the roles of any of these genes in SCI.

In the subacute stage, from the ten most significantly upregulated genes, we found that urokinase-type plasminogen activator (Plau), hexokinase-3 (Hk3) and triggering receptor expressed on myeloid cells 2 (Trem2) had not previously been reported in SCI. Plau is a serine protease involved in migration and proliferation of some tumor cells and degradation of the extracellular matrix [48]. Hk3, an isoform of hexokinases which catalyzes the first step of glucose metabolism, can increase ATP levels, reduce oxidant-induced reactive oxygen species, preserve mitochondrial membrane potential, increase mitochondrial biogenesis, and exert protective effects against oxidative stress [49]. Trem2, a receptor belonging to the immunoglobulin superfamily, is mainly expressed on myeloid cells where it stimulates neutrophil- and monocyte/macrophage-mediated inflammatory responses [50]. A study by Takahashi et al. reported that expression of trem2 on microglia increased their scavenging capability for apoptotic neurons without inflammation [51]. However, the functions of these genes are poorly understood in traumatic CNS injuries.

In the chronic phase, from the ten most upregulated genes, we found that tumor necrosis factor alpha-induced protein 6 (Tnfaip6), neutrophil cytosol factor 1 (Ncf1), ankyrin repeat and suppressor of cytokine signaling box protein 15 (Asb15), cupredoxins (Cp), tartrate-resistant acid phosphatase 5 (Acp5) and C-type lectin domain family 4, member A3 (Clec4a3) had not previously been reported in SCI.

Ncf1 is a cytosolic subunit of neutrophil NADPH oxidase which can be activated to produce superoxide anion. Ncf1 and NADPH oxidase 2 complex-derived reactive oxygen species are important regulators of several chronic inflammatory disorders, such as multiple sclerosis, gout, psoriasis and psoriatic arthritis, rheumatoid arthritis, and lupus [52]. Although there are no reports of increasing Ncf1 expression in SCI, it may play an important role in the chronic phase of this pathological process.

Tnfaip6, a secretory protein induced by tumor necrosis factor alpha, is involved in extracellular matrix stability and cell migration, and is important in the protease network associated with inflammation [53]. Asb15, a protein predominantly expressed in skeletal muscle, has a role in the regulation of protein turnover and muscle cell development [54]. Cps are widespread copper-binding proteins that bind a mononuclear type 1 copper redox site and are mainly involved in electron transfer reactions, photosynthesis and respiration [55]. Acp5 is involved in osteopontin/bone sialoprotein dephosphorylation, which is essential for bone resorption and osteoclast differentiation [56]. Clec4a3, a member of the C-type lectin-like domain-containing family which was originally thought to be important for initiating and shaping immune responses, has been found to be upregulated in rat spinal cord following traumatic nerve root injury [57]. However, whether these genes contribute to the pathological process of SCI remains unknown.

A second important point regarding our study is that the genes identified in SCI are involved in many signaling pathways. For example, the TNF [27, 58, 59], Toll-like receptor [60], NF-kappa B [61], NOD-like receptor [62] and synapse formation signaling pathways [63] have been reported to be involved in the pathological process of SCI [63, 64].

However, there were also some interesting new terms in SCI found by KEGG pathway analysis. For example, the ribosome term was significantly enriched in the whole pathological process. Ribosomes consist of complexes of proteins and RNAs and are complex molecular machines responsible for protein synthesis in living cells [65]. Therefore, the upregulation of ribosomal protein genes may increase production of corresponding proteins and thus boost ribosome function after SCI. The boosting function of ribosomal proteins can promote ribosome biogenesis, which may contribute to cell growth and proliferation, stress responses and the potential for repair [66]. However, the mechanisms related to ribosomes in regulating neuronal growth, stress responses and rehabilitation after SCI remain to be identified. In addition, some pathways such as retrograde endocannabinoid signaling, oxytocin signaling, Rap1 signaling and the synaptic vesicle cycle were also enriched in the subacute and chronic phases of SCI. The roles of these pathways in SCI have never been reported to our knowledge. These new findings will thus supportfuture research and therapeutic development for SCI.

In summary, our study has not only characterized global changes in gene expression in injured rat spinal cords, but has also systematically identified the critical genes and signaling pathways in SCI pathology. Although the genes associated with damage and repair of spinal cord injury are still largely unknown, the RNA-Seq analysis presented in this study will expand our understanding of the complex molecular mechanisms involved and provide a foundation for future studies of tissue damage and repair after SCI.