Skip to main content
Download PDF

Transcriptome profile of rat genes in injured spinal cord at different stages by RNA-sequencing

BMC Genomics volume 18, Article number: 173 (2017) Cite this article

Abstract

Background

Spinal cord injury (SCI) results in fatal damage and currently has no effective treatment. The pathological mechanisms of SCI remain unclear. In this study, genome-wide transcriptional profiling of spinal cord samples from injured rats at different time points after SCI was performed by RNA-Sequencing (RNA-Seq). The transcriptomes were systematically characterized to identify the critical genes and pathways that are involved in SCI pathology.

Results

RNA-Seq results were obtained from total RNA harvested from the spinal cords of sham control rats and rats in the acute, subacute, and chronic phases of SCI (1 day, 6 days and 28 days after injury, respectively; n = 3 in every group). Compared with the sham-control group, the number of differentially expressed genes was 1797 in the acute phase (1223 upregulated and 574 downregulated), 6590 in the subacute phase (3460 upregulated and 3130 downregulated), and 3499 in the chronic phase (1866 upregulated and 1633 downregulated), with an adjusted P-value <0.05 by DESeq. Gene ontology (GO) enrichment analysis showed that differentially expressed genes were most enriched 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 and cytokine and/or chemokine activity. Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis showed that the most enriched pathways included ribosome, antigen processing and presentation, retrograde endocannabinoid signaling, axon guidance, dopaminergic synapses, glutamatergic synapses, GABAergic synapses, TNF, HIF-1, Toll-like receptor, NF-kappa B, NOD-like receptor, cAMP, calcium, oxytocin, Rap1, B cell receptor and chemokine signaling pathway.

Conclusions

This study has not only characterized changes in global gene expression through various stages of SCI progression in rats, but has also systematically identified the critical genes and signaling pathways in SCI pathology. These results will expand our understanding of the complex molecular mechanisms involved in SCI and provide a foundation for future studies of spinal cord tissue damage and repair.

The sequence data from this study have been deposited into Sequence Read Archive (http://www.ncbi.nlm.nih.gov/sra; accession number PRJNA318311).

Background

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.

Methods

Animals

We used 12 female Sprague-Dawley rats (eight-week-old, 220–250 g, RRID: RGD_70508) in this study. Animal care in surgical procedures and after operation were in accordance with Regulations for the Administration of Affairs Concerning Experimental Animals (Ministry of Science and Technology, China; revised in June 2004) and Guidelines and Policies for Rodent Survival Surgery provided by the Animal Care and Use Committees of Bengbu Medical College.

Contusive SCI

To perform the contusive SCI, a New York University Impactor was used as described [9]. Rats were anesthetized by using pentobarbital (50 mg/kg, intraperitoneal) and received a T9 vertebral laminectomy. Sham-operated (sham) rats only received a laminectomy without contusive injury. For making contusive injury, the spine was stabilized by clamping T7 and T11 spinous processes. A weight drop injury of the cord was made by using a 10 g rod (2.5 mm diameter) dropped at a height of 25 mm on the exposed dorsal surface. After injury, animals were placed in a humidity- and temperature-controlled chamber. In 3 consecutive days, buprenorphine (0.05 mg/kg, SQ; Reckitt Benckise, Hull, England) was used every 12 h to relieve pain. Bladder emptying was manually performed three times daily until autonomous urination was established. Chloramphenicol (50 ~ 75 mg/kg) was daily provided via drinking water to prevent infections.

Basso, Beattie, and Bresnahan (BBB) locomotor rating scale

BBB locomotor rating scale was used to assess the hind-limb movements based on observation of the rat’s freely moving in an open field [10, 11]. The evaluations were performed by 2 scorers while rats were walking freely on the open-field surface for 4 min. The scores were evaluated before operation, at 1 day, 3, 6, 10, 14, 21, and 28 days after injury.

RNA isolation, quantification, and qualification

Spinal cords (5 mm spinal cord segment containing the injury epicenter) of sham, acute, subacute, and chronic phases were harvested at 6 days, 1 day, 6 days, and 28 days after injury, respectively (n = 3 in every group). Rats were perfused with phosphate buffer saline (PBS) before removing the spinal cords. The total RNA from spinal cords was extracted using TRIzol reagent (Invitrogen, New Jersey, NJ, USA). Genomic DNA was removed by using DNase I. The 1% agarose gels were used to monitor RNA degradation and contamination. NanoPhotometer® spectrophotometer (IMPLEN, CA, USA) was used to check RNA purity. Qubit® RNA Assay Kit in Qubit® 2.0 Flurometer (Life Technologies, CA, USA) was used to measure RNA concentration. RNA Nano 6000 Assay Kit of the Bioanalyzer 2100 system (Agilent Technologies, CA, USA) was used to assess RNA integrity.

Library preparation for transcriptome sequencing

For the sample preparations, 3 μg RNA was used as input material for every sample. NEBNext® Ultra™ RNA Library Prep Kit for Illumina® (NEB, USA) was used to generate sequencing libraries. For each sample, the index codes were added to attribute sequences. The poly-T oligo-attached magnetic beads were used to purify mRNA from total RNA. In NEBNext First Strand Synthesis Reaction Buffer (5×), the mRNA was fragmented by divalent cations under elevated temperature. M-MuLV Reverse Transcriptase (RNase H) and random hexamer primers were used to synthesize the first strand cDNA. RNase H and DNA polymerase I were subsequently used to synthesize the second strand cDNA. The exonuclease/polymerase activities converted remaining overhangs into blunt ends. To prepare for hybridization, the 3′ ends of DNA fragments were adenylated and ligated to NEBNext Adaptor with hairpin loop structure. AMPure XP system (Beckman Coulter, Beverly, USA) was used to purify the library fragments to select cDNA fragments of preferentially 150–200 bp. Before PCR, size-selected, adaptor-ligated cDNAs were treated with 3 μl USER Enzyme (NEB, USA) at 37 °C for 15 min followed by 95 °C for 5 min. Index (X) Primer, Universal PCR primers and Phusion High-Fidelity DNA polymerase were used to perform PCR. Finally, AMPure XP system was used to purify PCR products, and the Agilent Bioanalyzer 2100 system was used to assess library quality.

Clustering and sequencing

TruSeq PE Cluster Kit v3-cBot-HS (Illumina) was used to cluster the index-coded samples on a cBot Cluster Generation System by Novogene (Beijing, China) according to the manufacturer’s instructions. Finally, the 125 bp/150 bp paired-end reads were generated and sequenced on an Illumina Hiseq platform.

Data analysis

Quality control

To obtain high-quality clean data for downstream analyses, the low quality reads and reads containing adapters or ploy-N in raw data were removed and the Q20, Q30 and GC content of those clean data were also calculated.

Reads mapping to the reference genome

From genome website, the files of reference genome and gene model annotation were downloaded. Bowtie v2.2.3 was used to build the reference genome index. TopHat v2.0.12 was used to align the paired-end clean reads to the reference genome.

Gene expression level quantification

The reads numbers mapped to each gene were counted by HTSeq v0.6.1. Fragments Per Kilobase of transcript sequence per millions base pairs sequenced (FPKM) was calculated. FPKM is the most commonly used method for evaluating the sequencing depth and gene length of the reads and estimating the levels of gene expression [12].

Analysis of differential expression

Differential expression of four groups (three biological replicates per group) was analyzed using the DESeq R package (1.18.0). To control the false discovery rate, Benjamini and Hochberg’s approach was used to adjust the P-values obtained from DESeq analysis. The adjusted P-value < 0.05 was adopted as the standard for judging statistically significant differences in gene expression.

Gene ontology (GO) and Kyoto encyclopedia of genes and genomes (KEGG) enrichment analysis of differentially expressed genes (DEGs)

The GOseq R package was used to perform GO enrichment analysis of DEGs. A corrected P-value < 0.05 was adopted as the standard for judging statistically significant enrichment of differentially expressed genes. KEGG enrichment analysis of DEGs was implemented by KOBAS software.

qRT-PCR

To validate RNA-Seq results, 13 selected DEGs were verified by qRT-PCR as described [13]. A reverse transcription system (Promega, Madison, WI) was used to reverse-transcribe the total RNA used in RNA-Seq into cDNA. SYBR Green PCR Master Mix and the ABI 7900 PCR detection system (Applied Biosystems, Foster City, CA) were used to perform real-time PCR. To normalize gene expression, the housekeeping gene glyceraldehyde-3-phosphate dehydrogenase (GAPDH) was parallel amplified. PCR primer sequences are listed in Table 1. The ΔΔCt method [14] was used to calculate the relative expression level of target mRNAs. The relative mRNA expression value in sham group was designated as 1.

Table 1 Real-time PCR primers used in the study
Full size table

Results

BBB scores at different stages of SCI

To observe the level of impairment produced by SCI, BBB scores were evaluated. As shown in Fig. 1, all animals scored 21 points before injury. At 1 day post-SCI, all nine spinal cord-injured rats received a score of 0, while the three sham-operated animals still received scores of 21 points. On the 3rd and 6th days, the mean BBB scores of the six surviving injured rats were 1 ± 0 and 5.33 ± 1.53, respectively. After 6 days, with only three injured rats remaining, BBB scores had improved substantially. The mean scores at 10 days post-injury (dpi), 14 dpi, 21 dpi and 28 dpi were 8.33 ± 1.53, 9.33 ± 1.53, 10.67 ± 0.58, and 11.67 ± 0.58, respectively. According to our previous experience, this was a moderate level of impairment.

Fig. 1
figure 1

BBB scores at different stages of the SCI. The BBB scores were shown at different stages of the SCI: before injury (n = 12); sham (n = 3); 1 dpi (n = 9); 3 dpi and 6 dpi (n = 6); 10 dpi, 14 dpi, 21 dpi, and 28 dpi (n = 3)

Full size image

Identification of expressed transcripts in the rat spinal cord transcriptome

In this study, we established 12 cDNA libraries with the following designations. Libraries with the A designation were from the sham-operated spinal cords and libraries B1, B2 and B3 were from injured spinal cords at 1 dpi, 6 dpi and 28 dpi, respectively, which represented of sham and three different pathological stages of SCI (acute, subacute, and chronic phases). RNA-Seq generated 41,724,966 to 58,619,478 raw reads for each library. After filtering the low-quality reads, the average number of clean reads was 42,505,205 (87.73%), 43,520,840 (82.37%), 43,336,187 (91.37%), and 42,483,112 (89.70%) for the A, B, C and D groups, respectively (Table 2). The clean reads were used for all further analyses, and from them 77.44% to 88.79% of clean tags from the RNA-Seq data mapped uniquely to the genome, while a small proportion of them (<3%) were mapped multiple times to the genome (Table 3).

Table 2 Summary of sequence assembly after Illumina sequencing
Full size table
Table 3 Summary of clean reads mapped to the reference genome
Full size table

Identification of the source of variance in the expressed transcripts by principal component analysis (PCA)

To demonstrate the source of variance in our data, PCA analysis with three principal components (PC1, 2, and 3) was performed. As shown in Fig. 2, PC score plots showed that the contribution of PC1, 2, and 3 was 56.89%, 11.90%, and 7.77%, respectively. The three individual samples collected at the same time points after injury were clustered closely together which validated the finding of low variance in the present analysis study and showed that the data could be used for the following analysis.

Fig. 2
figure 2

PCA analysis of the expressed transcripts. PCA analysis with three principal components (PC1, 2, and 3) was performed to demonstrate the source of variance in our data (n = 3)

Full size image

Differential gene expression in spinal cords at different injury stages

RPKM was used to estimate the level of gene expression, and DEGSeq was used to examine the differential gene expression profile. Comparing the SCI and sham animals, we observed 1,979 DEGs at 1 dpi, with 1,223 upregulated and 574 downregulated genes (Fig. 3a and Additional file 1: Table S1), 6590 DEGs at 6 dpi, with 3460 upregulated and 3130 downregulated genes (Fig. 3b and Additional file 1: Table S1) and 3499 DEGs at 28 dpi, with 1866 upregulated and 1633 downregulated genes (Fig. 3c and Additional file 1: Table S1). There were 2950 DEGs between the 6 dpi and 1 dpi groups, with 1914 upregulated and 1036 downregulated genes (Fig. 3d and Additional file 1: Table S1). There were 861 DEGs between the 28 dpi and 1 dpi groups, with 591 upregulated and 270 downregulated genes (Fig. 3e and Additional file 1: Table S1). There were 2215 DEGs between the 28 dpi and 6 dpi groups, with 986 upregulated and 1229 downregulated genes (Fig. 3f and Additional file 1: Table S1). A list of the most relevant genes that were over-expressed or under-expressed and how their expression changed with time is shown in Additional file 2: Table S2.

Fig. 3
figure 3

Volcano map of differentially expressed genes. Significantly upregulated and downregulated genes are shown as a red and green dot, respectively. The blue dot represents no significant difference between the expressions of genes. a 1 dpi vs sham, b 6 dpi vs sham, c 28 dpi vs sham, d 6 dpi vs 1 dpi, e 28 dpi vs 1 dpi, f 28 dpi vs 6 dpi

Full size image

Confirmation of differential gene expression by qRT-PCR

To verify the expression profiles obtained by RNA-Seq, ten genes were randomly selected to analyze their transcript expression levels by qRT-PCR. The selected genes included C1qb, Ccl2, Cxcl2, enpp3, IL6, Lcn2, Ncf1, Pf4, Plau, and Tspo. The results showed that these genes displayed similar expression patterns between RNA-seq and qRT-PCR (Fig. 4).

Fig. 4
figure 4

Quantitative RT-PCR validations of DEGs characterized by RNA-Seq. The relative expression level of target mRNAs was calculated using the ΔΔCt method and expressed relative to the value in the sham group (designated as 1). Data represent the mean ± SD (n = 3). Log2 fold change was the ratio of average Log2 folds between samples

Full size image

Hierarchical cluster analysis of gene expression in spinal cords at different injury stages

Based on the similarity of gene expression patterns, 7632 DEGs were classified into expression cluster groups by hierarchical clustering (Fig. 5). These clusters contained genes that were up- or downregulated throughout the whole course of SCI. Based on the expression patterns of these clusters, we can found that some showed marked differences between different stages of SCI, while others showed relatively minor differences. For example, most of the genes upregulated at 1 dpi (B), 6 dpi (C) and 28 dpi (D) compared to the sham group (A) were in the upper cluster, while downregulated genes were observed in the lower cluster.

Fig. 5
figure 5

Hierarchical cluster analysis of gene expression in the spinal cords at different injured stages. Based on similarity of gene expression patterns, 7632 DEGs were classified into many expression cluster groups. The blue to red gradation represented gene expressions from down to up. a: sham control; b: 1 dpi; c: 6 dpi; d: 28 dpi

Full size image

K-means clustering of DEGs

To further investigate the biological characteristics of the 7632 DEGs, K-means clustering was performed. In the analysis, DEGs were statistically grouped into eight subclusters based on their expression patterns in spinal cords at different injury stages. Figure 6 shows the trends of distinct significantly expressed subclusters of all DEGs. There were 3082 and 400 genes in subclusters 1 and 6, respectively. The gene expression levels showed a similar trend in both clusters, which were highest in sham (A), and lowest in SCI animals at 6 dpi (C). Subclusters 3, 4, 5, 7 and 8 included 482, 1612, 1477, 274 and 62 genes, respectively. The gene expression levels showed a similar trend in these five subclusters, with lowest expression in sham (A) and highest expression in SCI animals at 6 dpi (C). There were 243 genes in subcluster 2, which expressed at the lowest level in sham (A) and at the highest at 1 dpi (C) with a progressive decrease at 6 dpi (C) and 28 dpi (D). The list of gene IDs in eight subclusters are shown in Additional file 3: Table S3.

Fig. 6
figure 6

K-means clustering of DEGs. The 7632 DEGs were statistically grouped into 8 subclusters. The trends of distinct significant expression subclusters were analyzed. A: sham control; B: 1 dpi; C: 6 dpi; D: 28 dpi. (n = 3 in every group)

Full size image

GO enrichment of DEGs

GO enrichment analysis was conducted to characterize the DEG profiles presented above. GO enrichment results for all DEGs are provided in Additional file 4: Table S4.

In acute SCI (1 dpi), 35 significant GO terms were detected in the upregulated genes (Additional file 5: Figure S1 and Additional file 6: Table S5). No significant terms were found in the downregulated genes at 1 dpi (Additional file 7: Figure S2 and Additional file 8: Table S6).

In subacute SCI (6 dpi), 33 significant GO terms were detected in the up-regulated genes (Additional file 9: Figure S3 and Additional file 6: Table S5). Among the down-regulated genes at 6 dpi, 40 significant GO terms were detected (Additional file 10: Figure S4 and Additional file 8: Table S6).

In chronic SCI (28 dpi), 31 significant GO terms were detected in the up-regulated genes (Additional file 11: Figure S5 and Additional file 6: Table S5). Among the down-regulated genes at 28 dpi, 28 significant GO terms were detected (Additional file 12: Figure S6 and Additional file 8: Table S6).

KEGG enrichment analysis

KEGG enrichment analysis was conducted to identify pathways that play important roles in the pathologic process of SCI. KEGG enrichment results of DEGs are shown in Additional file 13: Tables S7, Additional file 14: Table S8, and Additional file 15: Table S9.

In acute SCI (1 dpi), ribosome, TNF signaling pathway, proteoglycans in cancer, malaria, and staphylococcus aureus infection were enriched pathways in all DEGs; Of the 13 enriched pathways in upregulated DEGs, the top ten were ribosome, TNF signaling pathway, malaria, salmonella infection, leishmaniasis, HIF-1 signaling pathway, proteoglycans in cancer, leukocyte transendothelial migration, Chagas disease (American trypanosomiasis), and staphylococcus aureus infection; the 2 enriched pathways in downregulated DEGs were circadian entrainment and axon guidance (Additional file 16: Figure S7 and Additional file 13: Table S7).

In acute SCI (6 dpi), glutamatergic synapse, circadian entrainment, GABAergic synapse and retrograde endocannabinoid signaling were pathways enriched in all DEGs. The top ten pathways of the 27 enriched pathways in upregulated DEGs were ribosome, osteoclast differentiation, lysosome, NF-kappa B signaling pathway, TNF signaling pathway, Toll-like receptor signaling pathway, NOD-like receptor signaling pathway, Chagas disease (American trypanosomiasis), pertussis and salmonella infection. The top ten pathways of the 37 pathways enriched in downregulated DEGs were glutamatergic synapse, circadian entrainment, GABAergic synapse, retrograde endocannabinoid signaling, calcium signaling pathway, dopaminergic synapse, neuroactive ligand-receptor interaction, adrenergic signaling in cardiomyocytes, insulin secretion and nicotine addiction (Additional file 17: Figure S8 and Additional file 14: Table S8).

In chronic SCI (28 dpi), 39 pathways were enriched in all DEGs; the top ten were ribosome, circadian entrainment, glutamatergic synapse, GABAergic synapse, retrograde endocannabinoid signaling, salivary secretion, platelet activation, cAMP signaling pathway, nicotine addiction and lysosome; In the 21 pathways enriched in upregulated DEGs, the top ten were ribosome, lysosome, hematopoietic cell lineage, staphylococcus aureus infection, osteoclast differentiation, tuberculosis, leishmaniasis, Fc gamma R-mediated phagocytosis, phagosome and NF-kappa B signaling pathway. Of the 25 pathways enriched in downregulated DEGs, the top ten were circadian entrainment, cAMP signaling pathway, glutamatergic synapse, nicotine addiction, calcium signaling pathway, GABAergic synapse, retrograde endocannabinoid signaling, adrenergic signaling in cardiomyocytes, axon guidance and neuroactive ligand-receptor interaction (Additional file 18: Figure S9 and Additional file 15: Table S9).

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.

Conclusions

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.

Abbreviations

dpi:

Day postinjury

GAPDH:

Glyceraldehyde-3-phosphate dehydrogenase

GO:

Gene ontology

KEGG:

Kyoto encyclopedia of genes and genomes

qRT-PCR:

Real-time quantitative reverse transcriptase polymerase chain reaction

RNA-Seq:

RNA-sequencing

SCI:

Spinal cord injury

References

  1. Devivo MJ. Epidemiology of traumatic spinal cord injury: trends and future implications. Spinal Cord. 2012;50(5):365–72.

    Article  CAS  PubMed  Google Scholar 

  2. Qiu J. China Spinal Cord Injury Network: changes from within. Lancet Neurol. 2009;8(7):606–7.

    Article  PubMed  Google Scholar 

  3. Liu NK, Wang XF, Lu QB, Xu XM. Altered microRNA expression following traumatic spinal cord injury. Exp Neurol. 2009;219(2):424–9.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  4. Aimone JB, Leasure JL, Perreau VM, Thallmair M, Christopher Reeve Paralysis Foundation Research C. Spatial and temporal gene expression profiling of the contused rat spinal cord. Exp Neurol. 2004;189(2):204–21.

    Article  CAS  PubMed  Google Scholar 

  5. Nesic O, Svrakic NM, Xu GY, McAdoo D, Westlund KN, Hulsebosch CE, Ye Z, Galante A, Soteropoulos P, Tolias P, et al. DNA microarray analysis of the contused spinal cord: effect of NMDA receptor inhibition. J Neurosci Res. 2002;68(4):406–23.

    Article  CAS  PubMed  Google Scholar 

  6. Tachibana T, Noguchi K, Ruda MA. Analysis of gene expression following spinal cord injury in rat using complementary DNA microarray. Neurosci Lett. 2002;327(2):133–7.

    Article  CAS  PubMed  Google Scholar 

  7. Marioni JC, Mason CE, Mane SM, Stephens M, Gilad Y. RNA-seq: an assessment of technical reproducibility and comparison with gene expression arrays. Genome Res. 2008;18(9):1509–17.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  8. Gamsiz ED, Ouyang Q, Schmidt M, Nagpal S, Morrow EM. Genome-wide transcriptome analysis in murine neural retina using high-throughput RNA sequencing. Genomics. 2012;99(1):44–51.

    Article  CAS  PubMed  Google Scholar 

  9. Lu HZ, Xu L, Zou J, Wang YX, Ma ZW, Xu XM, Lu PH. Effects of autoimmunity on recovery of function in adult rats following spinal cord injury. Brain Behav Immun. 2008;22(8):1217–30.

    Article  PubMed  Google Scholar 

  10. Basso DM, Beattie MS, Bresnahan JC. A sensitive and reliable locomotor rating scale for open field testing in rats. J Neurotrauma. 1995;12(1):1–21.

    Article  CAS  PubMed  Google Scholar 

  11. Basso DM, Beattie MS, Bresnahan JC, Anderson DK, Faden AI, Gruner JA, Holford TR, Hsu CY, Noble LJ, Nockels R, et al. MASCIS evaluation of open field locomotor scores: effects of experience and teamwork on reliability. Multicenter Animal Spinal Cord Injury Study. J Neurotrauma. 1996;13(7):343–59.

    Article  CAS  PubMed  Google Scholar 

  12. Trapnell C, Williams BA, Pertea G, Mortazavi A, Kwan G, van Baren MJ, Salzberg SL, Wold BJ, Pachter L. Transcript assembly and quantification by RNA-Seq reveals unannotated transcripts and isoform switching during cell differentiation. Nat Biotechnol. 2010;28(5):511–5.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  13. Hu JG, Shen L, Wang R, Wang QY, Zhang C, Xi J, Ma SF, Zhou JS, Lu HZ. Effects of Olig2-overexpressing neural stem cells and myelin basic protein-activated T cells on recovery from spinal cord injury. Neurotherapeutics. 2012;9(2):422–45.

    Article  CAS  PubMed  Google Scholar 

  14. Pfaffl MW. A new mathematical model for relative quantification in real-time RT-PCR. Nucleic Acids Res. 2001;29(9), e45.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  15. Ma SF, Chen YJ, Zhang JX, Shen L, Wang R, Zhou JS, Hu JG, Lu HZ. Adoptive transfer of M2 macrophages promotes locomotor recovery in adult rats after spinal cord injury. Brain Behav Immun. 2015;45:157–70.

    Article  CAS  PubMed  Google Scholar 

  16. Chen K, Deng S, Lu H, Zheng Y, Yang G, Kim D, Cao Q, Wu JQ. RNA-seq characterization of spinal cord injury transcriptome in acute/subacute phases: a resource for understanding the pathology at the systems level. PLoS One. 2013;8(8), e72567.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  17. Areschoug T, Gordon S. Pattern recognition receptors and their role in innate immunity: focus on microbial protein ligands. Contrib Microbiol. 2008;15:45–60.

    Article  CAS  PubMed  Google Scholar 

  18. Kawai T, Akira S. The role of pattern-recognition receptors in innate immunity: update on Toll-like receptors. Nat Immunol. 2010;11(5):373–84.

    Article  CAS  PubMed  Google Scholar 

  19. Byrnes KR, Fricke ST, Faden AI. Neuropathological differences between rats and mice after spinal cord injury. J Magn Reson Imaging. 2010;32(4):836–46.

    Article  PubMed  PubMed Central  Google Scholar 

  20. Sroga JM, Jones TB, Kigerl KA, McGaughy VM, Popovich PG. Rats and mice exhibit distinct inflammatory reactions after spinal cord injury. J Comp Neurol. 2003;462(2):223–40.

    Article  PubMed  Google Scholar 

  21. Jin L, Wu Z, Xu W, Hu X, Zhang J, Xue Z, Cheng L. Identifying gene expression profile of spinal cord injury in rat by bioinformatics strategy. Mol Biol Rep. 2014;41(5):3169–77.

    Article  CAS  PubMed  Google Scholar 

  22. Wen T, Hou J, Wang F, Zhang Y, Zhang T, Sun T. Comparative analysis of molecular mechanism of spinal cord injury with time based on bioinformatics data. Spinal Cord. 2016;54(6):431–8.

    Article  CAS  PubMed  Google Scholar 

  23. Chamankhah M, Eftekharpour E, Karimi-Abdolrezaee S, Boutros PC, San-Marina S, Fehlings MG. Genome-wide gene expression profiling of stress response in a spinal cord clip compression injury model. BMC Genomics. 2013;14:583.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  24. Fortune RD, Grill RJ, Jr., Beeton C, Tanner M, Huq R, Loose DS. Changes in gene expression and metabolism in the testes of the rat following spinal cord injury. J Neurotrauma 2016. [Epub ahead of print].

  25. Park KW, Lin CY, Lee YS. Expression of suppressor of cytokine signaling-3 (SOCS3) and its role in neuronal death after complete spinal cord injury. Exp Neurol. 2014;261:65–75.

    Article  CAS  PubMed  Google Scholar 

  26. Zong S, Zeng G, Fang Y, Peng J, Tao Y, Li K, Zhao J. The role of IL-17 promotes spinal cord neuroinflammation via activation of the transcription factor STAT3 after spinal cord injury in the rat. Mediat Inflamm. 2014;2014:786947.

    Article  Google Scholar 

  27. Wang CX, Nuttin B, Heremans H, Dom R, Gybels J. Production of tumor necrosis factor in spinal cord following traumatic injury in rats. J Neuroimmunol. 1996;69(1–2):151–6.

    Article  CAS  PubMed  Google Scholar 

  28. Hayashi M, Ueyama T, Tamaki T, Senba E. Expression of neurotrophin and IL-1 beta mRNAs following spinal cord injury and the effects of methylprednisolone treatment. Kaibogaku zasshi J Anatomy. 1997;72(3):209–13.

    CAS  Google Scholar 

  29. Moon C, Heo S, Sim KB, Shin T. Upregulation of CD44 expression in the spinal cords of rats with clip compression injury. Neurosci Lett. 2004;367(1):133–6.

    Article  CAS  PubMed  Google Scholar 

  30. Mocchetti I, Rabin SJ, Colangelo AM, Whittemore SR, Wrathall JR. Increased basic fibroblast growth factor expression following contusive spinal cord injury. Exp Neurol. 1996;141(1):154–64.

    Article  CAS  PubMed  Google Scholar 

  31. Liu N, Han S, Lu PH, Xu XM. Upregulation of annexins I, II, and V after traumatic spinal cord injury in adult rats. J Neurosci Res. 2004;77(3):391–401.

    Article  CAS  PubMed  Google Scholar 

  32. Hansen CN, Fisher LC, Deibert RJ, Jakeman LB, Zhang H, Noble-Haeusslein L, White S, Basso DM. Elevated MMP-9 in the lumbar cord early after thoracic spinal cord injury impedes motor relearning in mice. J Neurosci. 2013;33(32):13101–11.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  33. Antonsson B, Conti F, Ciavatta A, Montessuit S, Lewis S, Martinou I, Bernasconi L, Bernard A, Mermod JJ, Mazzei G, et al. Inhibition of Bax channel-forming activity by Bcl-2. Science. 1997;277(5324):370–2.

    Article  CAS  PubMed  Google Scholar 

  34. Woodcock JM, Zacharakis B, Plaetinck G, Bagley CJ, Qiyu S, Hercus TR, Tavernier J, Lopez AF. Three residues in the common beta chain of the human GM-CSF, IL-3 and IL-5 receptors are essential for GM-CSF and IL-5 but not IL-3 high affinity binding and interact with Glu21 of GM-CSF. EMBO J. 1994;13(21):5176–85.

    CAS  PubMed  PubMed Central  Google Scholar 

  35. Huang X, Kim JM, Kong TH, Park SR, Ha Y, Kim MH, Park H, Yoon SH, Park HC, Park JO, et al. GM-CSF inhibits glial scar formation and shows long-term protective effect after spinal cord injury. J Neurol Sci. 2009;277(1-2):87–97.

    Article  CAS  PubMed  Google Scholar 

  36. van Leeuwen BH, Martinson ME, Webb GC, Young IG. Molecular organization of the cytokine gene cluster, involving the human IL-3, IL-4, IL-5, and GM-CSF genes, on human chromosome 5. Blood. 1989;73(5):1142–8.

    PubMed  Google Scholar 

  37. Kiss-Toth E, Bagstaff SM, Sung HY, Jozsa V, Dempsey C, Caunt JC, Oxley KM, Wyllie DH, Polgar T, Harte M, et al. Human tribbles, a protein family controlling mitogen-activated protein kinase cascades. J Biol Chem. 2004;279(41):42703–8.

    Article  CAS  PubMed  Google Scholar 

  38. Sung HY, Guan H, Czibula A, King AR, Eder K, Heath E, Suvarna SK, Dower SK, Wilson AG, Francis SE, et al. Human tribbles-1 controls proliferation and chemotaxis of smooth muscle cells via MAPK signaling pathways. J Biol Chem. 2007;282(25):18379–87.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  39. Satoh T, Kidoya H, Naito H, Yamamoto M, Takemura N, Nakagawa K, Yoshioka Y, Morii E, Takakura N, Takeuchi O, et al. Critical role of Trib1 in differentiation of tissue-resident M2-like macrophages. Nature. 2013;495(7442):524–8.

    Article  CAS  PubMed  Google Scholar 

  40. Chen YJ, Zhu H, Zhang N, Shen L, Wang R, Zhou JS, Hu JG, Lu HZ. Temporal kinetics of macrophage polarization in the injured rat spinal cord. J Neurosci Res. 2015;93(10):1526–33.

    Article  CAS  PubMed  Google Scholar 

  41. Liebermann DA, Hoffman B. Gadd45 in the response of hematopoietic cells to genotoxic stress. Blood Cells Mol Dis. 2007;39(3):329–35.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  42. Grill M, Syme TE, Nocon AL, Lu AZ, Hancock D, Rose-John S, Campbell IL. Strawberry notch homolog 2 is a novel inflammatory response factor predominantly but not exclusively expressed by astrocytes in the central nervous system. Glia. 2015;63(10):1738–52.

    Article  PubMed  PubMed Central  Google Scholar 

  43. Bosma M, Hesselink MK, Sparks LM, Timmers S, Ferraz MJ, Mattijssen F, van Beurden D, Schaart G, de Baets MH, Verheyen FK, et al. Perilipin 2 improves insulin sensitivity in skeletal muscle despite elevated intramuscular lipid levels. Diabetes. 2012;61(11):2679–90.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  44. van der Flier A, Sonnenberg A. Structural and functional aspects of filamins. Biochim Biophys Acta. 2001;1538(2–3):99–117.

    Article  PubMed  Google Scholar 

  45. Richard M, Louahed J, Demoulin JB, Renauld JC. Interleukin-9 regulates NF-kappaB activity through BCL3 gene induction. Blood. 1999;93(12):4318–27.

    CAS  PubMed  Google Scholar 

  46. Hunter RB, Kandarian SC. Disruption of either the Nfkb1 or the Bcl3 gene inhibits skeletal muscle atrophy. J Clin Invest. 2004;114(10):1504–11.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  47. van Galen J, Olrichs NK, Schouten A, Serrano RL, Nolte-‘t Hoen EN, Eerland R, Kaloyanova D, Gros P, Helms JB. Interaction of GAPR-1 with lipid bilayers is regulated by alternative homodimerization. Biochim Biophys Acta. 2012;1818(9):2175–83.

    Article  PubMed  Google Scholar 

  48. Liotta LA, Goldfarb RH, Brundage R, Siegal GP, Terranova V, Garbisa S. Effect of plasminogen activator (urokinase), plasmin, and thrombin on glycoprotein and collagenous components of basement membrane. Cancer Res. 1981;41(11 Pt 1):4629–36.

    CAS  PubMed  Google Scholar 

  49. Wyatt E, Wu R, Rabeh W, Park HW, Ghanefar M, Ardehali H. Regulation and cytoprotective role of hexokinase III. PLoS One. 2010;5(11), e13823.

    Article  PubMed  PubMed Central  Google Scholar 

  50. Sessa G, Podini P, Mariani M, Meroni A, Spreafico R, Sinigaglia F, Colonna M, Panina P, Meldolesi J. Distribution and signaling of TREM2/DAP12, the receptor system mutated in human polycystic lipomembraneous osteodysplasia with sclerosing leukoencephalopathy dementia. Eur J Neurosci. 2004;20(10):2617–28.

    Article  PubMed  Google Scholar 

  51. Takahashi K, Rochford CD, Neumann H. Clearance of apoptotic neurons without inflammation by microglial triggering receptor expressed on myeloid cells-2. J Exp Med. 2005;201(4):647–57.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  52. Holmdahl R, Sareila O, Olsson LM, Backdahl L, Wing K. Ncf1 polymorphism reveals oxidative regulation of autoimmune chronic inflammation. Immunol Rev. 2016;269(1):228–47.

    Article  CAS  PubMed  Google Scholar 

  53. Bardos T, Kamath RV, Mikecz K, Glant TT. Anti-inflammatory and chondroprotective effect of TSG-6 (tumor necrosis factor-alpha-stimulated gene-6) in murine models of experimental arthritis. Am J Pathol. 2001;159(5):1711–21.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  54. McDaneld TG, Spurlock DM. Ankyrin repeat and suppressor of cytokine signaling (SOCS) box-containing protein (ASB) 15 alters differentiation of mouse C2C12 myoblasts and phosphorylation of mitogen-activated protein kinase and Akt. J Anim Sci. 2008;86(11):2897–902.

    Article  CAS  PubMed  Google Scholar 

  55. Roger M, Biaso F, Castelle CJ, Bauzan M, Chaspoul F, Lojou E, Sciara G, Caffarri S, Giudici-Orticoni MT, Ilbert M. Spectroscopic characterization of a green copper site in a single-domain cupredoxin. PLoS One. 2014;9(6), e98941.

    Article  PubMed  PubMed Central  Google Scholar 

  56. Xia L, Huang W, Tian D, Chen Z, Zhang L, Li Y, Hu H, Liu J, Chen Z, Tang G, et al. ACP5, a direct transcriptional target of FoxM1, promotes tumor metastasis and indicates poor prognosis in hepatocellular carcinoma. Oncogene. 2014;33(11):1395–406.

    Article  CAS  PubMed  Google Scholar 

  57. Lindblom RP, Aeinehband S, Parsa R, Strom M, Al Nimer F, Zhang XM, Dominguez CA, Flytzani S, Diez M, Piehl F. Genetic variability in the rat Aplec C-type lectin gene cluster regulates lymphocyte trafficking and motor neuron survival after traumatic nerve root injury. J Neuroinflammation. 2013;10:60.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  58. Lee YB, Yune TY, Baik SY, Shin YH, Du S, Rhim H, Lee EB, Kim YC, Shin ML, Markelonis GJ, et al. Role of tumor necrosis factor-alpha in neuronal and glial apoptosis after spinal cord injury. Exp Neurol. 2000;166(1):190–5.

    Article  CAS  PubMed  Google Scholar 

  59. Esposito E, Cuzzocrea S. Anti-TNF therapy in the injured spinal cord. Trends Pharmacol Sci. 2011;32(2):107–15.

    Article  CAS  PubMed  Google Scholar 

  60. Impellizzeri D, Ahmad A, Di Paola R, Campolo M, Navarra M, Esposito E, Cuzzocrea S. Role of Toll like receptor 4 signaling pathway in the secondary damage induced by experimental spinal cord injury. Immunobiology. 2015;220(9):1039–49.

    Article  CAS  PubMed  Google Scholar 

  61. Brambilla R, Bracchi-Ricard V, Hu WH, Frydel B, Bramwell A, Karmally S, Green EJ, Bethea JR. Inhibition of astroglial nuclear factor kappaB reduces inflammation and improves functional recovery after spinal cord injury. J Exp Med. 2005;202(1):145–56.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  62. Kigerl KA, de Rivero Vaccari JP, Dietrich WD, Popovich PG, Keane RW. Pattern recognition receptors and central nervous system repair. Exp Neurol. 2014;258:5–16.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  63. Jacobi A, Loy K, Schmalz AM, Hellsten M, Umemori H, Kerschensteiner M, Bareyre FM. FGF22 signaling regulates synapse formation during post-injury remodeling of the spinal cord. EMBO J. 2015;34(9):1231–43.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  64. Pallier PN, Poddighe L, Zbarsky V, Kostusiak M, Choudhury R, Hart T, Burguillos MA, Musbahi O, Groenendijk M, Sijben JW, et al. A nutrient combination designed to enhance synapse formation and function improves outcome in experimental spinal cord injury. Neurobiol Dis. 2015;82:504–15.

    Article  CAS  PubMed  Google Scholar 

  65. Korostelev A, Trakhanov S, Laurberg M, Noller HF. Crystal structure of a 70S ribosome-tRNA complex reveals functional interactions and rearrangements. Cell. 2006;126(6):1065–77.

    Article  CAS  PubMed  Google Scholar 

  66. Hetman M, Pietrzak M. Emerging roles of the neuronal nucleolus. Trends Neurosci. 2012;35(5):305–14.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

Download references

Acknowledgements

We thank SX Wang and XY Lv for helping us to train rats and perform BBB locomotor rating scale.

Funding

This study was supported by grants from the National Natural Science Foundation of China (Nos. 81271363, 81571194) and Key Program of Anhui Province for Outstanding Talents in University (gxbjZD2016071).

Availability of data and materials

The sequence data of this study have been deposited into Sequence Read Archive (http://www.ncbi.nlm.nih.gov/sra; accession number PRJNA318311).

Authors’ contributions

HZL and JGH participated in the design and coordination of the study, and drafted the manuscript. NZ, XMX and YJC, performed animal experimental procedures and RNA sample collection. LLS, RW, LS and JSZ conducted data analysis and interpretation. All authors read and approved the final manuscript.

Competing interests

The authors declare that they have no competing interests.

Consent for publication

Not applicable.

Ethics approval

All surgical procedures and post-operative animal care were in accordance with the Regulations for the Administration of Affairs Concerning Experimental Animals (Ministry of Science and Technology, China; revised in June 2004) and the Guidelines and Policies for Rodent Survival Surgery provided by the Animal Care and Use Committees of Bengbu Medical College.

Author information

Author notes

    Authors and Affiliations

    1. Clinical Laboratory, the First Affiliated Hospital of Bengbu Medical College, Anhui, 233004, People’s Republic of China

      Ling-Ling Shi, Xiu-Mei Xie, Yue-Juan Chen, Lin Shen, Jian-Guo Hu & He-Zuo Lü

    2. Anhui Key Laboratory of Tissue Transplantation, the First Affiliated Hospital of Bengbu Medical College, 287 Chang Huai Road, Anhui, 233004, People’s Republic of China

      Nan Zhang, Xiu-Mei Xie, Yue-Juan Chen, Rui Wang, Lin Shen, Jian-Sheng Zhou, Jian-Guo Hu & He-Zuo Lü

    3. Department of Immunology, Bengbu Medical College, Anhui, 233030, People’s Republic of China

      Ling-Ling Shi & He-Zuo Lü

    Authors
    1. Ling-Ling Shi
      View author publications

      You can also search for this author in PubMed Google Scholar

    2. Nan Zhang
      View author publications

      You can also search for this author in PubMed Google Scholar

    3. Xiu-Mei Xie
      View author publications

      You can also search for this author in PubMed Google Scholar

    4. Yue-Juan Chen
      View author publications

      You can also search for this author in PubMed Google Scholar

    5. Rui Wang
      View author publications

      You can also search for this author in PubMed Google Scholar

    6. Lin Shen
      View author publications

      You can also search for this author in PubMed Google Scholar

    7. Jian-Sheng Zhou
      View author publications

      You can also search for this author in PubMed Google Scholar

    8. Jian-Guo Hu
      View author publications

      You can also search for this author in PubMed Google Scholar

    9. He-Zuo Lü
      View author publications

      You can also search for this author in PubMed Google Scholar

    Corresponding authors

    Correspondence to Jian-Guo Hu or He-Zuo Lü.

    Additional files

    Additional file 1: Table S1.

    Differentially expressed genes in the spinal cords at different injured stages. The downregulated genes are shown in green, the upregulated are shown in Red. (XLSX 634 kb)

    Additional file 2: Table S2.

    The list of most relevant genes that are over-expressed or under-expressed in the different injured stages. (XLSX 44 kb)

    Additional file 3: Table S3.

    The list of gene IDs in 8 subclusters for K-means clustering of DEGs. (XLSX 360 kb)

    Additional file 4: Table S4.

    The significant GO terms of differentially expressed genes. (XLSX 109 kb)

    Additional file 5: Figure S1.

    GO enrichment analysis of up-regulated genes in acute SCI (1 dpi). The 30 most enriched GO terms are shown. The asterisks (*) represent the significantly enriched (P ≤ 0.05) biological process, cellular component, and molecular function categories. (TIF 1475 kb)

    Additional file 6: Table S5.

    GO enrichment result of the up-regulated genes. (XLSX 64 kb)

    Additional file 7: Figure S2.

    GO enrichment analysis of down-regulated genes in acute SCI (1 dpi). The 30 most enriched GO terms are shown. No significantly enriched biological process, cellular component, and molecular function categories are detected. (TIF 1595 kb)

    Additional file 8: Table S6.

    GO enrichment result of the down-regulated genes. (XLSX 58 kb)

    Additional file 9: Figure S3.

    GO enrichment analysis of up-regulated genes in subacute SCI (6 dpi). The 30 most enriched GO terms are shown. The asterisks (*) represent the significantly enriched (P ≤ 0.05) biological process, cellular component, and molecular function categories. (TIF 626 kb)

    Additional file 10: Figure S4.

    GO enrichment analysis of down-regulated genes in subacute SCI (6 dpi). The 30 most enriched GO terms are shown. The asterisks (*) represent the significantly enriched (P ≤ 0.05) biological process and molecular function categories. No significantly enriched cellular component is detected. (TIF 646 kb)

    Additional file 11: Figure S5.

    GO enrichment analysis of up-regulated genes in chronic SCI (28 dpi). The 30 most enriched GO terms are shown. The asterisks (*) represent the significantly enriched (P ≤ 0.05) biological process, cellular component, and molecular function categories. (TIF 1644 kb)

    Additional file 12: Figure S6.

    GO enrichment analysis of down-regulated genes in chronic SCI (28 dpi). The 30 most enriched GO terms are shown. The asterisks (*) represent the significantly enriched (P ≤ 0.05) biological process, cellular component, and molecular function categories. (TIF 1699 kb)

    Additional file 13: Table S7.

    KEGG enrichment analysis of up DEGs in acute SCI (1 dpi). (XLSX 192 kb)

    Additional file 14: Table S8.

    KEGG enrichment analysis of up DEGs in subacute SCI (6 dpi). (XLSX 384 kb)

    Additional file 15: Table S9.

    KEGG enrichment analysis of up DEGs chronic SCI (28 dpi). (XLSX 267 kb)

    Additional file 16: Figure S7.

    KEGG enrichment analysis of DEGs in acute SCI (1 dpi). The 20 most enriched KEGG pathways in 1 dpi vs sham. “Rich factor” means that the ratio of the DEGs number and the number of genes has been annotated in this pathway. The greater of the Rich factor, the greater of the degree of enrichment. (TIF 1729 kb)

    Additional file 17: Figure S8.

    KEGG enrichment analysis of DEGs in subacute SCI (6 dpi). The 20 most enriched KEGG pathways in 6 dpi vs sham. “Rich factor” means that the ratio of the DEGs number and the number of genes has been annotated in this pathway. The greater of the Rich factor, the greater of the degree of enrichment. (TIF 1827 kb)

    Additional file 18: Figure S9.

    KEGG enrichment analysis of DEGs in chronic SCI (28 dpi). The 20 most enriched KEGG pathways in 28 dpi vs sham. “Rich factor” means that the ratio of the DEGs number and the number of genes has been annotated in this pathway. The greater of the Rich factor, the greater of the degree of enrichment. (TIF 1748 kb)

    Rights and permissions

    Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated.

    Reprints and Permissions

    About this article

    Verify currency and authenticity via CrossMark

    Cite this article

    Shi, LL., Zhang, N., Xie, XM. et al. Transcriptome profile of rat genes in injured spinal cord at different stages by RNA-sequencing. BMC Genomics 18, 173 (2017). https://doi.org/10.1186/s12864-017-3532-x

    Download citation

    • Received:

    • Accepted:

    • Published:

    • DOI: https://doi.org/10.1186/s12864-017-3532-x

    Keywords

    • Spinal cord injury
    • RNA-Seq
    • GO enrichment
    • Pathway analysis
    • Sprague-Dawley rats (RRID:RGD_70508)

    BMC Genomics

    ISSN: 1471-2164

    Contact us