Cluster wise analysis of different functional groups.

We intended to perform functional enrichment analysis on the clusters to study the functional groups that are largely activated during the process of regeneration. We have used Ingenuity Pathway Analysis (IPA) to identify functional groups. We provided both up and down-regulated genes of each cluster as input and selected 12 significantly enriched groups that are related to regeneration. All these groups are highly enriched in the clusters. Some of the enriched functions that we identified are ‘cell death’, ‘cellular growth and proliferation’, ‘cellular development’, ‘cell cycle’, ‘DNA replication and repair’, ‘nervous system development’ etc (Fig. 2A; Figure S1 and Table S2).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Figure 2. IPA analysis of biological functional groups and enriched pathways.

A) Pie charts showing the percentage of genes that are differentially expressed (both up and down regulated) at different time points for the various enriched functional groups in IPA analysis. B) Pie charts showing the percentage of genes that are differentially expressed (both up and down regulated) for the various enriched signaling pathways in IPA analysis at different time points after SCI.

https://doi.org/10.1371/journal.pone.0084212.g002

Functional group analysis in five different time points after injury.

The set of differentially expressed genes, both up and down-regulated, are also grouped on the basis of their highest expression at a particular time point during regeneration. For instance, when a differentially expressing gene shows highest expression at 3 dpi, it is assigned to the 3 dpi group. Similarly, the whole data is divided into 5 groups (1, 3, 7, 10 and 15 dpi groups) and each group is subjected to IPA analysis. We observed that same set of functional groups are significantly enriched in these day-wise groups (Fig. 2A; Figure S2 and Table S2). To analyze the data in more detail and relate information to the different phases of regeneration, we have analyzed the functional enrichments during the stages of regeneration (Fig. 2A). It is evident that different functional groups are enriched in different time points suggesting their involvement with particular event(s) in regeneration. Early time points like 1 and 3 dpi showed enrichment of similar functional groups like ‘cell death’, ‘cell growth and proliferation’ and ‘cellular development’. Whereas in 7 dpi, which is the proliferative phase of regeneration, other functional groups like ‘cell cycle’, ‘cell growth and proliferation’, ‘DNA replication and repair’ are significantly enriched. Similarly, in late regenerative phase that is related to tissue repatterning, enriched functional groups are ‘nervous system development and function’, ‘tissue development’ and ‘organ development’.

IPA analysis also allowed us to identify enrichment of genes in canonical pathways during the process of regeneration. Some of the important canonical pathways that are enriched at different time points of regeneration are Wnt, axonal guidance, interleukin, purine-pyrimidine biosynthesis, G-protein coupled receptor signaling, human embryonic stem cell pluripotency signaling and several other signaling pathways (Fig. 2B; Figure S3 and Table S3). Interleukin (IL) signaling molecules are highly enriched in very early and late stages of regeneration and involved all different signaling pathways like Il-8, Il-6, Il-22, Il-10 and Il-k. Molecules related to axonal guidance signaling are highly enriched in 3 dpi, 10 dpi and 15 dpi cord, where highest enrichment is observed in 3 dpi cord. Embryonic stem cell pluripotency is enriched at several time points 3 dpi, 7 dpi and 10 dpi cord. We find a high predominance of molecules related to purine and pyrimidine metabolism in 7 dpi, when cell proliferation is at its highest peak. Similarly, cAMP mediated signaling is enriched in late time points, in 10 dpi and 15 dpi cord where axonal regeneration is initiated and maintained. Interestingly we observed that many signaling pathways are enriched in the initial time points (Day 1 and Day 3) and also in late time points (Day 10 and Day 15). Another salient observation is enrichment of few but very interesting pathways like ‘One Carbon Pool by Folate’ and ‘N-Glycan Biosynthesis’ at 7 dpi cord (Figure S3).

Control of inflammatory response.

Inflammatory response is known to be controlled by members of cytokine, chemokine signaling pathways in mammalian SCI and other diseased state [17]–[19]. Here we found 27 differentially expressed genes (Fig. 3A, Table S4), among them, 21 genes were expressed their maximum level at 1 dpi or 3 dpi cord. Some of the differentially expressed important cytokines and chemokines are ccl1, ccrl1a, cmklr1, crfb2, crfb8, cxcl12b, cxcr3.2, socs3a, socs3b, il1b, il4r, il22, irf10, irf11, irf8 and irf9. Molecules like il4r, ifnα, tgfβ are upregulated both in mammalian CNS [20] and in zebrafish spinal cord after injury as studied by us, but their temporal pattern and level of expression varies. For example, tgfβ1 in mammalian brain after hypoxic injury shows 0.54 fold changes in 8 hr compared to 4 fold changes in 1 dpi zebrafish cord. On the contrary, il10, which is a neuroprotective interleukin, is downregulated in mammalian SCI [21] but upregulated in 1 dpi zebrafish cord (il10 data not shown due to low fold change 1.36). The tgfβ1 and il10 play critical role in suppressing immune response [22], [23]. We observed that socs3, stat3 and tgfβ1 are highly upregulated in 1 dpi zebrafish cord that is in early phase of regeneration (Fig. 3 and 4, Table S4).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Figure 3. Differentially expressed genes related to inflammatory response and validation of stat3 expression in uninjured and injured zebrafish spinal cord.

A) Dendrogram representing genes related to regulation of inflammation. Each horizontal line indicates the expression pattern of each gene and the vertical columns indicate the uninjured control and time points after SCI. The color chart indicates mean fold change of gene expression in each time points. Red and green colors represent increased and decreased expression respectively. B) Quantitative RT-PCR of stat3 expression showing fold change (red graph) and pattern of expression at different time points after injury. The temporal expression pattern of qRT-PCR (red graph) was compared with results of microarray analysis (blue graph). Error bars represent SEM, n = 3, p<0.01. C) A longitudinal section of uninjured cord stained with STAT3 and DAPI. D) A longitudinal section of 1 dpi cord showing many STAT3 positive cells (white arrows) close to the injury epicenter (double white stars) counter stained with DAPI. E–H) A longitudinal section of 3 dpi cord double stained with STAT3 and radial glial marker- GFAP, showing individual panels of STAT3 (E), GFAP (F), merge with DAPI (G) and DIC (Differential interference contrast) (H). A few STAT3⁺ cells are colocalized with GFAP⁺ cells (white arrowheads) close to injury epicenter. I–L) A longitudinal section of 3 dpi cord showing colocalization of STAT3 and newly formed neuronal marker, HuC/D with individual panels of STAT3 (I), HuC/D (J), merge with DAPI (K) and DIC (L). Some of the STAT3⁺ cells are colocalized with HuC/D⁺ cells (yellow arrowheads) close to the injury epicenter. M–P) A longitudinal section of 3 dpi cord colocalized with STAT3 and proliferating cell marker, BrdU shows individual panels of STAT3 (M), BrdU (N), merge (O) and DIC (P). Proliferating BrdU⁺ cells are also colocalized with STAT3⁺ cells (blue arrowheads) close to the injury epicenter. White dashed line in panel C to P mark the boundary of spinal cord tissue. Q–T) A representative higher magnification of 3 dpi cord section shows individual panels like STAT3 (Q), GFAP (R), DAPI (S) and merge (T). A STAT3⁺ cell (white arrowhead; the boundary of the cell nucleus is markedby white dashed line) is colocalized with GFAP⁺ radial glia (white arrowhead; the boundary of the cytoplasm of radial glia is marked by white dashed line). U–X) A representative higher magnification of 3 dpi cord section shows individual panels like STAT3 (U), HuC/D (V), DAPI (W) and merge (X). A STAT3⁺ cell (yellow arrowhead; the boundary of the cell nucleus is marked by white dashed line) is colocalized with HuC/D⁺ newly formed neuron close to ependyma (yellow arrowhead; the boundary of the cell nucleus is marked by white dashed line). Y) Quantitative analysis of STAT3⁺, STAT3⁺/GFAP⁺, STAT3⁺/HuC/D⁺ and STAT3⁺/BrdU⁺ cells in 3 dpi cord in longitudinal sections. The value represented as Mean±SEM of individual longitudinal section, n = 5 cord, p<0.01.Scale bar = 200 µm (C–D); 50 µm (E–P), 20 µm (Q–X).

https://doi.org/10.1371/journal.pone.0084212.g003

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Figure 4. TGFβ expression in uninjured and injured zebrafish spinal cord.

A) A longitudinal section of 1 dpi cord showing TGFβ positive cells in both grey matter (GM, yellow arrows) and white matter (WM, white arrows) of the cord at injury epicenter (double white star) and adjacent part. B) Same section merged with DAPI. C) Same section merged with DIC. D–G) Higher magnification of the boxed area of the section in B showing many TGFβ positive cells in both grey matter (GM, yellow arrows) and white matter (WM, white arrows) regions of the cord represented in individual panels of TGFβ (D), DAPI (E), merge of both (F) and merge with DIC (G). H) Quantitative RT-PCR of tgbβ1 expression showing fold change (red graph) and pattern of expression at different time points after injury. The temporal expression pattern of qRT-PCR (red graph) was compared with results of microarray analysis (blue graph). Error bars represent SEM, n = 3, p<0.05. I) Quantification of TGFβ expression in uninjured and injured cord at different time points by ELISA. Error bars represent SEM, n = 5, p<0.01.Scale bar = 50 µm (A–C); 10 µm (D–G).

https://doi.org/10.1371/journal.pone.0084212.g004

The array data of stat3 and tgfβ1 molecules have been further validated by qRT-PCR, immunohistochemistry and ELISA analysis (Fig. 3 and 4). Both stat3 and tgfβ1 show highest upregulation in 1 dpi cord and array hybridization data is confirmed by qRT-PCR, immunohistochemistry and ELISA analysis. STAT3 expression can be localized in injured cord, predominantly in the grey matter both in ependymal cells and in neuronal cells. It is important to mention that expression was observed exclusively at the injury site of 1 dpi and 3 dpi cord (Fig. 3D, E, I, M) and absent in the normal part of same cord. Colocalization study of STAT3 with GFAP (Fig. 3E–H, Q–T) and HuC/D (Fig. 3I–L, U–X) in 3 dpi cord suggests expression both in radial glia and in neurons, although we found that a higher number of STAT3 positive cells are HuC/D positive newly formed neurons [6] (Fig. 3K, X, Y; Table S5). We have also analyzed STAT3 colocalization with BrdU in 3 dpi cord, which is the beginning of proliferation stage and found that many STAT3⁺ cells are BrdU⁺ proliferating cells (Fig. 3M–P). There is also a population of STAT3 positive cells which are not colocalized with HuC/D and GFAP and they are probably inflammatory cells present in the injury epicenter of both 1 and 3 dpi cord. Our expression analysis of tgfβ1 expression showed highest level of expression in 1 dpi cord and it was present in both grey and white matter of adjacent part of the epicenter (Fig. 4A–G). TGFβ, a multi-potent cytokine, plays critical role in suppressing the immune response after CNS injury in mammals [23]–[26]. Furthermore, it also plays a role in epithelial-mesenchymal transition (EMT) after SCI, expressed in M2 type macrophages [27] and may induce initiation of cell proliferation during tail regeneration in Xenopus [28].

Other analysis in mammalian CNS injury refers to importance of complement factors and MHC class I and II molecules expression [20], [29]. We observed upregulation of two complement factors c1qc and c1a13b early after injury (Table S4). Between these two factors prominent role of c1qc in activation of phagocytosis of macrophage/microglia in mammalian CNS is well known [30] which may suggest a similar role in zebrafish SCI where we observe increased phagocytosis by macrophage in 3 dpi cord [6]. Unlike mammalian CNS injury, where both MHC class I and II molecules are upregulated, we only found upregulation of MHC class I molecules like mhc1uba, mhc1uaa and mhc1ze in early time points in zebrafish cord regeneration (Table S4).

Pro- and anti-inflammatory macrophages.

Our previous study highlighted importance of macrophage during early phase of regeneration in zebrafish spinal cord. In mammalian SCI presence of both pro- and anti- inflammatory macrophages are related to neuroprotective and neurodegenerative response [27]. In this context, it is worth mentioning that M1 macrophages are pro-inflammatory as characterized by pro-inflammatory cytokine release [31] and M2 macrophages as characterized by presence of scavenger receptors, IL-1 receptor antagonist [32] resulting in decreased production of pro-inflammatory cytokines like il-1β and upregulation of anti-inflammatory cytokines like il-4 and il-13 [27], [33]. We observed presence of a higher number of M2 type macrophage related genes (anti-inflammatory) compared to M1 type (pro-inflammatory) (Table S6). M2 type macrophage related genes include scarb1, scarb2, il4r, tgfβ1, vegfa, tgif1, arg2 and many of them are upregulated in 1 dpi and 3 dpi cord.. Only two of the genes related to M1 type macrophage (caspa and nos2) are expressed. Suppression of nos2 limit self-damage of phagocytic macrophages [34]. Similarly, we also found nos2 expression as a consequence of inflammatory response.

Activated macrophage are known to play a role in wound healing, allergic reaction and proliferation and arginase activity controls the role of activated macrophages after injury [35]–[37]. Inflammatory response is complex as we observe both pro- and anti-inflammatory types of molecules, although a bias towards anti-inflammatory response was observed as we see strong upregulation of many M2 type macrophage related molecules. Macrophage may also enhance axonal growth by secreting several molecules like cytokines and neurotrophic factors and allowing Schwann cell migration as it happens in peripheral nervous system (PNS) [38]. Macrophage along with Schwann cell helps in clearing myelin and axonal debris as reported by others and us [6], [39]–[41]. Macrophage may also induce MMP-9, which is playing important role in debris clearance and preventing collagen scar formation [42], [43]. Our present data highlight the importance of macrophage as a key player in controlling inflammation, debris clearance and creation of permissive environment for axonal regrowth during regeneration.

Cell death.

In mammalian SCI cell death contributes to huge neuronal and glial loss and involves distinct necrotic and apoptotic pathways. Apoptotic cell death is common to both mammalian and zebrafish SCI but in a different spatio-temporal pattern [6], [44]. We have identified many cell death related genes, their regulators and anti-apoptotic genes with differential expression during regeneration (Table S7). Two members, tnfα induced protein and fasl, which are upregulated in different time points and may have dual role in controlling apoptotic and necrotic pathway. Apoptotic cell death in mammalian CNS injury is regulated by Caspase and Calpain family of genes [45], [46] and we observed upregulation of caspa, caspb, casp7, capn9 in 1 dpi and 3 dpi zebrafish cord. Upregulation of 12 apoptotic genes in early time points like 1 dpi and 3 dpi are probably induced as a consequence of injury, as these genes are absent in control, and are significantly downregulated in 7 dpi cord. Among these 12 genes, capn9 and traf4a are upregulated in 10 dpi and 15 dpi cord. Significance of expression of these molecules in late time points during regeneration is not well understood although these are implicated in zebrafish CNS development [47]. Mammalian anti-apoptotic protein Bcl-X_L is known to control cell death and cell survival [48]. We also observed that many anti-apoptotic genes like bcl2113, mcm5, bag3, tnfaip6 and igf1ra [49]–[52] are upregulated in different time points after injury.

Cell migration and epithelial-mesenchymal transition.

Cell migration is one of the key phenomena in early regeneration time points as we showed that accumulation of cells in the injury epicenter occur as early as Day 3 and continues at least until Day 7 post injury [6]. Molecules like fibronectin, integrin, laminin, cadherin, matrix metalloproteinase(s) (MMPs) and collagen are known to be involved in cell adhesion, cell migration in developing and regenerating organs [5], [14], [15], [53]–[57]. Our data showed that most of the genes related to cell migration are up-regulated either in 1 dpi (fn1b, mmp13, mmp9) or in 3 dpi (itgav, vim, itgb4, itgb3a, cdh2, cdh23, lamb4, lamb1, cxcl12b) cord (Table S8). Many genes are also upregulated in different regenerating organs probably indicate their common role in tissue remodeling [15], [54]–[56]. Epithelial-mesenchymal transition (EMT) or interaction during regeneration is well known [54], [58], [59]. A similar EMT is probably operating during early phases of regenerating zebrafish spinal cord although it has not been characterized yet. We see expression of 8 EMT related genes which are upregulated in zebrafish cord after SCI (Table S8). A separate group of six genes represents Adam and Integrins of which adam8a and adam9 show their upregulation in early time point (Table S8). Genes like cdh2, cdh23, lamb4, lamb1, cxcl12b are highly upregulated in 3 dpi cord. Further upregulation of these genes at 10 dpi or at 15 dpi cord, may suggest their role in different types of cell migration. In this context, it is important to refer that both neuronal migration and Schwann cell migration [60], [61] is known to occur at later time points after injury.

Epithelial-mesnchymal (EM) interaction is a phenomenon during dedifferentiation phase. We see expression of some of the genes related to dedifferentiation process (like msx-c, msx-e and vim) in 3 dpi zebrafish cord (Table S9).

Proliferative response after injury in spinal cord.

Injury induced proliferation has been documented previously and proliferation is necessary to replace the lost neural cells. The proliferative response is controlled by many genes, which have been categorized as genes related to cell cycle and other regulators, involved in proliferation (Fig. 5; Table S10). As many as 48 genes associated with cell cycle are differentially regulated. Interestingly almost all genes expressed in uninjured control are also represented in 7 dpi cord, but are not expressed in other injury time points. Exception among these patterns is hsp90a.1 and hsp90a.2, which are upregulated in uninjured cord and 1 dpi but down regulated in other injury time points (Fig. 5A). A high number of genes are Cyclins and Cdc/Cdks (Fig. 5A; Table S10). Based on cell counts of colocalized BrdU/H3P along with DAPI indicate that in uninjured cord there are a low number of cells in S-phase (2%) and M-phase (1%) compared to Go-G1 phase (97%). Number of cells in both S-phase and M-phase, increases after injury, although percentage of cells entering into S-phase are higher in both 3 dpi and 7 dpi cord (8% and 12% respectively) when compared to percentage of cells in M-phase (2% and 5% respectively) (Figure S4). Our present data on array analysis corroborated with the cell count data and showed that genes involved in G1-S phase transition are selectively upregulated either in 3 dpi (ccnd1, ccni and myca) and 7 dpi (cdk2, cdk7, ccne and ccnh) cord and all 3 genes associated with S-phase (ccna2, pcna, uhrf1) are upregulated only in 7 dpi cord when highest number of proliferating cells are in S-phase. Both the uninjured and injured cord showed expression of cdk2, ccnh, ccne and ccna2, but fold change varies, expression level is higher in injured than uninjured cord. Among 12 genes involved in G2-M transition, 10 genes (ccnb1, ccnb2, cdc20, klf11a, mcm6, mcm2, mad2l1, ttk, plk1, kifc1) which are M-phase related are expressed in both uninjured control and higher in 7 dpi cord. Expression of four genes have been further validated by qRT-PCR which confirms upregulation of ccnd1 within 1 dpi to 3 dpi cord and ccne, cdk2 and ccnb1 in 7 dpi cord (Fig. 5C–F). Cell cycle regulators like myca and ccnd1 both are down regulated in uninjured cord but upregulated early after injury and myca is known to upregulate its downstream molecules ccnd1 and cdk2 [62]. Cell cycle regulators like cdc20 have a role in degradation of ccnb1 during mitosis [63]; both cdc20 and ccnb1 are expressed in uninjured cord and highly upregulated in 7 dpi cord and subsequent downregulation of ccnb1 in 10–15 dpi cord probably refers to reduction of cell proliferation and down regulation of mitosis markers at late stage of regeneration [6]. Cell cycle regulators also include members of Mcm families, members of Hsp families, plk1, ttk, klf11a, ppp1r3b and ppp1r7 (approximately 19 genes) are upregulated both in uninjured control and 7 dpi cord and many of them are M-phase regulators.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Figure 5. Differentially expressed genes related to cell proliferation in zebrafish spinal cord after injury and validation of cell cycle regulators.

A) Dendrogram representing a group of cell cycle regulators directly involved with different cell cycle phases. B) Same representing group of genes indirectly control cell proliferation which include both negative and positive regulators. Each horizontal line indicates the expression pattern of each gene and the vertical columns indicate the uninjured control and time points after SCI. The color chart indicates mean fold change of gene expression in each time points. Red and green colors represent increased and decreased expression respectively. C–F) Quantitative RT-PCR analysis of ccnd1 (C), cdk2 (D), ccne (E) and ccnb1 (F) showing fold change and pattern of expression at different time points after injury. The temporal expression pattern of qRT-PCR (red graph) was compared with results of microarray analysis (blue graph). Error bars represent SEM, n = 3, p<0.05.

https://doi.org/10.1371/journal.pone.0084212.g005

As mentioned earlier that proliferative response are being controlled by cell cycle regulators along with a second category of genes, which include cell proliferation related genes (Fig. 5B; Table S10). 30 genes are further separated into four different subgroups based on their temporal expression pattern and known function. There are 11 genes, which are known to be negative regulators of cell cycle progression, mostly upregulated in 1 dpi and 3 dpi cords. Among them anxa1a, cebpd and cebpb are all upregulated in 1 dpi cord and may play a role in arresting of G1-S transition [64], [65]. Another important functional group of genes involved in positive regulation of cell cycle progression includes 12 genes. Atf-3 is one of the highly upregulated genes in 1 dpi and 3 dpi cord compared to uninjured cord. A similar elevated expression of atf-3 was demonstrated after injury in optic nerve [66], in dorsal root ganglia (DRG) and in motor neurons [67].

Differentially expressed genes involved in neurogenesis and neuronal differentiation.

Many genes related to neurogenesis are differentially regulated during the process of regeneration. Previously we showed that generation of new neurons is related to cell proliferation and Hu⁺/Brdu⁺ cells are continuously present in 3 dpi, 7 dpi and 10 dpi cord. Some of these cells become NeuroD positive and have undergone full neuronal differentiation [6]. Molecular analysis showed at least 54 genes are differentially expressed (Fig. 6Q; Table S11) and among them 41 genes are transcription factors (Figure S5). All the 54 genes are categorized in different groups based on their temporal pattern of expression. Group-I includes her2, dab2, pou5f1, emx3, bmi, paxip, sox19b and sox21a that are expressed in uninjured control and are upregulated in 7 dpi cord, when we also observe a very high level of proliferation, which contributes significantly to neurogenesis [6]. Among these genes, expressions of her2 and dab2 have been reported during retinal neurogenesis and in proliferating neuroepithelium of developing CNS [68], [69]. Group-II consists of large number of genes (36 genes) that are upregulated in 3 dpi, 10 dpi and 15 dpi cord and down regulated in uninjured control and 7 dpi cord. These are injury induced since these are not expressed in uninjured control. Genes like notch1a, jagg1a, deltab and deltad are all upregulated in 3 dpi and 10 dpi cord whereas proneural genes like neurog1, neurod2, neurod4 and olig2 are all upregulated in 10 dpi cord. Many of the genes in this group are members of Notch signaling pathway. Similar to developing CNS the above mentioned genes may promote differentiation of progenitors selectively to different neural fate and maintenance of reserve pools. Group-III consists of 10 genes neurod2, nrarpa, hes5, nkx3.2, ascl1a, numb1, hesx, shhb, bdnf and gadd45a all of these are highly upregulated in 1 dpi but downregulated in uninjured cord.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Figure 6. Differentially expressed genes involved in neurogenesis and neuronal differentiation after injury in zebrafish spinal cord and validation of Sox2 expression by immunohistochemistry and qRT-PCR analysis.

A–C) A cross section of uninjured cord showing SOX2 positive cells (white arrowhead) in ependyma (ep) around the central canal (cc), counter stained with DAPI (B) and merge (C). D) A cross section of 3 dpi cord showing higher number of SOX2 positive cells compared to uninjured (A) in grey matter (white arrowhead) mostly in ependyma (ep) around the central canal (cc). E–F) Same 3 dpi cord section counter stained with DAPI (E) and stained with radial glial marker GFAP (F). G) Same 3 dpi cord section showing both SOX2 positive cells (arrowheads) and their colocalization with GFAP (arrow). H) Same section (G) in higher magnification showing colocalized SOX2 and GFAP positive cells (arrow) around the central canal (cc), I) DIC of same section in H. J–L) 3 dpi cord section showing SOX2 positive cells (J, arrowheads) and BrdU positive cells (K, yellow arrow) in the grey matter and many of them are SOX2⁺/BrdU⁺ (L, yellow arrowheads) in ependyma (ep) around the central canal (cc). M–O) A 3 dpi cord section showing many HuC/D positive neuronal cells (M, red arrow) and SOX2 positive cells (N, white arrowheads) in the grey matter. Same section showing colocalization of SOX2 and HuC/D (O, red arrowheads). P) qRT-PCR of sox2 in uninjured cord and in injured cord at different time points. Error bars represent SEM, n = 3, p<0.05. Q) Dendrogram showing differentially expressed genes involved in neurogenesis and neuronal differentiation after injury. Each horizontal line indicates the expression pattern of each gene and the vertical columns indicate the uninjured control and time points after SCI. The color chart indicates mean fold change of gene expression in each time points. Red and green colors represent increased and decreased expression respectively. Scale bar = 50 µm (A–G and M–O), 20 µm (H–L).

https://doi.org/10.1371/journal.pone.0084212.g006

Many members of Sox gene family are differentially controlled during neurogenesis. There are two distinct pattern of expression within Sox family members. Sox19b and sox21a are upregulated both in uninjured cord and in 7 dpi cord. Others Sox family members (sox2, sox4a, sox9a, sox9b, sox10, sox11a, sox14, sox21b and sox32) are injury induced since they are upregulated in 3 dpi, 10 dpi and 15 dpi cord but down regulated in uninjured and 7 dpi cord. The highest expression values are different for different sox genes suggesting multiple functions during regeneration. Sox2, one of the key transcription factors is expressed in the neuroepithelium of the developing mammalian CNS and adult zebrafish brain. Its function in neural stem cell (NSC) maintenance, proliferation and specifying their identity [70]–[73] is well known. Our real-time PCR analysis showed that there is upregulation of sox2 after injury at Day 3 compared to uninjured cord (Fig. 6P). The immunohistochemical analysis showed that in uninjured cord SOX2 positive cells are present both in ventricular region and in subependymal region of grey matter (Fig. 6A–C). A higher number of SOX2 positive cells have been found in the grey matter of 3 dpi cord (Fig. 6D–H) compared to uninjured cord. These cells are situated mostly in ependyma and few are in subependyma. A subpopulation of SOX2 positive cells are radial glia since they co-localize with GFAP with appropriate morphology (Fig. 6G–H; Table S12). When we co-localize SOX2 positive cells with a neuronal marker (HuC/D), we found that many newly formed HuC/D⁺ neuronal cells (as they are morphologically smaller than specified mature neuron) in the grey matter are SOX2 positive (Fig. 6M–O; Table S12). The colocalization of SOX2 with BrdU identified many cells present in both ventricular and subventricular region (Fig. 6J–L; Table S12) thus confirming presence of proliferating neural progenitors in regenerating cord.

Differentially expressed genes involved in repatterning.

Regeneration is a complex process of rebuilding diverse tissue that involves patterning process and the mechanism substantially represents recapitulation of developmental pattern formation. The repatterning event has been studied in amphibian limb and spinal cord [74]–[77]. We have identified 84 pattern-forming genes that are differentially regulated (Fig. 7; Figure S6 and Table S13). Based on the temporal expression pattern we have separated 84 genes in 4 different groups. Highest numbers of genes (63 genes) are upregulated in 3 dpi, 10 dpi and 15 dpi cord that are downregulated in uninjured cord. Another category of genes (8) is upregulated in 1 dpi cord and gradually downregulated in later time points. All these genes are also down regulated in uninjured cord. The third and fourth groups, which include eight and five genes respectively, are upregulated in 7 dpi and 10/15 dpi cord. Genes belonging to these two groups are also upregulated in uninjured control unlike the first two categories (Figure S6 and Table S13).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Figure 7. Differential expression of pattern forming genes during zebrafish spinal cord regeneration and validations of pax6a, shhb and dbx2.

A) Represents known genes related to anterio-posterior patterning in CNS that are differentially expressed at different time points after SCI. The color chart indicates mean fold change of gene expression in each time points. Red and green colors represent increased and decreased expression respectively. B–D) Quantitative RT-PCR analysis of pax6a, shhb and dbx2 respectively showing fold change and differential pattern of expression in various time points after injury and compared with the temporal expression patterns of microarray data. Error bars represent SEM, n = 3, p<0.01. E–F) Quantitative analysis of PAX6 and SHH expression by ELISA in uninjured and injured cord at different time points. Error bars represent SEM, n = 5, p<0.001. G) Represents known genes related to dorso-ventral patterning in CNS that are differentially expressed at different time points after SCI. The color chart indicates mean fold change of gene expression in each time points. Red and green colors represent increased and decreased expression respectively. H–K) In situ hybridization of dbx2 in uninjured (H), 3 dpi (I) and 10 dpi (J) cord section shows presence of dbx2 transcripts in few neurons (black arrowheads) in uninjured cord but complete absence of dbx2 transcripts in 3 dpi cord and reappearance of dbx2 transcripts in many neurons in 10 dpi cord (black arrowheads). Significantly, higher numbers of dbx2 positive cells are present in the ventral side of the 10 dpi cord compared to uninjured cord. Yellow arrows indicate neurons of ventral progenitor domain. Section in K represents sense control of dbx2 in a 10 dpi cord. The mark ‘cc’ denotes central canal of the cord. ‘D’ and ‘V’ indicates dorsal and ventral side of the cord respectively. L) Representative transverse section of 10 dpi cord stained with HuC/D and DAPI showing the distribution of regenerated neurons in the dorso-ventral axis of the cord. Yellow arrows indicate neurons of ventral progenitor domain. M) A schematic diagram of various progenitor domains in dorso-ventral axis of spinal cord on the basis of spatial distribution and gene expression, adapted from Wilson and Maden [81]. ‘D’ and ‘V’ indicates dorsal and ventral side of the cord respectively. RP = Roof plate, FP = Floor plate, NC = Notochord, DP = Dorsal progenitor, VP = Ventral progenitor. Scale bar = 50 µm.

https://doi.org/10.1371/journal.pone.0084212.g007

Based on the studies on different species, combinatorial expression of different genes involved in dorso-ventral patterning of spinal cord have been identified [78]–[81]. Most importantly gradient of Shh induces five distinct classes of ventral neurons namely V₃, V_MN, V₂, V₁, V₀ from neural progenitors [78], [80] and ventral cell types are determined by expression of different transcription factors. A gradient of SHH-N establishes a dorso-ventral domain by modulating expression of transcription factors of class-I and class-II proteins. SHH gradient is ventrally high and dorsally low, it represses class-I proteins creating an opposite gradient of pax6, pax7, dbx1, dbx2, irx. Similarly same SHH gradient activates class-II proteins namely nkx6.1 and nkx2.2 [79], [80]. Our present array analysis data indicates that upregulation of shh during regeneration probably creating gradient similar to development. Most of the genes involved in specification of different neuronal subtypes in ventral cord (Fig. 7G; Figure S6 and Table S13). For example dbx2, irx3 and pax6, all involved in V₀ and V₁ neuronal patterning, are upregulated in regenerating cord similar to developing zebrafish spinal cord [82]. Similarly, olig2, pax6, hlxb9 and isl2, involved in motor neuron progenitor domain, are also upregulated in regenerating cord. The neuronal populations identified in the dorsal domain is specified by expression of Mash, Math, Neurogenin and LIM homeobox proteins in developing vertebrate cord [81]. In regenerating cord we also observed upregulation of neurog1, lim1 and lim3b (Fig. 7). We have validated three genes dbx2, pax6a and shhb by ELISA (Fig. 7E–F), qRT-PCR (Fig. 7B–D) and in situ hybridization analysis (Fig. 7H–K). Our qRT-PCR analysis revealed highest upregulation of all three genes in 10 dpi cord and significantly down regulated in control (Fig. 7B–D) and it suggests involvement of these genes in late events like repatterning of spinal cord. ELISA results of PAX6 and SHH confirms the expression pattern observed in array hybridization and qRT-PCR. Dbx2 is another transcription factor involved in patterning of ventral spinal cord and specification of V₀ and V₁ domain interneurons during development [81]. Our array analysis and qRT-PCR data showed up-regulation of dbx2 in 10 dpi and 15 dpi cord (fold change 8 and 6 respectively) and down regulation in other injured time points and uninjured control (Fig. 7D). In situ hybridization data reveals presence of few dbx2 positive cells in the grey matter of uninjured cord and whereas 3 dpi cord showed complete absence of dbx2 positive cells. The numbers of dbx2 positive cells have increased significantly in 10 dpi cord compared to uninjured cord (Fig. 7H and J). Dbx2 positive cells can be localized in the grey matter more ventrally, suggesting its possible role in respecifying ventral neuron during regeneration.

We have identified 52 genes, which are thought to be related to A–P patterning during CNS development in different species [83]–[85]. These genes are again differentially regulated at different time points after SCI. Based on the temporal expression pattern we observed three different groups. Majority of the genes (38 genes) showed upregulation in 3 dpi, 10 dpi and 15 dpi cord (21 genes are showing highest expression in 10 dpi and 15 dpi) but downregulated in uninjured control. Another set of genes showed upregulation both in uninjured control and 7 dpi cord but down regulated in other time points (Fig. 7A; Figure S6 and Table S13). A high number of genes (approximately 24) implicated in conferring A–P polarity of CNS belongs to four hox clusters. In regenerating zebrafish cord we observed six genes, viz. hoxb7a, hoxb8a, hoxc8a, hoxd10a, hoxa11b, hoxd11a which are all upregulated in 3 dpi, 10 dpi and 15 dpi cord. Upregulation of above mentioned hox genes in injured cord at later time points might suggest that they are probably involved in respecification of caudal spinal cord during regeneration. Furthermore, other patterning genes like cdx1a and gdf11 are also upregulated in regenerating 3 dpi cord which are known to be involved in caudal spinal cord patterning during development [86], [87]. Thus, different hox genes with overlapping rostro-caudal domain of expression in the developing cord establish different region along the A–P axis [88], [89]. We inflicted injury at approximately 15^th–16^th vertebrae level of adult spinal cord that is more caudal and any repatterning during regeneration would require respecification of spinal cord by combinatorial expression of relevant hox codes. More experimental validations involving different hox genes family members need to be done to corroborate our expression analysis.

Differentially expressed genes involved in axonal regrowth and guidance.

Time course analysis of regeneration demonstrated axonal regrowth and substantial functional recovery after SCI in zebrafish [4], [6], [10]. We can see reconnected axons at 15 dpi cord although molecular events related to initiation of axonal regrowth may start early [6]. In the present study, we have identified many genes involved in axonogenesis and axonal guidance. There are at least 80 genes involved in axonal regrowth and guidance that also include many transcription factors (Figure S7 and Table S14). Based on temporal expression patterns, we have found 3 different groups. Group-I consists of 11 genes and they are upregulated both in uninjured control and 7 dpi cord. Group-II consists of 13 genes and they are particularly downregulated in control but upregulated in 1 dpi cord and gradually downregulated in later injury time points. In Group-III, highest number (56) of genes showed upregulation in 3 dpi, 10 dpi and 15 dpi cord but downregulated in uninjured control. It is important to note that IPA analysis also showed high enrichment of axonal guidance signaling genes in all three above mentioned time points (Fig. 2; Figure S2 and Table S3). Similar to developing nervous system there is also a complex mechanism of axonal pathfinding during axonal regrowth that include different axonal guidance molecules with attractive, repulsive and adhesive properties [90]. For example molecules like robo1, robo2, slit1b, slit3, sema3ab and sema3h all are with repellent properties and are all upregulated in 3 dpi, 10/15 dpi cord. Axonal attractant molecules like netrin1a, netrin1b, plexina4, are also upregulated in 3 dpi and 10 dpi cord. Cell adhesion molecules like ncam show slightly different expression pattern, as they are upregulated in 3 dpi and 10/15 dpi cord. Ephrins and receptor tyrosine kinases are known to be involved in organization of axons in different topographic planes during development of CNS [91], [92]. Efnb1 and efnb3 are upregulated in 3 dpi and 10 dpi cord whereas rtk6 are upregulated in 1 dpi and 3 dpi cord but all remained downregulated in uninjured control. In developing zebrafish, motor guidance and target recognition molecules include CAMs [85], semaphorins [93], [94], netrins [93], robo, slit [95], ephrin and rtk [91]. Another gene atf-3 is also highly upregulated in 1 dpi, 3 dpi and 10 dpi cord and may be involved in axonal regrowth as indicated in various reports [66], [96]. Induction of atf-3 in both PNS injury and optic nerve injury is implicated in creation of permissive environment for axonal regeneration. Interestingly in zebrafish SCI, there is also creation of PNS like environment in CNS. As many as six genes (like atf3, fgf2, mmp-9, stat3, socs3 and sox11), which are related to PNS regeneration, are upregulated after injury. Similar to the previous studies on axonal regeneration in zebrafish CNS [5], [97] contactin1b and gefiltin are also upregulated in either 10/15 day after injury. All these molecules are upregulated in regenerating cord during axonogenesis suggesting similar mechanism of axonal growth and guidance are probably being employed in directing axonal regrowth. In adult, functional analysis of identified genes during axonogenesis would provide some clue to understand regrowth of axonal tracts in regenerating cord.

Wnt signaling.

Many members of Wnt signaling pathway are differentially regulated at different time points suggesting their diverse role in cellular events of regeneration (Fig. 8H; Figure S8 and Table S15). Several Wnt members like wnt4a, wnt4b, wnt7a, wnt8a, wnt9a, wnt11 and β-catenin are upregulated at different time points after injury. Other members of Wnt signaling pathway represent Frizzled like (fzd2, fzd3i and celsr2; zebrafish homolog of mouse flamingo) which are mostly upregulated in 3 dpi and 10 dpi cord. Both wnt4a and wnt9a are upregulated in 3 dpi and either in 10 and or 15 dpi cord and downregulated in other time points of regeneration. On the other hand, wnt8a showed upregulation in control and 7 dpi cord and down regulated in other time points. Regulators of Wnt signaling genes like dixdc1, gsk3b, tcf711 and tcf712 are all differentially upregulated in 3 dpi, 10 dpi and 15 dpi cord and down regulated in 7 dpi cord. There are different inducible targets of different Wnt, β-catenin or Frizzled, which include myca, mycn, myc11a, jun, ccnd1, retinoic acid receptors and many members of sox family genes. Role of different Wnts like wnt3a, wnt5a, wnt5b, wnt7a, wnt8a, wnt10a are well documented in different regenerating tissues [56], [98]–[101]. In regenerating fin, over expression of wnt8a caused activation of β-catenin and which leads to increased proliferative response and accelerated regeneration [100]. We have validated array expression data of wnt8a and β-catenin by qRT-PCR and in situ hybridization experiments (Fig. 8). The qRT-PCR analysis shows both wnt8a and β-catenin are upregulated in 7 dpi cord, where we see highest rate of proliferation (Fig. 8A–B). In situ hybridization analysis with wnt8a transcript showed high expression in 7 dpi cord and absent in uninjured and 3 dpi cord. Many cells with wnt8a mRNA transcript are present around the central canal and are ependymal cells (Fig. 8F) which are thought to be precursors [12]. These data confirm our previous observation where we have showed the presence of proliferating BrdU⁺ precursors in ependyma [6]. Few wnt8a expressing cells are also present in white matter of 7 dpi cord (Fig. 8F and F.1). The spatio-temporal distribution of wnt8a transcript in regenerating cord suggests that it is upregulated predominantly in ependymal cells, which actively proliferates following injury. Another member of Wnt family wnt7a is upregulated in wound epidermis, cartilage, spinal cord and muscle in regenerating tail in urodele [102]. Both wnt7a and wnt4b are highly upregulated in regenerating cord in early time points and validated by qRT-PCR (Fig. 8C–D. Apart from the positive regulators there are many negative regulators in Wnt signaling pathway like dixdc1, gsk3, sfrp1, sfrp1b and sfrp5 suggesting both presence of pro and negative Wnt signaling during spinal cord regeneration.

BMP/TGF beta signaling.

Members of BMP/TGF beta signaling pathways are differentially up or down regulated during the different events of regeneration (Figure S8A and Table S15). This list contains both TGFβ signaling genes and their target genes. We observe different group of genes based on temporal expression pattern. Among members of TGFβ signaling family, temporal expression pattern of tgfβ1 gene is distinctly different being upregulated immediately after injury (1 dpi cord) and down regulated in subsequent time periods. This gene is associated with inflammatory process early during regeneration phase as discussed previously in section 3A. Members of BMPs and growth and differentiation factor 11 (Gdf11) are known to control neurogenesis in olfactory neuroepithelia and help in maintenance of progenitor population by replacing dying neurons [103], [104]. We have also found bmp5 and gdf11 both are upregulated in 3 dpi and 10 dpi cord and these genes are known to play key role in this pathway.

JAK-STAT signaling.

Three members of JAK-STAT pathway like stat1b, stat5.2 and stat3 are upregulated during regeneration of zebrafish spinal cord (Figure S8B and Table S15). Stat5.2 and stat1b are both upregulated at 3 dpi and 10 dpi cord and down regulated in uninjured control and other injury time points. Stat3 showed different temporal expression pattern than others and uprgulated in 1 dpi and 3 dpi cord. Similarly, stat3 interacting genes like rtk6 and target genes like socs3a and socs3b, mapk2l and ptpn6 all are upregulated in 1 dpi and 3 dpi cords. Stat3 expression has been validated and described previously (Fig. 3), which is an important regulator of inflammation and may promote cell survival and proliferation. The temporal expression of stat3 and its downstream targets like socs3a and socs3b are similar and may indicate the immunomodulatory role of socs3 by regulating JAK-STAT signaling as suggested earlier [105].