Identification of key genes involved in recovery from spinal cord injury in adult zebrafish.

Zebrafish are an effective vertebrate model to study the mechanisms underlying recovery after spinal cord injury. The subacute phase after spinal cord injury is critical to the recovery of neurological function, which involves tissue bridging and axon regeneration. In this study, we found that zebrafish spontaneously recovered 44% of their swimming ability within the subacute phase (2 weeks) after spinal cord injury. During this period, we identified 7762 differentially expressed genes in spinal cord tissue: 2950 were up-regulated and 4812 were down-regulated. These differentially expressed genes were primarily concentrated in the biological processes of the respiratory chain, axon regeneration, and cell-component morphogenesis. The genes were also mostly involved in the regulation of metabolic pathways, the cell cycle, and gene-regulation pathways. We verified the gene expression of two differentially expressed genes, clasp2 up-regulation and h1m down-regulation, in zebrafish spinal cord tissue in vitro. Pathway enrichment analysis revealed that up-regulated clasp2 functions similarly to microtubule-associated protein, which is responsible for axon extension regulated by microtubules. Down-regulated h1m controls endogenous stem cell differentiation after spinal cord injury. This study provides new candidate genes, clasp2 and h1m, as potential therapeutic intervention targets for spinal cord injury repair by neuroregeneration. All experimental procedures and protocols were approved by the Animal Ethics Committee of Tianjin Institute of Medical & Pharmaceutical Sciences (approval No. IMPS-EAEP-Q-2019-02) on September 24, 2019.

Chinese Library Classification No. R459.9; R363; R364

Introduction

Spinal cord injury (SCI) causes serious complications and permanent disability, thus placing a heavy burden on society and families (Bai et al., 2014; Kang et al., 2018; Kim et al., 2018; Huo et al., 2021). Effective repair of nerve injury has always posed difficulty for neuroscience researchers (Hutson and Di Giovanni, 2019; Wang et al., 2019; Yang et al., 2019). In stark contrast to mammals, adult zebrafish can recover to pre-injury levels after complete spinal cord transection, showing their superior spinal cord regeneration capacity (Pan et al., 2013; El-Daher and Becker, 2020). Therefore, zebrafish are an ideal model to investigate the mechanisms underlying spinal cord regeneration (Cigliola et al., 2020).

For the past 20 years, zebrafish spinal cord regeneration studies have shed light on mammalian SCI regeneration. For example, semaphorin 4D contributed to the successful regeneration of adult zebrafish axons and the recovery of swimming ability after SCI (Peng et al., 2017; Li et al., 2020). Similar successes involve L1 (Becker et al., 2004), L1.2 (Chen et al., 2016), major vault protein (Pan et al., 2013), E2f4 (a member of the DREAM complex) (Sasagawa et al., 2016), recombinant tenascin c (Yu et al., 2011a), syntenin-a (Yu and Schachner, 2013), activating transcription factor 6 (Ji et al., 2020), and miR-133b (Yu et al., 2011b; Danilov et al., 2020). In addition, Liu et al. (2016) suggested that Mcam promotes neurogenesis and angiogenesis to repair SCI. In another study, Fang et al. (2014) reported that Hmgb1 promotes recovery from SCI by stimulating neural regeneration and angiogenesis. Guo et al. (2011) and Ogai et al. (2012) showed that Sox11b affects neural regeneration and the differentiation fate of endogenous neural stem cells after SCI. Mokalled et al. (2016) showed that Ctgfa is involved in the formation of bridge-like structures in spinal cord regeneration. Glial bridges are also required for natural spinal cord regeneration in rats (Tran et al., 2018). These results indicate that zebrafish and mammals have both common and unique mechanisms related to spinal cord regeneration. Thus, studying the mechanisms behind zebrafish SCI regeneration provides new research targets for mammalian SCI regeneration.

The subacute phase, which is approximately 2 weeks post-injury, is a key time point for adult zebrafish spinal cord regeneration, in which the key processes, such as tissue bridging and proliferating, are involved (Cigliola et al., 2020). According to Reimer et al. (2008), newly proliferated motor neurons migrate to the damaged area to replace lost neurons. The peak for motor neuron proliferation, induced by injury, occurs on the 14th day after SCI (Reimer et al., 2008). Mokalled et al. (2016) conducted RNA sequencing of spinal cord cells 2 weeks after SCI in zebrafish. Therefore, we used this dataset to analyze differences in gene expression between the control and injury groups. We also performed real-time quantitative polymerase chain reaction (RT-qPCR) verification on the key genes that were screened out. The purpose of this study was to find key genes and pathways for axon regeneration after SCI in zebrafish and to provide new ideas and targets for the treatment of SCI in mammals, using bioinformatics analysis.

Materials and Methods

Animals

We obtained one hundred 6-month-old adult female zebrafish from Shanghai FishBio Co., Ltd, Shanghai, China. Zebrafish were kept in 3.5-liter tanks on a 14-hour light and 10-hour dark cycle at 28°C and fed twice a day. All zebrafish experiments complied with the standards of the Animal Ethics Committee of Tianjin Institute of Medical & Pharmaceutical Sciences (approval No. IMPS-EAEP-Q-2019-02) on September 24, 2019.

SCI modeling

One hundred zebrafish with body lengths of 2.5–3.0 cm were divided into two equal groups using a random-number table and a double-blind procedure (n = 50 for each group). The zebrafish were subjected to spinal cord transection or sham surgery, as previously described (Fang et al., 2012; Peng et al., 2017). Briefly, zebrafish were anesthetized by immersion in 0.033% ethyl 3-aminobenzoate methanesulfonate (MS-222, MilliporeSigma, St. Louis, MO, USA) in phosphate-buffered saline (PBS) for 3–5 minutes. A longitudinal incision was made on the zebrafish’s side with fine scissors (15000-00, Fine Science Tools GmbH, Heidelberg, Germany) to expose the spine. For the SCI group, the spinal cord was cut with scissors (15000-03, Fine Science Tools, Inc.) between two vertebrae about 4 mm caudal to the brainstem/spinal cord transitional junction. The SCI fish had a straight body posture and were not able to swim freely after surgery. For the sham group, an incision was made to expose the spinal cord, but no SCI was performed.

Swim tracking

Swim tracking analysis was performed on zebrafish at 3, 7, and 14 days after SCI in a 22 cm × 10 cm × 15 cm tank containing 5 cm-deep water. Each fish’s swim path was recorded with a digital camera (Ordro HDV-V7, Shenzhen Boya Times Technology Co., Ltd., Shenzhen, China) for 5 minutes, and the total swimming pattern was quantified using the Animal Tracker program of ImageJ software (Schneider et al., 2012).

RNA sequencing dataset

We used an online dataset of RNA sequencing results from spinal cord tissue of zebrafish 2 weeks after SCI to analyze changes in subacute-phase locomotion and to perform RT-qPCR on genes of interest. The original GSE77025 dataset was downloaded from the NCBI Gene Expression Omnibus (http://www.ncbi.nlm.nih.gov/geo/). The dataset was submitted by Mokalled et al. (2016) and is based on the GPL14875 platform (Illumina HiSeq 2000 [Danio rerio]). The GSE77025 dataset includes two sham-injured spinal cord samples and two transected SCI samples (Mokalled et al., 2016). According to a study by Mokalled et al. (2016), adult zebrafish underwent complete spinal cord transection, and the spinal cord tissue was collected 2 weeks post-injury. Total RNA from the spinal cord tissue of the zebrafish was extracted 2 weeks later for transcriptome sequencing. The sham-operated group was used as the control.

Quality assessment of RNA sequencing data

We used STAR (Dobin et al., 2013) version 2.2.1 to align reads and quantify transcripts. We then used quantile normalization of transcripts to reduce batch effects and systematic biases from sequencing. Differentially expressed gene (DEG) analysis was performed by DESeq2 (Love et al., 2014) version 1.26.0, followed by principal component analysis (PCA) and Pearson correlation analysis.

Identification of differentially expressed genes

DESeq2 version 1.26.0 was used to identify DEGs (|log2 FoldChange| > 0.5 and false discovery rate [FDR] < 0.001). The volcano plot of the DEGs was drawn using the volcano plot package (Chan, 2018). Then, the DEGs were sorted from high-to-low expression, and the top 410 higher expressed genes (DHEGs) were obtained that play the most important role after SCI. Finally, we determined the top 50 genes in the DHEGs and drawn the heat map (Chan, 2018).

Function and pathway enrichment analysis of differentially expressed genes

Gene ontology (GO) analysis is a common method of gene-function annotation that can be used to analyze transcriptome data (Ashburner et al., 2000). Kyoto Encyclopedia of Genes and Genomes (KEGG) is a database resource for studying advanced functions and biological systems at the molecular level, and it is especially useful for genomic sequencing data (Ogata et al., 1999). DEGs were characterized by GO and KEGG enrichment using Metascape (http://metascape.org) (Zhou et al., 2019).

Protein-protein interaction network and transcription regulation

To further understand the specific molecular mechanisms by which genes regulate neural-repair-related pathways after SCI, we constructed a protein-protein interaction (PPI) network using the STRING database (Szklarczyk et al., 2017). A combined score of > 400 was set. The darker the color of the edge connecting two nodes, the greater the combined score. Cytoscape software (Shannon et al., 2003) version 3.7.1 was employed to analyze the network. The NetworkAnalyzer tool in Cytoscape was used to compute node degree value, which was mapped to node height and width. Therefore, larger nodes resulted in greater node degree values. The fill color and border paint of the nodes represent log2 FoldChange and gene expression after SCI, respectively.

RT-qPCR analysis

At 2 weeks after SCI or sham surgery, total RNA was extracted from the 2-mm segment of spinal cord directly caudal to the lesion site using an Eastep Total RNA extraction kit (Promega, Madison, WI, USA) in accordance with the manufacturer’s protocol. Then, 500 ng of total RNA was reverse transcribed using a TransScript First-Strand cDNA synthesis kit (TransGen Biotech, Beijing, China). RT-qPCR was subsequently performed with TransStart Green qPCR SuperMix (TransGen Biotech), and products were detected with a LightCycler 96 (Roche, Basel, Switzerland). The primers used are shown in . Relative expression of gene transcripts was assessed using the 2–ΔΔCt method, for which ΔΔCt = (Cttarget gene – CtGAPDH) injury – (Cttarget gene – CtGAPDH) sham. Three independent experiments were conducted for each condition, with 10 zebrafish in each group.

Statistical analysis

Swim tracking and RT-qPCR data were expressed as the mean ± SEM. GraphPad Prism version 8.02 software (GraphPad Software, San Diego, CA, USA) was used for statistical analysis. The Student’s t-test was used to compare two independent groups. P < 0.05 was considered statistically significant.

Results

Zebrafish partially recover locomotor ability 2 weeks post-spinal cord injury

The subacute phase is critical to zebrafish spinal cord recovery, when glial bridge formation and spinal cord regeneration occur (Mokalled et al., 2016). Therefore, the 2-week time point was used to investigate the mechanisms controlling locomotor recovery. The study design is presented in . First, a zebrafish SCI model was produced, and sham surgery was used as the control (). Zebrafish with “C” and “S” shapes were considered to have damaged spines and were excluded on day 3 after SCI to ensure the consistency of SCI surgery results. Motor function of zebrafish was observed on days 3, 7, and 14 after SCI (Figure and ). The results show that zebrafish did not recover at 3 days or 7 days after SCI. They recovered 44% of their swimming ability 14 days after SCI (P < 0.001).

Establishment of zebrafish SCI model and evaluation of motor function after SCI.

(A) Experimental design scheme. Zebrafish SCI surgery was performed on day 0 in a dark environment. On day 3, the zebrafish were returned to normal feeding and light/dark cycle conditions, and abnormal fish in the shape of “C” or “S” were removed. Zebrafish motor ability was measured on days 3, 7, and 14 after SCI. The RNA-sequencing datasets relating to samples from zebrafish on day 14 after SCI were downloaded for bioinformatics analysis, and the RT-qPCR validation was performed. (B) Zebrafish model of SCI. A longitudinal incision was made with scissors on the side of the zebrafish to expose the vertebral column. We cut the spinal cord about 4 mm caudal to the brainstem/spinal cord transitional junction between two vertebrae. A sham control lesion was performed by making the incision to expose the spinal cord. (C, D) Swim tracking. The movement of zebrafish was tracked and recorded with a camera for 5 minutes in a quiet environment. The total swimming distance of zebrafish in the tank in 5 minutes was calculated using the ImageJ’s Animal Tracker program (Schneider et al., 2012). ***P < 0.001 (Student’s t-test). Data were expressed as the mean ± SEM (n = 10). RT-qPCR: Real-time quantitative polymerase chain reaction; SCI: spinal cord injury.

RNA sequencing data information and quality

The RNA sequencing dataset GSE77025 (Mokalled et al., 2016) was used to investigate gene-expression changes in the adult zebrafish spinal cord 2 weeks after complete transection. The dataset included two sham-injured spinal cord samples and two SCI samples. PCA analysis showed that the two groups are distributed in different regions, indicating an obvious difference in transcriptomes between the control and injury groups (data not shown). Pearson correlation analysis was used for data quality control. The transcriptomes of the two groups (SCI and sham) are similar when the correlation coefficient is close to 1. Pearson correlation analysis showed that there are significant differences in the transcriptomes of the control and injury groups (). Our results show that the data have good repeatability, are reliable, and can be used for subsequent bioinformatics analysis.

To target the genes that play the most important role after SCI, we first obtained the DEGs after SCI by using DESeq2 (|log2 FoldChange| > 0.5 and FDR < 0.001). In total, 7762 DEGs were identified after SCI, including 2950 up-regulated and 4812 down-regulated DEGs. We sorted the DEGs from highest-to-lowest-expressed genes in the injury group and obtained the top 410 higher expressed genes (DHEGs). Finally, we obtained the top 50 up-regulated and down-regulated DHEGs (), which are also strongly expressed. We predicted these genes would play important roles in SCI repair.

Differentially expressed genes 2 weeks after SCI in zebrafish.

(A) Volcano plot of DEGs between control and injury groups. We screened DEGs based on FDR < 0.001 and |log2FoldChange| > 0.5. The green, red and black points represent down-regulated genes, up-regulated genes, and no significant difference, respectively. (B) Heat map of the top 50 up-regulated and down-regulated DHEGs. The hierarchical clustering heat map of DHEGs was screened based on |log2FoldChange| > 0.5 and FDR < 0.001. The top 50 genes were strongly expressed and were also DHEGs in the injury group. Green indicates down-regulation and red indicates up-regulation. The two columns on the left are the control samples. The two columns on the right are the injury samples. DEG: Differentially expressed gene; DHEG: differentially higher expressed gene; FDR: false discovery rate; SCI: spinal cord injury.

GO and KEGG pathway enrichment of differentially expressed genes

DEGs were characterized by GO and KEGG enrichment using Metascape. The top 10 pathways that reached a statistically significant value (P < 0.05) are shown in and Tables and . For the GO analysis, the up-regulated genes are mainly enriched in respiratory chain, synapse, neuron part, cellular component morphogenesis, and Ca2+-binding pathways. The up-regulated DEGs relate to neuronal repair and regeneration. The down-regulated genes are mainly enriched in RNA processing, DNA metabolic process, chromosome organization, and cell cycle. For the KEGG pathway analysis, the signaling pathways of up-regulated genes are mainly enriched in oxidative phosphorylation, NADH dehydrogenase (ubiquinone) pathway, and endocytosis. The signaling pathways of down-regulated differential genes are mainly enriched in RNA transport, protein processing in endoplasmic reticulum, and cell cycle. The GO and KEGG analyses of down-regulated DEGs show that these genes are mainly involved in cellular replication.

GO and KEGG pathway enrichment of differentially expressed genes.

(A) Heat map of top 10 GO terms of DEGs up-regulated after injury. (B) Heat map of top 10 KEGG pathways of DEGs up-regulated after injury. (C) Heat map of top 10 GO terms of DEGs down-regulated after injury. (D) Heat map of top 10 KEGG pathways of DEGs down-regulated after injury. The vertical and horizontal axes represent the names of pathways and enrichment significance, respectively. The color gradient from green to red indicates the gene ratio from low to high. DEG: Differentially expressed gene; GO: gene ontology; KEGG: Kyoto Encyclopedia of Genes and Genomes.

Distribution of the top 50 genes into the top 10 GO pathways

Through GO analysis, we found that the up-regulated genes are enriched in synapses and neurons. Hence, we assumed that up-regulated genes are beneficial for SCI repair. We distributed the top 50 genes into the top 10 GO pathways to determine how the top 50 up-regulated genes perform certain functions and to identify the most important genes after SCI (). Most of the top 30 of the top 50 genes fall into at least one pathway. Furthermore, calb2b, sptan1, map4l, stmn2a, clasp2, and ache fall into at least three pathways, which implies that these genes may play key roles in their respective pathways. These genes are related to the components and repair of neurons, such as respiratory chain, synapse, neuron part, and Ca2+ binding. We also distributed the top 50 down-regulated genes into the top 10 GO pathways. The results showed that h1m and ddx3xb fall into four pathways , and all the genes are involved in cellular replication.

Distribution of the genes into the GO pathways.

(A) The distribution of the top fifty up-regulated genes into the top ten GO terms. (B) The distribution of the top fifty down-regulated genes into the top ten GO terms. Blue indicates that the gene is enriched in one key pathway. The pink-to-red gradient indicates that the gene is enriched in two or more key pathways. GO: Gene ontology.

Key genes were verified by RT-qPCR

Key genes in and are involved in multiple pathways, suggesting that they play an important role in spinal cord repair. Therefore, we performed RT-qPCR to verify the expression changes of calb2b, sptan1, map41, stmn2a, clasp2, ache, h1m, and ddx3xb after SCI. The mRNAs encoding calb2b, sptan1, map41, stmn2a, clasp2, and ache are up-regulated, and the mRNAs encoding h1m and ddx3xb are down-regulated in the RNA sequencing dataset for SCI. According to Mokalled et al. (2016), ctgfa is up-regulated 2 weeks after zebrafish spinal cord crush injury (P < 0.05). Ctgfa was used as a positive control for this dataset, and it was also up-regulated consistently, as previously reported (Mokalled et al., 2016) (). RT-qPCR results verified that clasp2 is up-regulated (1.46 times, P < 0.05) and h1m is down-regulated (0.46 times, P < 0.001) after SCI (). The expression changes are consistent with the RNA sequencing (RNAseq) transcriptome data. The RT-qPCR results show that sptan1 and stmn2a were unchanged and calb2b, map4l, and ache were downregulated to 0.48, 0.33, and 0.62 times the sham group, respectively, and ddx3xb was up-regulated to 1.69 times the sham group ().

RT-qPCR verified partial DEGs and PPI.

(A) RT-qPCR of the spinal cord tissue was conducted 14 days post-injury. Relative expression was normalized to gapdh levels and to the average expression of uninjured controls. The red border indicates up-regulation in RNA sequencing, and the green border indicates down-regulation in RNA sequencing. Three independent experiments were conducted for each condition, with 10 zebrafish in each group. *P < 0.05, ***P < 0.001 (Student’s t-test). Data were expressed as the mean ± SEM (n = 10). (B) PPI network based on Clasp2. (C) PPI network based on H1m. Green and yellow colors indicate proteins with up-regulated and down-regulated expression, respectively. The node fill color shows the changes in gene expression. The more green the color, the greater the up-regulation post-injury; the more yellow the color, the greater the down-regulation post-injury. The node border paint shows the fragments per kilobase million mapped reads of the gene post-injury. The darker the red, the higher the expression post-injury. The degree of network connectivity is mapped to the size of the node shape. The combined score of the PPI network is mapped to the color of the network connection line. The darker the color, the greater the possibility of interaction. DEG: Differentially expressed gene; PPI: protein-protein interaction; RT-qPCR: real-time quantitative polymerase chain reaction; SCI: spinal cord injury.

Protein-protein interaction networks of Clasp2 and H1m

To further predict the specific molecular mechanisms of Clasp2 and H1m for regulating nerve injury repair-related pathways after SCI, protein-protein interaction networks were created for Clasp2 protein or H1m protein using the STRING database version 11.0 (https://string-db.org/). PPI network analysis showed the interactions among the DEGs. There were nine DEGs interacting with Clasp2 (). Pafah1b1a (Drerup et al., 2010) is involved in the development of the nervous system, and Ckap5 is associated with microtubule dynamics (Cassimeris et al., 2009; Chanez et al., 2015). There were 18 DEGs interacting with H1m (); most of their functions are related to cell cycle and embryonic development, such as Cdc6 (Yao et al., 2017), Ccna1 (Beaudoin et al., 2018), and Nanog (He et al., 2020; Palfy et al., 2020).

Discussion

Adult zebrafish spinal cord regeneration is multi-faceted. We chose to identify key genes by using an established bioinformatics database. Mokalled et al. (2016) generated RNAseq data from zebrafish with spinal cord injury and found ctgfa to be the key to zebrafish spinal cord regeneration by promoting formation of glial bridge-like structures. By further analysis of the same dataset, we identified that clasp2 and h1m may also play a key role in spinal cord regeneration.

We screened 2950 up-regulated genes and 4812 down-regulated genes by genetic differential analysis. Similar to previous studies (Fu et al., 2019; Li et al., 2020), we found that some up-regulated genes were also highly expressed in regenerated neurons 2 weeks after SCI. In the GO analysis, the up-regulated genes were mainly enriched in respiratory chain, synapse, neuron part, cellular component morphogenesis, and Ca2+ binding. The respiratory chain is in the inner membrane of mitochondria and undergoes oxidative phosphorylation via the electron transport chain to supply ATP required for metabolism (Mazat et al., 2020). Mitochondria are damaged after SCI, resulting in electron transport chain dysfunction, mitochondrial-membrane-potential disappearance, and ATP synthesis obstruction (DiMauro and Schon, 2003; Li et al., 2019). There is evidence that axon regeneration requires the consumption of 50% of the intracellular ATP (Smith and Gallo, 2018), and hence it relies heavily on the complete respiratory chain to provide energy for axon regeneration (Weng et al., 2017). In this study, proteins expressed by the up-regulated genes participate in the respiratory chain pathway by contributing to the synthesis of ATP and providing energy for axon regeneration. Ca2+ plays a significant role in signal transduction, but its overload activates a variety of protein lysosomes and lipid lysosomes, leading to plasma membrane damage, decreased ATP production, and ultimately skeletal proteolysis (Hall and Springer, 2004; Bandura and Feng, 2019). Ca2+ is also involved in nerve conduction between synapses, and Ca2+ imbalance hinders neurological recovery after SCI (Pivovarova and Andrews, 2010). The proteins expressed by the up-regulated genes in this study participate in the Ca2+-binding pathway, which helps in eliminating the neuron damage caused by Ca2+ overload (Hall and Springer, 2004; Sullivan et al., 2004). Cell proliferation also plays an important role in secondary injury after SCI (Ren et al., 2018). During the acute stage of injury, the glial scar (caused by the proliferation of astrocytes) and the inflammatory response (caused by microglia cells) can limit the scope of injury and remove toxic metabolites over time (Sofroniew, 2018; Tran et al., 2018; Li et al., 2019). However, the inflammatory role maintained by microglia during late axonal growth is detrimental (Totoiu and Keirstead, 2005; Witcher et al., 2015; Tsarouchas et al., 2018). Chondroitin sulfate proteoglycans in the glial scar inhibit axon outgrowth through the protein tyrosine phosphatase sigma (PTPσ) (O’Shea et al., 2017; Tran et al., 2018). Therefore, it is important to regulate the proliferation of neural cells according to the spatiotemporal changes of the microenvironment after injury (Fouad et al., 2005; Duncan et al., 2018). The microenvironment of SCI regions can be regulated by target genes or pathways to promote SCI repair (Fan et al., 2018). In addition, the regulation of these target genes or pathways has been applied to transplanted exogenous stem cells, so that the cells have good nerve-regeneration ability (Lu et al., 2017; Duncan et al., 2018).

In the KEGG pathway analysis, the signaling pathways of DEGs were mainly enriched in oxidative phosphorylation, NADH dehydrogenase (ubiquinone) pathway, RNA transport, protein processing in endoplasmic reticulum, and cell cycle. These pathways focus on the cell cycle and oxidative stress. The results of GO enrichment analysis were also similar. Results from the GO and KEGG analyses suggest that these differentially expressed genes promote regeneration in two ways. On one hand, these genes promote cell metabolism, such as oxidative phosphorylation and NADH dehydrogenase (ubiquinone) pathway, while on the other hand, the genes influence the fate of cells, such as by controlling RNA transport, DNA degradation, and protein processing. Acute inflammation is intertwined with axon regeneration, as described previously in zebrafish (Bollaerts et al., 2017). The acute inflammatory phase is beneficial to central nervous system regeneration in zebrafish. The anti-inflammatory drug dexamethasone significantly inhibits the regeneration of motor neurons in zebrafish (Ohnmacht et al., 2016). However, the inflammatory response is suppressed in the first week post-injury without any chronic consequences in zebrafish (Ceci et al., 2018). Interestingly, in this RNAseq data, DEGs were not significantly enriched in a specific inflammatory pathway. This could be due to the time period used (14 days post-injury), which was beyond the acute phase of inflammation. Because RNAseq data were obtained in the whole spinal cord tissue homogenate, we could not figure out the distinctive expression pattern of inflammatory cells such as microglia or macrophages. Therefore, single cell sequencing may be used in the future to better understand the function of inflammation in specific neuronal or immune cell types.

We distributed the top 50 up-regulated genes into the top 10 GO pathways. The results showed that calb2b, sptan1, map4l, stmn2a, clasp2, and ache fall into at least three pathways. Moreover, the top 50 down-regulated genes, including h1m and ddx3xb, fall into at least four pathways. The RT-qPCR results show that the expression of clasp2 and h1m is consistent with the sequencing results. These genes may play key roles in their respective pathways.

Microtubules play an important role in the formation and maintenance of axons and dendrites in neurons. Clasp2 is a microtubule plus-end tracking protein that promotes the stabilization of dynamic microtubules (Maki et al., 2015). Because it is induced in zebrafish upon SCI, we hypothesized that it has a pro-regenerative role. Clasp2 mediates growth-cone microtubules and is in exquisite balance with the co-regulation of GSK3 (Hur et al., 2011). It has also been reported that Clasp2 regulates neuronal polarity and synaptic function, and it has been found that overexpression of Clasp2 can lead to an increase in presynaptic-terminal circumference and total synapse number (Beffert et al., 2012). Dillon et al. (2017) demonstrated that Clasp2 regulates neuronal production and controls neuronal migration in mice. Through PPI network analysis, we found that Clasp2 interacts with microtubule-related proteins. Ckap5 plays an important role in axon growth and regulates the growth cone by promoting interaction between microtubules and F-actin (Lowery et al., 2013; Slater et al., 2019). Pafah1b1a (also known as lissencephaly-1 or LIS1) activates dynein, which facilitates microtubule transportation that is critical for axon regeneration (Elshenawy et al., 2020; Htet et al., 2020; Marzo et al., 2020; McKenney, 2020). We hypothesized that Clasp2 may promote the growth of neuronal axons and growth cones (). Neurogenesis is a key process in zebrafish spinal cord regeneration (Ceci et al., 2018). Neurogenesis may come from transdifferentiation as well as from direct differentiation from neural stem cells or progenitor cells. Notwithstanding transcription factors, non-coding RNAs, such as microRNAs, play an emerging role in transdifferentiation (Sabater et al., 2020). Therefore, further mining of the RNAseq data for the microRNAome or investigating single-cell RNAseq techniques are needed to study neurogenesis.

Schematic map of the involvement of key genes in the two facets of spinal cord regeneration.

The first facet of spinal cord regeneration is axon regeneration promoted by Clasp2 and its related proteins, Ckap5 and Pafah1b1a. The second facet of spinal cord regeneration is endogenous neural stem cell differentiation, regulated by H1m and its related proteins, Nanog and Ubb. The up-regulated genes are orange in color, and the down-regulated genes are blue in color.

In this study, h1m is down-regulated after SCI in zebrafish. We hypothesized that H1m may block spinal cord regeneration. Hence after injury, reduced H1m expression may trigger the spontaneous regeneration programmed in adult zebrafish. The exact mechanism needs further elucidation. H1m is a linker histone that is in the nucleus of primordial germ cells and plays an important role in zebrafish embryonic growth. When the zebrafish embryo develops to the gastrulation stage, the expression of h1m in primordial germ cells and somatic cells declined rapidly (Muller et al., 2002). In Xenopus, Nap-1 was shown to be the main chaperone of H1m. It mediated chromosome aggregation during mitosis (Muller et al., 2002). These studies show that H1m may modulate regeneration at the chromatin level. PPI network analysis has also revealed that H1m interacts with ubiquitin b (Ubb) and several cell cycle-related proteins. The studies found that knocking down UBB in mice affected the differentiation of neural stem cells into neurons (Jung et al., 2018; Lim et al., 2019). Therefore, in this study, upregulation of Ubb may contribute to the differentiation of neural stem cells into neurons by promoting ubiquitin-related degradation of H1m. There are many down-regulated proteins that interact with H1m, such as Cdc6, Ccna1, and Nanog. Cdc6 is crucial for the initiation of DNA replication and regulation of the early stage of DNA replication (Yao et al., 2017). Ccna1 is expressed in the testis, brain, and leukemia cell lines, and is thought to play an important role in controlling the meiotic cell cycle in the germ line (Beaudoin et al., 2018). Nanog plays an important role in embryonic development by blocking embryonic stem cell differentiation (He et al., 2020; Palfy et al., 2020). Downregulation of Nanog after SCI may promote the differentiation of stem cells, such as the differentiation of neural stem cells into neurons. Therefore, down-regulation of H1m after SCI promotes differentiation of endogenous neural stem cells into neurons by altering the rate of interactions with these proteins ().

There are limitations to our study. The dataset we used contains only two sequenced samples for each group, so the evidence may not be strong enough to draw meaningful conclusions. However, we verified the candidate genes, identified using the dataset, by RT-qPCR, which ensures the credibility of the experimental results. In the RT-qPCR experiment, we performed three independent replicates with 10 zebrafish in each group. The mechanism for Clasp2 promotion of axon regeneration after SCI has not been verified, and different neuronal cell types may have different regulation patterns (Hur et al., 2011). There is evidence that Clasp2 plays a role in promoting axon regeneration in vivo, but it remains unclear whether Clasp2 promotes motor-function recovery from SCI. We need to discover the temporal and spatial expression patterns of these genes as well as the specific cell-types that express these genes. With our established swim assay, the effects of key genes’ knockdown and overexpression on swim ability should be conducted to investigate the role of these genes in spinal cord regeneration. Furthermore, we need to clarify whether these genes are the key to boosting regeneration in mammalian species by first studying these genes in rodents or even non-human primates.

Zebrafish can recover 44% of their motor ability within 2 weeks after SCI. We analyzed RNA sequencing data from the GEO database for zebrafish 2 weeks after SCI. Through genetic differential analysis, GO, and KEGG, we identified the key genes for acute-phase SCI recovery. RT-qPCR confirmed that clasp2 was up-regulated after SCI and h1m was down-regulated after SCI, both of which are involved in multiple signaling pathways. Up-regulated Clasp2 is a microtubule-binding protein that is associated with Ckap5 and is responsible for axon elongation modulated by microtubules. Pafah1b1a is a dynein-associated protein responsible for microtubule transportation. H1m down-regulation may be attributed to modulation of endogenous stem cell differentiation upon SCI. Taken together, findings from this study reveal that up-regulation of Clasp2 promoted axonal regeneration and down-regulation of H1m promoted endogenous neural stem cell differentiation. These changes may exist as key contributors to zebrafish spinal recovery in the acute phase (). Clasp2 and H1m are potential therapeutic targets for spinal cord regeneration. Our study provides new candidate genes for the repair of SCI.

Additional files:

: Pearson correlation of samples.

Pearson correlation of samples.

Pearson correlation analysis showed the correlation coefficient in the data quality control. When the correlation coefficient was close to 1, the two sets of samples were similar in transcriptome. Each coefficient in the correlation matrix was represented as a square, with the color of the square representing the degree of correlation. The darker the blue, the closer the correlation coefficient was to 1. The above results showed that the data had good repeatability and reliable, and could be applied to the subsequent bioinformatics analysis.

RT-qPCR verified partially differentially expressed genes (DEGs)

RT-qPCR verified partially differentially expressed genes (DEGs).

RT-qPCR of the spinal cord tissue was conducted at 14 days post injury. Relative expression was assayed in 3 pools of 10 fish for each and normalized to GAPDH levels and to the averaged expression of uninjured controls. The mRNA relative expression levels of calb2b (A), sptan1 (B), map4l (C), stmn2a (D), ache (E) and ddx3xb (F). The red border indicates up-regulation in RNA sequencing, and the green border indicates down-regulation in RNA sequencing. Three independent experiments were conducted for each condition with 10 zebrafish in each group. ***P < 0.001 (Student's t-test). Data were expressed as the mean ± SEM.

Table 1

Primers for real-time qPCR used in this study

Gene	Sequence (5’–3’)	Product size (bp)
gapdh	Forward: GTG TAG GCG TGG ACT GTG GT	20
	Reverse: TGG GAG TCA ACC AGG ACA AAT A	22
ctgfa	Forward: CTG GTG TAA GCC TAG TTC TGG	21
	Reverse: GGG TCG CAA ACA TCT CGT	18
calb2b	Forward: GGA AAC TGG GAC TGT CTG	18
	Reverse: TTG CCA TCC TTG TCG TAA T	19
sptan1	Forward: AGG AAC AAG CCG ACT ACT G	19
	Reverse: AAG AGC GAG CGG GTG AAT	18
map4l	Forward: ATA CCC AAA CAC TCA CCC	18
	Reverse: CAC AGA CTG GCA TCT CAC	18
stmn2a	Forward: GCT CTG GCT GAG AAA CGT	18
	Reverse: CCT TGT TCT GCT CCA TCT T	19
clasp2	Forward: TGG AGG CAC ATA AAG ACC	18
	Reverse: TGA CTG GAT GAT GGG ACA	18
ache	Forward: TTT CCT TTG TGC CTG TCG	18
	Reverse: CCT CAT CCT GAT TTA CGC	18
h1m	Forward: AGC TCC TGA AGA CTC CAA C	19
	Reverse: TAG CCA CGA ATA GCC TGA	18
ddx3xb	Forward: CAG AAC AGC AGC AGT GGC	18
	Reverse: TTT CCG AAT CCT CCC TTA C	19

Table 2

GO pathways enrichment of DEGs in 2 weeks after injury

GO term	Description	Count	q-value
Upregulated
GO:0070469	Respiratory chain	43	2.26E-16
GO:0045202	Synapse	109	1.69E-15
GO:0097458	Neuron part	146	1.77E-14
GO:0032989	Cellular component morphogenesis	152	6.59E-11
GO:0043235	Receptor complex	69	8.89E-08
GO:0022839	Ion gated channel activity	77	9.56E-07
GO:0008137	NADH dehydrogenase (ubiquinone) activity	14	1.24E-06
GO:0007155	Cell adhesion	109	1.71E-06
GO:0030029	Actin filament	101	2.16E-06
GO:0005509	Calcium ion binding	124	9.1E-06
Downregulated
GO:0006396	RNA processing	238	4.82E-63
GO:0006259	DNA metabolic process	201	2.04E-59
GO:0051276	Chromosome organization	217	1.35E-56
GO:0022402	Cell cycle process	194	5.86E-53
GO:0005694	Chromosome	197	8.15E-49
GO:0030684	Preribosome	56	4.37E-34
GO:0000003	Reproduction	107	2.18E-29
GO:0005654	Nucleoplasm	169	2.04E-27
GO:0140098	Catalytic activity, acting on RNA	109	2.08E-24
GO:0098687	Chromosomal region	61	4.19E-24

There are top 10 significantly enriched gene GO terms of up-regulated and down-regulated DEGs in spinal cord injury. DEG: Differentially expressed genes; GO: gene ontology.

Table 3

KEGG pathways enrichment of DEGs in 2 weeks after injury

KEGG Pathway	Description	Count	q-value
Upregulated
dre00190	Oxidative phosphorylation	57	2.67E-16
M00147	NADH dehydrogenase (ubiquinone) 1 beta subcomplex	10	3.7E-06
M00101	Cholesterol biosynthesis, squalene 2,3	9	1.62E-05
M00146	NADH dehydrogenase (ubiquinone) 1 alpha subcomplex	10	0.00027
M00142	NADH: ubiquinone oxidoreductase, mitochondria	6	0.000849
dre04261	Adrenergic signaling in cardiomyocytes	43	0.002158
dre04260	Cardiac muscle contraction	24	0.01052
dre_M00148	Succinate dehydrogenase (ubiquinone)	4	0.029785
dre04514	Cell adhesion molecules (CAMs)	28	0.048195
dre04144	Endocytosis	58	0.048195
Downregulated
dre03013	RNA transport	71	4E-22
dre03008	Ribosome biogenesis in eukaryotes	47	7.96E-22
dre03030	DNA replication	28	4.76E-16
dre04141	Protein processing in endoplasmic reticulum	64	3.81E-14
dre04110	Cell cycle	54	5.66E-14
dre03460	Fanconi anemia pathway	28	1.72E-11
dre03040	Spliceosome	43	4E-08
dre00240	Pyrimidine metabolism	34	2.71E-06
dre_M00693	Cell cycle	8	1.32E-05
dre03018	RNA degradation	28	1.77E-05

There are KEGG terms of up-regulated and down-regulated DEGs in spinal cord injury. DEG: Differentially expressed genes; KEGG: Kyoto Encyclopedia of Genes and Genomes.