Genome-wide gene expression profiling of stress response in a spinal cord clip compression injury model

Feature analysis of affymetrix GeneChip Rat genome 230 2.0 array

Analysis and filtering on the resulting file of 31,042 ProbeSets revealed that 10,791 ProbeSet IDs had no annotations, i.e. no Entrez IDs or official gene symbols, which were flagged out. This reduced the number of workable ProbeSet IDs to 20,251. In addition, there were duplicate or multiple ProbeSet IDs which represented a single gene. Conversely, there were ProbeSet IDs with multiple annotations (EntrezID/Gene Symbol) due to sequence identity across more than one gene segment in the genome. This issue could not be easily resolved as the level of uniqueness of the oligonucleotide sequence is not high enough to allow annotation to one gene exclusively. This feature requires manual curation of the data based on Affymetrix instructions to use the latest annotation, which is also the most relevant. Taking the above two features into consideration results in 14,324 gene symbols on the GeneChip RG230 2.0 array. The resulting data file still contains the ProbeSet IDs that have “LOC” or “RGD” identifiers instead of actual gene symbols. These identifiers are applied to genes that are less well characterized and usually belong to similar or orthologous proteins in other species. They may also belong to non-coding regions of the genome. Software platforms developed for GO enrichment and pathway analysis rarely map the LOC and RGD identifiers. In the original Affymetrix GeneChip 230 2.0 array annotation file, the number of LOC and RGD annotated ProbeSets sets are 1163 and 1135, respectively. The same issue of duplicate/multiple entries also applies to these ProbeSets, hence the numbers of “LOC” and “RGD” identifiers in the final output file with 14,324 entries were less: 939 for the LOC and 829 for the RGD identifiers. This means that the total number of annotated ProbeSets in the Affymetrix GeneChip Rat Genome 230 2.0 array annotation file that were mapped to known gene candidates was 12,557, which is equivalent to 62.4% and 71.2% of the total number of genes annotated and listed in the Rat Genome Database (RGD) and European Bioinformatics Institute (EBI) association files, respectively.

Analysis of ProbeSet data

We used a divisive hierarchical clustering algorithm (DIANA) to identify the strongest trends within the dataset, in all pair-wise comparisons (see Methods). The results were visualized using the Heatplus package of BioConductor (Figure 1A). Interestingly, three main clusters were observed, corresponding to the sham controls (green cluster), day 1 animals (blue cluster), and day 3 animals (yellow cluster). The later time-points are further subdivided into a cluster composed primarily of day 7 animals (red cluster) and another composed primarily of day 14 and day 56 animals (brown cluster). Principal component analysis of the ProbeSet data was performed to assess variability between individual animals in each group and also at different time points and resulted in clustering of the transcripts belonging to each group of sham or injured animals shown in Figure 1B. There are inter-individual differences but eclipses show that there are no outliers in our experiment. Additionally, the eclipses of day 7, day 14 and day 56 cross each other, which indicate some level of commonality between these time points, as was evidenced and shown in the tree view of the heat map.

Data normalization and expression/signal value determination resulted in a list of all 31,099 ProbeSets, their fold change values relative to sham (in Log₂ scale), and associated ANOVA t test p-values across the time points. Volcano plots of the corresponding fold change values against transformed (−log₁₀) p-values for every time point are displayed in Figure 1C. As shown, all volcano plots display a normal distribution of ProbeSets with fold change values from −8.7 to 11.2 for down- and up-regulated genes, respectively. The shape of the volcano plot changes as time post-injury goes by. Thus, day 3 ProbeSet data plots are not as populated, especially on the down-regulated area and are less similar to other data points. The day 1 plot, on the other hand, looks more similar to the day 7 volcano plot. The more chronic data points of day 14 and day 56 look more similar to each other than to earlier data points.

Examination of the number of ProbeSets with marginal ANOVA t test p-values gave an estimate as to the reliability of data obtained. Thus, we analyzed our data for the number of ProbeSets with ANOVA t test p-values higher than 0.05 at different fold change values (Figure 1D). We found that the majority of changes in gene expression with significance levels of p > 0.05 generally belong to ProbeSets with lower fold change values. For example, the number of ProbeSets with ANOVA t test p > 0.05 did not exceed 6% of the total number of ProbeSets, irrespective of the fold change values. At a more stringent significance level of p ≤ 0.001, however, it would be necessary to filter out the ProbeSets with expression values less than 2 fold changes in order to keep the number of filtered ProbeSets around 10% or less across the time points (data not shown). Thus, filtering the data with higher fold change values automatically targets for transcripts with smaller t test p-values. Based on the results presented in Figure 1D, we performed the functional analysis on the ProbeSet data with fold change values of  ≥ 1.5 and p ≤ 0.05.

Analysis of gene set data

We next examined the nature of deregulated transcripts at different time-points relative to each other. Figure 2 (A-D) depicts the Venn diagrams with overlapping regions demonstrating the number of common genes showing changes at various time points as well as time-point specific genes. In terms of gene contents, the day 1 pattern of gene expression is more similar to day 7 as is evidenced by 760 common genes between day 1 and day 7, compared to 317 between day 1 and day 3, 186 between day 1 and day 14 and 113 between day 1 and day 56. On the other hand, each time-point has unique genes, whose expressions do not appear to change at other time points. This observation supports the notion that, although some processes that are invoked early after SCI may stay active throughout the acute or chronic phase, there are unique features to the early response genes that are dramatically different from the response in the following days or weeks post-injury. Additionally, deregulated transcripts on day 14 and day 56 were found to be very similar to each other with approximately 82% of the genes showing changed expression being identical at these two time-points. This result was also predicted from the heat map (Figure 1A). This indicates that the biological processes in response during the chronic phase of SCI remain constant.

Time-series expression profile clustering by STEM

To simplify the graphical presentation of the data, fold changes in expression values for all genes associated with only the statistically significant profiles were averaged and plotted against the post-injury observation time points (Figure 3A-H). Two broad profile classes become apparent from these data: "up-down/down-up genes” (Figures 3A-E) and “fluctuating genes” (Figures 3F-H). The “up-down” category comprises the cluster of profiles that share a pattern of up-regulation early after injury over the course of 24 hours (profiles 45, 46 and 48). In these profiles, the early response is subsequently followed by down-regulation of the genes clustered in these profiles. This occurs gradually, approaching a steady level at normal or higher than normal values. Although genes in profiles 45 and 48 were allocated to two separate expression profiles by STEM, visual examination indicates striking similarities between the two profiles. Patterns of late phase expression in these profiles look similar to each other as the late expression values of almost all gene transcripts in both profiles approach values comparable to those of the sham animals. For many genes in profile 46 the pattern is different as the transcript levels of about 50% of genes in this profile remain up-regulated throughout the course of the study and also at the end of 8 weeks, hence their functions seem essential not only in the acute and subacute phase but also during the chronic stage of the injury. Profiles 1 and 0 are quite similar to each other as they display down-regulations of many genes on day 1, which stay at lower than normal levels even 8 weeks after injury.

Class II profiles represent fluctuating profiles (44, 41 and 6), with a surprising but more complex pattern of gene expression, most notably during the 24–72 hours post-injury (Figure 3G-H). This results in a bi-phasic expression pattern, which falls into two main clusters. The first cluster comprises profiles 44 and 41 (Figure 3F-G) and is characterized by an initial up-regulation of gene transcript levels early on day 1, followed by a sharp decrease in gene expression on day 3. More than 53% and 83% of the genes in profiles 44 and 41, respectfully, displayed at least a 1.0 (Log₂ scale) fold change reduction in transcript levels on day 3 compared to day 1. For profiles 44 and 41, this bi-phasic pattern of gene expression is further followed by escalation of gene expression, which peaks at day 7 and stabilizes on day 14 onward. The second cluster only includes profile 6, which is essentially the mirror of profile 44 and comprise down-regulated genes (Figure 3H). It is characterized by an early and substantial down-regulation on day 1. Next, a period of recovery to normal transcript levels is observed that peaks on day 3 post-injury and then switches direction again and remains low through today 56. Finally, detailed information in Figure 3 indicates that the majority of transcripts belong to profile 44 and 6 with up-regulated transcripts clustering in the former and down-regulated transcripts in the latter.

In summary, the following conclusions can be drawn from the cluster analysis of transcripts both at the ProbeSet and gene level following clip-compression injury of the spinal cord in rats:

GO enrichment analysis of deregulated genes

Fold change and p-value criteria affect the number of enriched terms

GO enrichment analysis can result in numerous enriched GO terms with overlapping or redundant terms making the reduction and prioritization task difficult. Depending on the fold change value criteria, which determine the number of deregulated transcripts, the number of enriched GO terms can vary. In order to rationalize an approach where meaningful numbers of GO terms are achieved, we first examined how the parameter of fold-change in expression can affect the number of enriched terms. We chose to perform this preliminary analysis in a time-series fashion, meaning that deregulated transcripts with significant fold change values (ANOVA t test p ≤ 0.05) across all time-points were considered and only one missing value was permitted. Thus, deregulated transcripts at 1–4 fold change values in at least one time point were separately subjected to GO Biological Process (BP) enrichment analysis (minimum GO level of 3, minimum number of genes of 5) and the number of GO terms were plotted as a function of fold change in expression (Figure 4A). As shown, this time-series analysis at 1.0, 1.5 and 2.0 fold change values resulted in very high numbers of enriched terms, e.g. 698, 649 and 720 GO terms (with the adjusted p-values ≤ 0.001), respectively. Setting the fold change criteria to higher values did not limit the number of GO terms as it only gradually declined (Figure 4A). For example, GO enrichment analysis on the list of transcripts with 2.5, 3.0, 3.5 and 4.0 fold change values resulted in 625, 590, 487 and 276 terms (p ≤ 0.001), respectively. It is interesting to note that significant reduction in the number of GO terms is not achieved until much higher fold change values are considered, e.g. 4 fold in Log₂ scale (64 fold in normal scale), which also significantly reduces the number of deregulated transcripts included in the analysis (data not shown). However, examining the number of terms at lower p-values remarkably reduced the number of GO terms, although similar trends across different fold change values were observed (Figure 4A). For example, changing the p-value parameter from 10^-3 to 10^-4 reduced the number of enriched GO terms to about half or even less, at all fold change values (Figure 4B). Therefore, it seems that the logical approach for enrichment analysis in regard to the fold change and p-value criteria is to include all transcripts with even lower fold change values to include higher number of genes and especially to avoid excluding potential key regulatory genes in the analysis, which may not have displayed dramatic changes in expression. Finally, although the actual size p-value cutoff seems to be a much more important parameter in limiting the number of enriched GO terms to a workable value, it may not help with prioritization of the enriched terms and their specificity.

Time-series GO enrichment

Based on the results obtained from analysis of the effects of fold change, p-value cut off and GO level criteria, the pool of deregulated transcripts throughout all time-points was analyzed by setting the GO level at different values with the intention of obtaining more specific categories. The enrichment analysis at p-value cutoff of 10^-4 on transcripts with a minimum of 1.5 fold changes in expression (ANOVA t test p ≤ 0.05) resulted in a significant reduction of the number of enriched GO terms to 329 at GO level 3 and higher. Within this collection of enriched GO terms, there are 267 terms whose GO levels are 5 and higher. The 329 and 267 terms along with their p-values were further summarized independently by the REViGO reduction analysis tool that condenses the GO description by removing redundant terms [40]. The results of these further reductions are visualized in Figures 5. Only categories with lower dispensability and frequency, i.e. more uniqueness are shown. As shown, more general terms such as adaptive and innate immune response, immune effector process, cell-cell signaling, cell communication, cell adhesion, cell activation and phagocytosis are significant at level 3 and higher (Figure 5A). Other, more specific terms are also visualized such as regulation of apoptotic process, I-kappaB kinase/NF-kappaB cascade, glial cell migration, and synaptic transmission and synapse organization. However, to achieve more specific terms, one needs to look at the terms with higher GO levels (Figure 5B). The majority of terms visualized in Figure 5B are specific terms such as mitotic cell cycle G2/M transition checkpoint, apoptotic cell clearance, hydrogen peroxide biosynthesis, signal transduction by p53 class mediator and regulation of Toll-like receptor signaling pathway.

Most enriched GO terms visualized in Figure 5 represent deregulated transcripts with the different expression profiles shown in Additional file 1: Figure S1. As expression clustering is performed prior to GO enrichment, the time-series GO enrichment may produce enriched terms that are also significantly represented by the expression profiles shown in Additional file 1: Figure S1. Indeed, we found that the majority of enriched GO terms for up-regulated transcripts represent the expression profiles 44 and 46 and in some cases, profile 48.

Time-point GO enrichment

We next analyzed the temporal pattern of each GO term in a time-point fashion in order to examine the order of events after SCI. To accomplish this, we made multiple comparisons of the enriched GO terms obtained for deregulated transcripts (minimum fold change value ≥ 1.5; ANOVA t test p-value ≤ 0.05) at individual time-points. Figure 6 depicts the Venn diagram of this analysis on the multiple lists of enriched GO terms (p-value cut off of 0.05) obtained for each individual time-point. A less stringent condition was selected to allow for all possible similarities and differences to be observed. As shown, this analysis resulted in 736 common GO terms between all time-points, of which 284 had a p-value ≤ 0.00001 throughout the course of the study (data not shown). Additionally, some biological processes were shown to be significantly up- or down-regulated at only a certain time window as their respective GO terms were uniquely specified to one time point only. For example, 278, 359 and 170 terms were uniquely specified to day 1, day3 and day 7 post-injury, respectively. There are fewer unique terms detected at the chronic stage of the injury i.e. 69 and 67 for day 14 and day 56 post-injury, respectively. The fact that there are significantly higher numbers of common GO terms (i.e. 736) relative to the number of unique terms at each time-point indicates the complexity of the many common processes involved following moderate to severe SCI and that the significance of contribution of these processes is diminished within our time window. Having determined the GO levels for all categories in the previous steps, we then determined the most specific terms and positioned them in their GO tree hierarchy. A summary of significantly enriched BP terms that were found to be commonly up-regulated across all time-points is presented in Table 3. Some general terms such as the “response to external stimulus”, “response to mechanical stimulus” and “inflammatory response” possess the highest significance but belong to lower levels of GO hierarchy. Other categories, however, present more specific functions and are positioned at higher GO levels. As shown, significant up-regulation of GO terms corresponding to “response to extracellular and mechanical stimuli”, “inflammatory response” as well as “response to lipid” and “response to lipopolysaccharide” was observed across all time points, day 1-day 56 post-injury. The “angiogenesis” term is also consistently enriched in the day 1-day 56 time points. Genes involved in “blood coagulation” were up-regulated to high levels at day 1, and at day 7 through today 14. The “complement activation” process is only enriched at later time points, i.e. day 14-day 56 post-injury. Although a significant “response to hypoxia” was observed from day 1-day 3, the peak of response to “reactive oxygen species” and “hydrogen peroxide” was observed exclusively on day 1. The “response to glucocorticoid stimulus” was observed on day 1 and day 14 post-injury. “Production of IL-6” process peaked on day 1 post-injury but continued to be enriched on day 7 and at later time points, while the peaks for “tumor necrosis factor production” and “IL-8 production” were on day 7. Previous studies have shown a marked increase in TNF-alpha production immediately after injury. Since this study did not include time-points earlier than 24 hours, this result may imply the second wave of TNF-alpha production following the acute phase of the injury. Interestingly, the responses to cytokines such as tumor necrosis factor and IL-1 were observed only on day 3. “Neutrophil chemotaxis” process initiated on day 1 but peak activity was observed on day 14. From Figure 7H, it is evident that the level of “phagocytosis”-related transcripts are up-regulated early after injury but peak at one week post-injury and stay up-regulated. Likewise, the expression of genes that belong to the “Toll-like receptor signaling pathway” term alters with the same profile (Figure 7K).

A significant finding is the occurrence of the “apoptotic process” on day 1-day 7 post-injury. This process is accompanied by events whose peak of response also corresponds exclusively to day 7, such as “interleukin-6 production”, “tumor necrosis factor production”, “macrophage activation involved in immune response”, “phagocytosis” and “engulfment and “apoptotic cell clearance” (Figure 7E and 7H).

On both days 1 and 7, a significant up-regulation in induced-innate immunity related GO terms such as “pattern recognition-mediated signaling’ , “Toll-like receptor signaling” and “integrin-mediated signaling pathways” was detected. “Leukocyte cell-cell adhesion” was observed from day 1 to day 7. While genes involved in activation of “innate immune response”, “B and T cell activation”, “cytokine biosynthetic process”, and “phagocytosis” were up-regulated at day 1 and from day 7 onwards to day 56; ‘T cell differentiation” and “B cell mediated immunity” up-regulation is only observed during the chronic phase of injury, i.e. day 14-day 56 (Figure 7O and 7P). Thus, it is not surprising that the “B and T cell proliferation” and the “B cell receptor signalling pathway” peaks of response were on day 7-day 14. Day 14 also marks a peak response to “ossification” (data not shown). Importantly, the peak response to interferon-gamma and the immunoglobulin-mediated immune response is observed on day 56. These two mark the late response biological processes induced after injury to spinal cord.

Our analysis also showed that “programmed cell death” and its related child terms “apoptotic process” and “positive and negative regulation of apoptotic process” are commonly enriched only during day 1- day 7 post-injury. Apoptotic processes significantly increase early after injury on day 1 post-injury and reach a peak between day 3 and day 7 post-injury, after which the contribution of apoptotic processes is diminished but stays significantly enriched (p ≤ 0.001-0.0001). Both positive and negative regulations of apoptosis are significantly enriched, which indicates the fact that the injured cells struggle for survival. However, activation of apoptosis seems to be more predominant than its suppression, as the positive regulation of apoptosis becomes activated earlier than negative regulation and its peak of activity is on day 1 post-injury, although it stays continuously up-regulated up to 1 week post-injury. In contrast, the only significant activity of negative regulation of apoptosis (p ≤ 0.00001) is on day 3 (Figure 7I-J).

We can summarize the biological processes listed in Table 3 into three main categories: stress response including processes such as blood coagulation, complement activation, response to hypoxia and reactive oxygen species, angiogenesis and inflammation. The second category consists of induced innate immune response processes such as activation of macrophages and microglia by Toll-like receptor signaling, cytokine production and secretion, chemokine production and neutrophil chemotaxis, IL6 and tumor-necrosis factor production and their responses. A significant set of induced innate immune-related biological processes involve “Phagocytosis” and “Toll-like receptor signaling pathway”. The third category is mainly the components of adaptive immune response processes such as T cell activation, migration and proliferation, B cell activation and immunoglobulin-mediated immunity. Both humoral and cell-mediated elements of “adaptive immune response” processes seem to actively participate in the pathology of SCI. The pattern of change in mRNA levels for many genes in the above GO biological processes follow the expression profiles observed in profiles 44, 46, 48 and 45, which have been discussed earlier (Figure 3A-H). The continuous up-regulation of the “immunoglobulin mediated immune response” and its parent term “B cell mediated immunity” is striking and may imply that these processes should be categorized as chronic phase responses to SCI as their peak of activity appears after 2–8 weeks, although initiated at early timepoints post-injury.

Our GO analysis also resulted in enrichment of many regulatory processes, the majority of which are positive regulations of the enriched GO biological processes listed in Table 3, as illustrated for apoptotic processes shown in Figures 7I-J.

The injury model

Since its introduction in 1978 as the first SCI model in rats [20], the clip-compression model has become a standard injury model in animals as it mimics the human primary mechanism of injury to the spinal cord as well as the histo-pathological and behavioural outcomes of human SCI. Our lab has previously characterized this mode of SCI [42–44]. The clip compression of the spinal cord results in central cavitation and axonal loss in the white matter of spinal cord [45]. Rats that receive the clip-compression injury have a very similar pathological progression to humans with SCI including the formation of a cystic cavity surrounded by a glial scar [44]. In addition, animals injured by clip compression will have the same functional recovery profile as is observed in humans [20, 23, 28].

Previous studies have shown that the response of the spinal cord tissue to injury consists of a complex series of cellular responses and events. These cellular events are reflected in a more complex change in temporal and spatial pattern of molecular events at the mRNA level, which, in turn, depends on the type and severity of the primary injury and the following cascade of secondary events [1]. Earlier reports on high throughput gene expression analysis after SCI in animals have been almost exclusively performed in contusion-based models of injury using weight drop method [36–38, 46]. As no such study on the clip-compression injury model has been reported, we aimed to investigate the rat transcriptome dynamics after a moderate to severe injury using hemorrhagic SCI by clip model, similar to most human SCIs. Additionally, the primary injury in the clip model consists of both impact and persistent compression. Therefore, we hypothesized that both similarities and differences between the two models of injury would be evident by examining how the changes in transcriptome occur. Moreover, unlike the majority of earlier studies that chiefly examined the acute and subacute events, we extended the time-frame of our study to 8 weeks post-injury to allow examination of the acute, subacute and chronic phases of the injury. The chosen time-points were based on previous behavioural and immunohistochemical analyses, which showed that following SCI by clip-compression, the first 24 hours post injury would represent a very acute stage and possible involvement of most immediate early stress genes. Days 3 and 7 represent a time during which the peak of delayed apoptotic cell death for the neural cells occurs. Days 10–14 are considered the subacute stage, as the inflammation appears to subside. Finally day 56 is considered the chronic stage as it is the time when the BBB motor recovery test for the spontaneous recovery/improvement in the rat animal model reaches a plateau.

GO enrichment analysis as a tool for biological process inference

Functional analysis of microarray data is a challenging task as the result of initial analysis is only the fold change values representing deregulations in the expression of thousands of transcripts. There are different approaches to analyzing the results of a microarray experiment in order to make efficient biological inferences. Various platforms share a common feature in that they perform an overrepresentation analysis on the list of deregulated genes and statistically analyze if the pool of up- and/or down-regulated transcripts is significantly enriched compared to the list of genes previously annotated to be part of a defined Biological Process, Molecular Function or Cellular Component, as is the case with GO enrichment, or to a certain metabolic or signaling pathway as is observed in pathway analysis platforms. Various pathway analyses are currently in practice for microarray data analysis and there are different approaches to accomplish this. KEGG pathway [47–49], Wikipathways [50–52], and Ingenuity (http://www.ingenuity.com) are amongst the currently available platforms for pathway analysis. A recent analysis showed that among the above three pathway databases, (KEGG, Ingenuity and Wikipathways) there is a low level of consistency, comprehensiveness and compatibility [53] and the level of consistency varies significantly when different pathways are compared. Due to these limitations, and because GO is considered to represent a relatively current, comprehensive, and, more importantly, a controlled vocabulary for gene function [54], we analyzed our microarray data using GO enrichment analysis. However, we are also aware of the limitations of GO enrichment analysis [55]. For example, prior to GO enrichment analysis in this study, we determined the number of annotated genes in the list of deregulated transcripts and found that only 55% and 75% of the 14,327 genes on the Rat GeneChip 230 2.0 are annotated in the EBI and RGD association files, respectively (data not shown). The above percentages of annotated genes in Rat genome are similar to the number of annotated genes in all other organisms whose genome has been sequenced and only a subset of known genes are annotated for each of the three domains of GO tree, i.e. BP, MF and CP components [56].

An advantage in using a controlled vocabulary of gene function such as GO on the SCI microarray data comes from the challenging nature of such analysis due to the inherent complexity of the spinal cord tissue and also the type and level of injury itself. Spinal cord tissue is composed of an array of highly specialized neurons, astrocytes, oligodendrocytes, microglia, and pericytes. Another specialized and complex structure within the cord tissue whose permeability is highly compromised [57] upon injury is the blood spinal cord barrier (BSCB), which is composed of neurovascular unit (NVU), that maintains the integrity of BSCB and is again comprised of endothelial cells, neurons, astrocytes, and pericytes [58]. Additionally, SCI is generally categorized as a severe injury that leads to loss of normal physiological functions. Thus, the development of a complex series of secondary damage [1] to the spinal cord after the primary injury is due both to the vast array of cell types affected as well as the injury severity that sets many processes in motion. Such an injury model demands a nonbiased and yet comprehensive coverage of annotations such as GO for clustering of deregulated genes into relevant processes and events. The reliability of this approach is shown by its successful conjecturing of previously known biological processes as well as their dynamic of contribution to the pathology of spinal cord injury as explained below.

Blood coagulation and blood protein signaling

The supply of blood and nutrients is crucial for normal functioning of neural cells. It is well-documented that an early and progressive development of hemorrhage is a common feature of all experimental models of SCI and this includes the clip-compression model [59, 60]. Shearing of the blood vessels and disruption of the vascular architecture within the lesion epicenter by mechanical force leads to hemorrhage, a progressive process which extends to the rostral but more towards the caudal regions of the grey matter [24, 61–63]. As post-traumatic ischemia develops [1, 59], further vasospasm [64] and loss of autoregulation of blood flow [65, 66] exacerbate the condition. Therefore, the earliest event following compression injury to the spinal cord is a profound damage to the local vasculature (capillaries and venules), hemorrhage (especially in the grey matter) and disruption of cord microcirculation by mechanical, thrombotic or vasospasm mechanisms. Consequently, the normal blood flow to the spinal cord is significantly reduced, which leads to a marked ischemia in the gray and white matter [60].

The results of our microarray data analysis clearly confirm the outcome of the primary impact and persistent compression injury to the spinal cord, which is disruption of the vasculature and hemorrhage as the major and initial result of the primary injury. Our data indicate that representative genes in the blood coagulation cascade are up-regulated (Figure 7A). For example, the transcript levels of the integral membrane protein tissue factor (coagulation factor III, F3), coagulation factors VIII (F13A1, F8), platelet factor (PF4) and V (F5) are up-regulated, the latter being elevated only on day 1 post-injury (data not shown). Permanent binding of tissue factor F3 to membrane surface is thought to be crucial for the speed of enzymatic reactions in coagulation processes [67]. Additionally, we found that platelet factor (PF4) mRNA levels were increased upon injury. PF4 (CXCL4) is a chemokine released from activated platelets to bind heparin and inhibit its anticoagulant activity. ANO6 is a transmembrane protein that may have a calcium-activated chloride channel activity but it is thought to be essential for calcium-dependent exposure of phosphatidylserine on the surface of activated platelets. Importantly, ANO6 (anoctamin 6 or TMEM16F) transcript level is also elevated early after injury and is continues to be up-regulated up to 8 weeks post-injury. Higher than normal transcript levels of ANO6 during both acute and chronic phases of SCI may explain why the coagulation process is up-regulated even at 8 weeks post-injury. Regulatory proteins such as protein C, a serine protease that is activated in the blood coagulation cascade, along with its receptor (PROCR) are up-regulated as well. Activated protein C has potent anticoagulant activity due to its ability to inactivate factor Va and VIIIa (Yesilirmak et al., 2008) and seems to alleviate the secondary SCI by reducing the ischemia/reperfusion effect by inhibiting neutrophil activation (Hirose et al., 1999) and or leukocyte activation [68], inducing insulin growth factor-1 and its receptor leading to an increased number of motor neurons [69].

The GO enrichment analysis identified another 30 coagulation-related genes whose transcripts were up-regulated throughout the course of the study. Amongst these were regulatory proteins with anticoagulant properties such as tissue factor pathway inhibitor 2 (TFPI), which is released by endothelial cells and binds factor VIIa complexes, inhibiting them to generate factor Xa. TFPI function regulates the extrinsic coagulation pathway. Additionally, we found that thrombomodulin (THBD) transcripts were elevated upon SCI up to 2 weeks post-injury. THBD binds thrombin and promotes its interaction with protein C. The resulting complex inactivates factors VIIIa and Va. Elevated levels of these regulatory proteins indicate the importance of endogenous signaling mechanisms to limit excessive spreading of clot formation.

A serious side effect of hemorrhage is the infiltration of blood components such as hemoglobin and fibrinogen to the spinal cord tissue which have been shown to be toxic to CNS tissue [70–73]. Infiltration of hemoglobin creates a hostile environment that is rich in reactive oxygen species and other toxic materials, which induces the cellular response to these toxic mediators of cell death and apoptosis. Hemoglobin, released from red cells after trauma, can promote tissue injury through iron-dependent mechanisms such as inhibiting the Na/K ATPase activity and catalyzing substantial peroxidation of CNS lipids [71]. In our study, the majority of Na/K ATPase enzymes such as ATP1A2, ATP1A3, ATP1B1 and ATP1B2 were down-regulated during the acute as well as the subacute phase of the injury (data not shown). Fibrinogen has been shown to trigger an inhibitory signal transduction pathway in neurons by acting as a ligand for beta-3 integrin, which induces the transactivation of EGF receptor (EGFR) in neurons, thereby inhibiting neurite outgrowth [73]. It also triggers astrocyte scar formation through TGF-beta signaling [72]. The microarray data in our study confirms that genes in the TGF-beta signaling cascade are up-regulated. For example, TGFB1, its receptor and SMAD2 transcripts were up-regulated throughout the 8 weeks post-injury study period (data not shown).

Complement activation

Along with the blood coagulation cascade, a concomitant increase in the complement activation system is observed, whose temporal pattern is not the same as blood coagulation but rather develops in a more delayed fashion. The blood coagulation cascade peak of activity is on day 7 post-injury but stays up-regulated until 8 weeks. The complement activation, however, is turned on with a lag time in the first few days with activity increasing at later time points in the experiment (Figure 7C). Whether the late activation of complement system is due to the effect of reperfusion after ischemia needs further investigation. The complement system can be activated by three different but overlapping classical, lectin and alternative routes [74].

Representative genes in the complement activation system were deregulated following clip injury to spinal cord. For example, the transcript level of the main activator of the classical pathway of complement activation (C1S) is down-regulated one day after injury. However, it returns to normal values by day 3 and is further up-regulated by day 7 remaining at higher than normal levels even at day 56 post-injury (data not shown). C1S catalyzes the consecutive conversion of C4 to C4a and C4a to active C4b2a (C3 convertase), whose main function is to cleave parental C3 into C3a and C3b. As shown the mRNA levels of C1qa, C1qb, C1qc, Cfd and Cr1l are increased relative to sham un-injured animals. The transcript level of Factor H (CFH), a negative regulator of the alternative pathway for complement activation, is decreased after injury but fluctuates back to higher than normal levels by day 7 post-injury. The elevated level of CFH in our study is in agreement with previous reports that complement inhibitor proteins such as factor H were expressed at elevated levels on neurons and oligodendrocytes after SCI in rats [75, 76].

Using inhibitor approaches, both classical and lectin pathways of complement activation have been shown to participate in SCI pathology [77–79]. C1q Knockout mice showed improved recovery and thus the classical complement activation via C1q is thought to be detrimental to the injured spinal cord [80]. Our data show that the mRNA level of C1 inhibitor (C1-INH, SERPING1), an inhibitor of the lectin pathway, is also increased in a similar profile as observed in CFH mRNA deregulation. C1-INH inhibits complement activation through binding and inactivating MASP1 and MASP2 [74]. Up-regulation of C1-INH has been shown to be protective and independent of C1q and the classical pathway [81].

Ischemia, response to hypoxia and reactive oxygen species

The decrease in the local blood-flow leads to ischemic-hypoxic damage to the spinal cord tissue. Ischemia generally leads to a decrease in cytoplasmic levels of ATP, cellular swelling through malfunctioning of Na/K ATPases and also the mitochondrial membrane permeability transition [82]. Additionally, hypoxia induces certain transcription factors such as hypoxia inducible factor 1 (Hif-1) heterodimer which is composed of the inducible Hif-1a and the constitutive Hif-1b subunits [83, 84]. The induction of Hif-1a is under the control of NF-kB transcription factor which serves to link hypoxia to innate immune response [85]. This is reflected in an increase in the mRNA level for the genes that function in response to hypoxia. We found that, following clip-compression injury to the spinal cord, the transcript levels of Hif-1a were up-regulated. Negative control of Hif-1 transcriptional activity is under the control of EGLN3, a propyl hydroxylase that, in the presence of oxygen molecule, permits ubiquitination and proteosomal degradation of Hif-1a monomer and Hif1an, which blocks Hif-1 transcriptional activity by preventing Hif-1 association with p300 [83–88]. In this study we found that the transcript level of EGLN3 is decreased upon injury to the spinal cord. EGLN3 acts as the cellular oxygen sensor and is the most important enzyme in promoting Hif-1a degradation. This may explain why its down-regulation causes a positive regulation of the response to hypoxia. EGLN3 has other functions such as NGF-induced proapoptotic effect in neurons, probably through regulating CASP3 activity [89].

Hif-1a induction and activation under hypoxic condition induces NF-kB and its inhibitor at the same time [83, 86]. In this study, we found that NF-kB related transcripts were all up-regulated. For example, the transcript levels of NFKB2 and of the inhibitors NFKBIA, NFKBIE, NFKBIZ are all up-regulated during the first week after injury (data not shown). Another complication of disruption of blood supply is the phenomenon of ischemia/reperfusion injury causing necrotic injury to oligodendrocytes, neurons, astrocytes, and endothelial cells in the epicenter [90, 91]. This involves many events such as hypoxia, reactive oxygen species (ROS) and lipid peroxidation, cytokines, complement activation, and pro- and anti-apoptotic signaling cascades [91–93]. The ischemia/reperfusion injury is mainly under the regulation of the NF-kB signaling cascade and NF-kB transcription and its signaling cascade are, in turn, responsible for positive regulation of many immune-related responses, anti-apoptotic and equally important but opposing and controversial pro-apoptotic pathways [94].

Induced innate-immune response and Toll-like receptor signaling: a biphasic process

The inflammatory response to injury is initiated within minutes after SCI [95]. Our enrichment analysis scored inflammation as the most significant process starting immediately after injury and transcription activation of many immune-related genes. Many cytokines and chemokines are produced and secreted by various cells in the spinal cord tissue. It has been shown that IL-1B is produced immediately by astrocytes and neurons [95, 96]. Similar to other studies, our data indicates an up-regulation of IL-1B and TNF-alpha after injury. Most notably, we observed that the inflammatory response, in general, and specifically the cytokines’ expression pattern follow profile 44 (Figure 3F). Profile 44 represents the change in transcript levels of many genes with the first wave of up-regulation on day 1. The early up-regulation then disappears on day 3 and comes back to high levels for many transcripts from day 7 onwards. Such a phenomenon has been reported previously in mice with a contusion injury [95], although the cessation of primary up-regulation occurred after 24 hours and returned to an increased state on day 14 post-injury. In line with this observation, a biphasic model of cellular inflammatory response has been shown when various immune cells were analyzed using flow cytometry after SCI [97]. Various categories of processes are depicted in Figure 7A-P, which confirm such studies. Accordingly, we can extrapolate our findings and assume the same mechanism of expression or secretion for transcripts with the same profile of expression. This biphasic mode of expression was observed in other enriched GO biological processes such as activation of innate immune response, response to lipopolysaccharide, response to interferon-gamma, tumor necrosis factor superfamily, cytokine production, interleukin-6 production, interleukin-8 production, cytokine secretion, neutrophil chemotaxis, endocytosis and phagocytosis Toll-like receptor signaling pathways, integrin-mediated signaling pathway, T and B cell activation, and immunoglobulin-mediated immune response. The simplest explanation for this observed biphasic response is that the first wave of transcription activation of these genes originates primarily from neurons, astrocytes and microglia cells within the injured area of spinal cord, which subside by day 3. By day 7, post-injury immune cells such as neutrophils, macrophages, T and B lymphocytes have infiltrated the injured cord and amplify the production and secretion of related cytokines and chemokines as the secondary response tends to be at a higher magnitude.

The synthesis of IL-1B in neurons was shown to be dependent on NALP1 inflammasome [96]. In astrocytes, however, overexpression of inflammatory cytokines such as CCL2, CCL3, CXCL1 and CXCL2 is triggered by the IL-1 receptor and not the Toll-like receptor signaling proteins TLR2 and TLR4. Our data indicate that the mRNA levels of IL-1B, IL-1R2 and its accessory protein IL-1RAP, were up-regulated especially on day 3. Central proteins in the Toll-like receptor signaling such as TLR2, TLR4 and MYD88 were all up-regulated, the expression pattern of which follow profile 44. It has been shown that for neutrophils to enter the damaged zone in the spinal cord, the expression of IL-1R and MYD88 are essential [98]. Additionally, the cellular extravasation of neutrophils and other leukocytes into the injured area of spinal cord also requires up-regulation of matrix metalloproteinases (MMPs) [99, 100]. MMPs up-regulations are, in turn, dependent on the interaction of Fas and its ligand and on the peripheral myeloid cells and activation of Syk kinase to trigger recruitment to the injury sites [101, 102]. In our injury model, we observed an increase in the mRNA levels of MMP2, MMP9 and MMP12. We did not observe an increase in transcript levels of Fas or its ligand, but the Sky mRNA was up-regulated on day 1 and afterwards up to 8 weeks post-injury.

Toll-like receptor signaling is initiated after pattern recognition receptors (PRRs) detect pathogen-associated molecular patterns (PAMPs) or danger-associated molecular patterns (DAMPs), which are endogenously generated from tissue and cellular damage. It is now thought that for induction of innate immune response, two signals are required, the first from Toll-like receptors (TLRs) and the second from Nod-like receptors (NLRs). NLRs are responsible for processing of pro-interleukin-1B to IL-1B and pro-IL-18 to IL-18 [103]. Following injury to the spinal cord, processing of pro-IL-1B and pro-IL-18 into the mature form requires NALP1, ASC (PYCARD), CASP11, and finally CASP1 action to cleave the pro- forms [96, 98, 104]. Activation by endogenous signals in response to SCI seems to be the mechanism of activation of inflammation after SCI. We observed the up-regulation of the NOD1 component early after SCI. We also found that, after clip-injury to the spinal cord, PYCARD and CASP1 transcripts are highly up-regulated until 8 weeks post-injury as well as IL-1B and IL-18 transcripts. In addition, the expression of purinergic receptor P2X, ligand-gated ion channel 4 (P2RX4), which has been shown to regulate the inflammasome activation after spinal cord injury [105] was persistently increased in our injury model.

Adaptive immune response and antibody production

Both IL-1 and IL-18, produced during the first phase of inflammation mediated through the two-signal model of TLRs and NLRs, can induce the cellular and humoral modes of the adaptive immune response. IL-18 affects natural killer (NK) cells, monocytes, dendritic cells, T cells, and B cells, thereby regulating not only the innate, but also the adaptive immune responses [106]. Administration of IL-18 promotes production of interferon-gamma by natural killer (NK) as well as T cells. In our study we observe a late interferon-gamma response, which could be part of the second wave of cytokine production by T cells. T cell migration and activation precede the response to interferon-gamma, but other developing adaptive immune responses such as immunoglobulin-mediated immune response run in parallel to the response to interferon-gamma, which may explain the timing of the two processes (Figure 7).

It has been shown that autoantibodies are generated and detected in patients with chronic SCI [107, 108]. These detected antibodies can recognize a variety of related and unrelated antigens to CNS tissue. Mice defective in production of B cells, and thus antibody production, exhibit reduced pathological symptoms and improved locomotor recovery [109]. The activation of B cells has been shown to be level dependent as T3 injury completely abolished B cell response and T9 injury level induces significant B cell activation and antibody production [110]. It has been postulated that delayed antibody production and accumulation of autoantibodies leads to complement activation through C1q, which triggers the enzymatic cascade of the classic complement activation pathway and recruitment of microglia and macrophages to the site of injury [111]. Our data clearly show that a delayed adaptive immune response is initiated through immunoglobulin-mediated signaling and that this response is consistently and increasingly up-regulated towards the chronic phase in parallel to activation of the complement cascade (Figure 7C and 7N-P). However, the initial events such as T cell migration, T cell and B cell activation and proliferation starts very early after the injury (Figure 7N-P). As shown, the B cell-receptor signaling pathway seems to be a much more significant process than T cell receptor signaling (Figure 7N) which implies that, compared to the cellular T-cell mediated immune response, B cell-mediated immunity and neutralizing antibody production is the dominant immune response during the chronic phase of the injury to the spinal cord (Figure 7).

Animal care and thoracic spinal cord injury

All experimental protocols were approved by the animal care committee of the University Health Network in accordance with the policies established in the guide to the care and use of experimental animals prepared by the Canadian Council of Animal Care. Female Wistar rats (250 g; Charles River Laboratories, 4 sham and 4 injured animals for each time point) were used for this study. Injuries by the aneurysm clip method were performed as previously described [20, 29, 45, 112]. Briefly, under halothane anesthesia (1-2%) and a 1:1 mixture of O2/N2O, the surgical area was shaved and disinfected with 70% ethanol and betadine. A midline incision was made at the thoracic area (T4-T9), and skin and superficial muscles were retracted. Rats underwent a T6-T8 laminectomy and then received a 35 g clip (Walsh) moderate to severe compression injury at T7 for 1 min. The surgical wounds were sutured, and the animals were given Clavamox (Amoxicillin plus Clavulanic acid) for 7d and standard postoperative analgesia treatments and saline (0.9%; 5 ml) to prevent dehydration. Animals were allowed to recover and remained housed under standard condition for the duration of the experiment.

RNA isolation, processing and microarray hybridization

Rats were sacrificed at 1, 3, 7, 14 and 56 days after injury, and a 5 mm sample of the spinal cord containing the epicenter of the injured tissue was extracted for RNA analysis. Total RNA from each individual sample was extracted using TRIzol reagent (Invitrogen, Burlington, ON, Canada). RNeasy mini spin columns (QIAGEN, Mississauga, ON, Canada) were used for purification of total RNA molecules larger than 200 bp, which excludes smaller RNAs such as miRNAs. RNA quality was assessed with a 2100 Bioanalyzer (Agilent). cRNA for microarray hybridization was prepared from 5 ug of starting RNA using the protocol supplied by Affymetrix (Santa Clara, CA). cRNA was hybridized to GeneChip Rat Genome 230 2.0 arrays (24 chips total) at the Centre for Applied Genomics, The Hospital for Sick Children, Toronto, Canada). Primary data sets were saved in a MIAME-compliant format and uploaded to GEO (series GSE45006).

Microarray data analysis

Data analysis was performed in R with the Affy package (v1.12.2) [113] in BioConductor [114]. Data were investigated for spatial and distributional homogeneity. Normalization was performed with the sequence-specific GCRMA algorithm (package v2.6.0) in BioConductor [115]. Significance testing of this dataset was performed using linear models and pair-wise comparisons [116]. Each set of animals from a given time point was analyzed and pre-processed separately. The pre-processed data were then significance-tested using a linear modelling implemented in the limma package (v2.9.10) of BioConductor. Each sub-group was fitted to a separate factor in the design matrix, and the pair-wise contrast corresponding to differential expression of injured animals relative to control (sham) animals was extracted using a contrast matrix. Empirical Bayes moderation of the standard error [117] and false-discovery rate correction for multiple testing [118] were employed, again as implemented in the limma package. ProbeSets were deemed differentially expressed at p < 0.001 in any given comparison. Significantly different ProbeSets were visualized using the Heatplus package (v1.4.0) of Bioconductor. Euclidean distance was used as the distance metric for unsupervised hierarchical clustering using the DIANA algorithm with the cluster package (v1.11.4) in R (v2.4.1), and scaling was performed across rows. Clustering was used as a tool for replicate visualization and contrast comparison, not for gene selection [119].

The resulting gene set data with fold change and associated ANOVA t test p-values were analyzed by Short Time-Series Expression Miner (STEM) (discussed below), which allows the temporal expression patterns to be examined and extracted from the pool of up- and down-regulated transcripts across all time-points. Alternatively, individual time-point data were analyzed separately for up- and down-regulated genes, protein classes and signaling pathways. Both approaches were combined with functional analysis of transcripts using gene ontology (GO) enrichment.

Time-series expression profile clustering

We used the non-parametric clustering algorithm of STEM (Short Time-series Expression Miner, version 1.3.7) that is specifically designed to analyze short time-series expression data [120]. STEM implements a novel clustering method that can differentiate between real and random patterns and clusters genes by assigning them to a series of pre-defined patterns, named expression profiles. A profile is considered significant if the number of genes assigned to it exceeds the number of genes that are expected to occur by chance. The statistical significance of the number of genes assigned to each profile versus the expected number was computed and corrected for false discovery rate at p ≤ 0.05.

GO enrichment analysis

STEM is a statistical technique based on unsupervised clustering to find cluster-centroids followed by assignment of genes using distance-classifications, with statistical analysis using enrichment-based techniques. The biological significance of a set of genes can be assessed by GO enrichment analysis. Deregulated transcripts with ANOVA t test p-values ≤0.05 and fold change values > 1.5 were analyzed by the GO enrichment analysis module of STEM. Temporal analysis of the list of deregulated genes was performed using both time-series and time-point approaches. Due to more comprehensive gene coverage of RGD annotation data source file, the enrichment analysis was performed with reference to the RGD association file. For GO analysis of various expression profiles, we applied the annotations of “Biological Process” (BP) domain and the minimum expression fold change (in log₂ scale) was set to different values from zero. Other parameters were set to different values as follows: “minimum GO level” to different values from 3 to 20, “minimum number of genes” to 5, and “multiple hypothesis correction method for actual size based enrichment” to Bonferroni. STEM also offers to run the GO enrichment analysis at different GO tree levels, which allows limiting the results to more specific terms in the directed acyclic graph (DAG) structure of the gene ontology hierarchy. In this study, the time-point GO enrichment analysis was also employed to discover common up- and down-regulated biological processes across the time-points as well as possible unique processes to each time-point. The output GO terms were used for inter-relationship analysis and visualization by Venn diagram tool and or visualized as a scatter plot or interactive graph using REViGO [40].