Abstract
Regenerative medicine offers hope for patients with diseases of the central and peripheral nervous system. Urodele amphibians such as axolotl display an exceptional regenerative capacity and are considered as essential preclinical model organisms in neurology and regenerative medicine research. Earlier studies have suggested that the limb regeneration ability of this salamander notably decreases with induction of metamorphosis by thyroid hormones. Metamorphic axolotl requires further validation as a negative control in preclinical regenerative medicine research, not to mention the study of molecular substrates of its regenerative abilities. In this study, we report new observations on the effect of experimentally induced metamorphosis on spinal cord regeneration in axolotl. Surprisingly, we found that metamorphic animals were successful to functionally restore the spinal cord after an experimentally induced injury. To discern the molecular signatures of spinal cord regeneration, we performed transcriptomics analyses at 1- and 7-days postinjury (dpi) for both spinal cord injury (SCI)-induced (experimental) and laminectomy (sham) groups. We observed 119 and 989 differentially expressed genes at 1- and 7-dpi, respectively, while the corresponding mouse orthologous genes were enriched in junction-, immune system-, and extracellular matrix-related pathways. Taken together, our findings challenge the prior notions of limited regenerative ability of metamorphic axolotl which exhibited successful spinal cord regeneration in our experience. Moreover, we report on molecular signatures that can potentially explain the mechanistic substrates of the regenerative capacity of the metamorphic axolotl. To the best of our knowledge, this is the first report on molecular responses to SCI and functional restoration in metamorphic axolotls. These new findings advance our understanding of spinal cord regeneration, and may thus help optimize the future use of axolotl as a preclinical model in regenerative medicine and integrative biology fields.
Introduction
Regenerative medicine offers hope for patients with diseases of the central and peripheral nervous system. Spinal cord injuries are quite common and stand to benefit from regenerative medicine science. The spinal cord accommodates neurons and supporting cells, forming a tubular conduit through which motor and sensory information are relayed between the brain and peripheral tissues, acting as a center for coordinating certain reflexes as well (Kirshblum et al., 2011; Harting et al., 2008; Silva et al., 2014).
Most vertebrates display a limited regeneration capacity restricted to certain tissues, while others, including amphibians, possess an astonishing capacity for tissue, organ, and appendages regeneration. Strikingly, salamanders have retained the functional regeneration capability after spinal cord injury (SCI) along with internal organs and extremities (Diaz Quiroz and Echeverri, 2013; Tazaki et al., 2017).
An aquatic salamander, axolotl (Ambystoma mexicanum), exhibits remarkable abilities to functionally restore the brain (Maden et al., 2013) and spinal cord (Mchedlishvili et al., 2007) following an injury, and therefore has emerged as an important experimental model to identify the molecular substrates of central nervous system (CNS) regeneration, and in preclinical studies in life sciences (Öktem et al., 2019). Axolotl does not undergo metamorphosis naturally, and this interesting characteristic results in lifelong-lasting neoteny that may account for retaining the regenerative capacity in adulthood.
Metamorphosis in axolotl can be induced by administration of thyroid hormones (McCusker and Gardiner, 2011); morphological and structural alterations with metamorphosis were documented previously (Demircan et al., 2016). Structural changes due to metamorphosis are accompanied by a significant reduction in limb regeneration rate (Demircan et al., 2018; Monaghan et al., 2014), which is well aligned with observations in anurans (Lin et al., 2013; Mescher and Neff, 2006).
Facultative metamorphosis in axolotl provides a unique experimental advantage to characterize the molecular basis of regeneration, such as those in limb or spinal cord regeneration. Decrease in regenerative capacity after metamorphosis for a large group of amphibians (Mescher and Neff, 2006), and reduction in capacity and fidelity of limb regeneration for metamorphic axolotl prompted us to test the regenerative competence of metamorphic stage at spinal cord level. In addition, the metamorphic axolotl requires further validation as a negative control in preclinical regenerative medicine research, not to mention the study of molecular substrates of its regenerative abilities.
In this study, we report new observations on the effect of experimentally induced metamorphosis on spinal cord regeneration in axolotl. To the best of our knowledge, this is the first report on molecular responses to SCI and functional restoration in metamorphic axolotls. These new findings contribute to the understanding of the regulation of spinal cord regeneration and may help optimize the use of axolotl as a preclinical model in regenerative medicine in the near future.
Materials and Methods
Animal care, experimental design, and sample collection
The experimental design of the presented study is shown in Figure 1a. Axolotls (6–8 months old, 46 siblings) used in this study were obtained from Istanbul Medipol University (IMU) Medical Research Center. The Institutional Ethics Committee of the IMU approved the following experimental procedures in this study with authorization number 38828770-E.2499.

Generation of sham (laminectomy only) and experimental (laminectomy followed by SCI) groups.
Animals were maintained at 18°C in 20% Holtfreter's solution as one animal per aquarium and were daily fed with staple food (JBL Novo) as described previously (Demircan et al., 2016). Axolotls were metamorphosed by administration of 50 nM T4 hormone prepared by dissolving
Ten metamorphic axolotls were selected for inducing SCI to be longitudinally observed over a timeline. The SCI model was established according to a previous protocol (Diaz Quiroz et al., 2014). In brief, axolotls first underwent laminectomy, whereby the spine corresponding to the position just above the hindlimb of animals was cracked and extracted following skin incision and muscle removal. Thereafter, SCI was induced by removing the spinal cord in between the two neighboring spines, followed by replacing the skin over the incision site.
Animals were then kept in Holtfreter's solution as described above. Afterward, reflex tests were conducted on these animals at disparate timepoints by which their hindlimbs were squeezed with forceps, and concomitantly, the animals were stimulated to walk to assess their restoration of mobility. Spinal cord regeneration of these animals was inspected through which the skin covering the injury sites was incised, allowing for a macroscopic investigation of the status of the spinal cord at the severed site.
The other 36 metamorphic axolotls were randomly categorized into 2 groups representing the sham (n = 18) and SCI-induced (n = 18) groups, from each of which 2 additional subgroups were formed based on the timepoints of sample collection (first and seventh day, n = 9 for each group). Nine animals were subgrouped into three biological replicates (R1, R2, and R3) to obtain repeatable accuracy, and each replicate was formed by pooling samples from three animals to eliminate individual variations.
The operation on the SCI-induced group and postoperation maintenance was conducted as described above. For sham group, only laminectomy was performed, and animals were placed back to their chambers without an induced SCI. The 1 and 7 days postinjury (dpi) samples were collected from 1 mm rostral and caudal to the injury site in addition to the newly formed tissue. All tissue samples were cryopreserved in liquid nitrogen after the collection and stored at −80°C until RNA isolation.
RNA isolation and next-generation sequencing
Total RNA isolation from the collected samples was performed using TRIzol Reagent (Cat. No. 15596018; Invitrogen™) according to the manufacturer's protocol. Isolated RNAs were quantified by Qubit assay (Cat. No. Q33140; ThermoFisher Scientific), and their quality was assessed using gel electrophoresis (1% agarose). By having sufficient and high-quality RNA samples for sequencing, as a next step, mRNA sequencing libraries were prepared using “TruSeq Stranded mRNA” (Cat. No. 20020594; Illumina) Kit following the producer's instructions. RNA sequencing was performed on Illumina NextSeq 500 platform using “NextSeq 500/550 High Output Kit v2.5 150 cycles” (Cat. No. 20024907; Illumina).
Raw data processing
Illumina BaseSpace platform was used for demultiplexing and removing of sequence adapters of the paired-end sequencing dataset. Merged datasets for each time point were uploaded to SRA database with bioproject ID: PRJNA624092.
Quality control and quality filtering of barcode-free raw sequencing was executed by using FASTQC (Andrews et al., 2015) and Cutadapt (Martin, 2011) softwares, respectively. Only high-quality reads (minimum average Phred quality score of 30) that passed through quality filtering were used for downstream analysis. For the assembly of the raw data, axolotl genome-guided approach was followed. The de novo assembly was generated in a genome-guided manner using Trinity protocol (Grabherr et al., 2011) and axolotl genome (Nowoshilow et al., 2018).
The generated transcriptome contained 156,104 transcripts, 91.3% of which were homologous to at least one transcript found in The National Center for Biotechnology Information (NCBI) nr database. Significantly differentially expressed genes (DEGs; adjusted p < 0.05, |FC| ≥ 2) between SCI-induced and sham groups at 1 and 7 dpi timepoints were identified using RNA-Seq by Expectation Maximization (RSEM) protocol (Li and Dewey, 2011). For further analyses, mouse orthologous of DEGs were found following the elimination of redundant genes, and nonmatching genes to mouse nr database were removed before downstream analyses.
Gene enrichment and pathway analyses
“clusterProfiler” R/Bioconductor package (Huber et al., 2015; Yu et al., 2012) was used first to map the accession numbers of mouse orthologs to EntrezID and subsequently to conduct gene ontology and pathway enrichment analyses using the overrepresentation statistical test. Network plot (cnetplot) from the package was used to visualize the resultant KEGG pathways. KEGG pathway viewer was generated by “pathview” R package (Luo and Brouwer, 2013). All analyses were implemented using R version 3.6.1.
Enrichment parameters were set to p- and q-value = 0.05 cutoffs, “org.Mm.eg.db” (mouse) as the organism database, and Benjamini and Hochberg as the adjusted p-value method. Physical and functional interactions among the DEGs between sham and injury groups at 7 dpi were retrieved from the Search Tool for the Retrieval of Interacting Genes online database (Szklarczyk et al., 2019). Interactions between upregulated and downregulated genes were visualized using Cytoscape (version: 3.7.1) (Shannon, 2003).
Results
Metamorphic animals regenerate the spinal cord after an induced injury
SCI induction on the randomly selected 10 animals revealed their inability to walk or move just upon the injury, indicting the successful implementation of the SCI model on metamorphosed axolotls. Interestingly, a responsive reflex (within 2 weeks) and walking capability (within 2–3 months) was regained within after the injury (Supplementary Movie S1), and this underlines the functional spinal cord restoration. Macroscopic observation of the site of injury consecutively further confirmed the successful regeneration of the spinal cord of the metamorphic axolotls (Fig. 1b).
Filtering steps of acquired DEGs for downstream analyses
A total of 119 genes and 989 genes were found differentially expressed after running RSEM analyses at 1 and 7 dpi, respectively (Supplementary Tables S1 and S2). Afterward, genes which had no matching mouse ortholog were filtered out, as were those having a fold change (FC) ranging between 2 and −2, followed by querying the remaining genes for their corresponding EntrezID to be used in downstream enrichment analysis. This series of filtering left 40 DEGs at 1 dpi and 224 DEGs at 7 dpi. Some of the accession numbers redundantly corresponded to the same gene ID, and so they were accounted for by taking the mean of the corresponding FC, yielding a total of 38 (20 upregulated, 18 downregulated) DEGs at 1 dpi and 219 (116 upregulated, 103 downregulated) DEGs at 7 dpi, which were used for downstream analyses.
Enrichment analysis of DEGs at 1 dpi: small differences between groups
We sought to obtain an overview of some putative enriched biological processes (BPs) by the combined list of up- and downregulated genes between sham and injury groups at 1 dpi, and thereby, overcoming the limitation that would manifest in less informative enrichment results arising from having a relatively small number of genes to begin with (119 DEGs at 1 dpi), many of which were further discarded due to nonmatching mouse orthologs and EntrezID as well as low-fold changes. The results showed favorable enrichment by upregulated genes in processes revealed by KEGG analysis; namely: “Ribosome,” “Apoptosis,” “Tight junction,” “Gap junction,” and “Fructose and mannose metabolism.”
Among the genes enriched in these pathways, two of them (ROCK1 and NFKBIA) are downregulated, and they are represented in “Tight junction” and “Apoptosis” pathways, respectively (Fig. 2a and Supplementary Table S3). Furthermore, the enrichment of gene ontologies by the DEGs at 1 dpi was investigated (Supplementary Table S4). Interestingly, no significant (adjusted p < 0.05) BPs were reported.

Pathway enrichment of DE genes between SCI-induced and sham groups at 1-dpi.
However, “structural constituent of ribosome” and “structural constituent of cytoskeleton” were the top enriched molecular functions (MFs) among a total of 7, and “cytosolic ribosome” and “ribosomal subunit” were found the top enriched cellular components (CCs) among a total of 16. The DEGs in the “Apoptosis” process at 1 dpi were visualized as part of the whole cascade (Fig. 2b) showing the downregulated NFKBIA and the upregulated TUBA1 genes and their interactions with the relevant genes, favoring a state in which prosurvival mechanisms are activated.
Enrichment analysis of DEGs at 7 dpi: marked divergence in molecular signatures
As the seventh dpi timepoint is considered an intrinsic timing for the course of regeneration (Demircan et al., 2017; Monaghan et al., 2006; Voss et al., 2015), and because the number of DEGs (989) is higher than those in 1 dpi, we directly opted to separately enrich the upregulated and downregulated genes between sham and injury groups at this timepoint to get more specialized enrichment by KEGG analysis. Upregulated genes enriched only two processes; “proteasome” and “beta-Alanine metabolism” (Fig. 3a).

Pathway enrichment of DE genes between SCI-induced and sham groups at 7-dpi. A gene-concept network of KEGG pathways enriched by
On the contrary, “IL-17 signaling pathway,” “Hypertrophic cardiomyopathy (HCM),” and “Legionellosis” were the only enriched pathways by the downregulated genes (Fig. 3b). The detailed results are documented in Supplementary Table S5. Gene ontology enrichment was similarly carried out on the separate lists of DEGs at 7 dpi. The top BPs enriched by upregulated genes included “mRNA cis splicing via spliceosome,” “ncRNA metabolic process,” and “ribonucleoprotein complex biogenesis,” among others (Supplementary Table S6). While no significant MFs were reported, “nucleoid,” “mitochondrial nucleoid,” and “peptidase complex” were the top reported significant CCs among others (Supplementary Table S6).
Downregulated genes, on the contrary, did not enrich any CCs, but enriched a single MF (“S100 protein binding”) along with a vast array of BPs, the top of which included “porphyrin-containing compound biosynthetic process,” “cofactor biosynthetic process,” and “tetrapyrrole biosynthetic process” (Supplementary Table S7).
DEG interactions at 7 dpi by STRING network analysis
So far, we have investigated the roles of the DEGs between our compared groups at the two distinct regeneration timepoints by detecting their enrichment in certain ontologies or pathways. We next thought of how the DEGs would interact among themselves at the intrinsic regeneration timepoint (7 dpi), which might give us certain clues as to which gene hubs could be considered crucial for spinal cord regeneration.
The interactions among the up- and downregulated genes were visualized using Cytoscape (Fig. 4) (Shannon, 2003), which yielded a total of 169 nodes (genes) and 424 edges (interactions). The top 5 scores of gene coexpression were the interaction between MCM5 and MCM7, IMPDH2 and FBL, COL1A1 and COL6A2, RPF1 and DCAF13, and FBL and NIP7. Moreover, the downregulated HSP90AA1, upregulated FBL, upregulated CDH1, and upregulated NIP7 were among the genes with the highest degrees. It is intriguing to see that the downregulated proinflammatory cytokine Interleukin 6 (IL-6) had the highest degree (23°) among all. Details on the results of the network can be found in Supplementary Tables S8 and S9.

STRING interactions of DE genes between SCI-induced and sham groups at 7-dpi. STRING interactions among DE genes between both comparisons at 7 dpi visualized by Cytoscape. Darker gray nodes represent upregulated genes while lighter gray nodes represent downregulated genes. White nodes represent genes which were introduced by the network analysis featuring certain interactions.
Discussion
Utilization of salamanders, particularly axolotl, has broadened our understanding of the regeneration concept due to its large regenerative capacity, including the ability to restore the CNS after an injury. As a commonly observed phenomenon among the vast majority of amphibians, based on current literature, regenerative potential in axolotl also significantly decreases following the induced metamorphosis (Demircan et al., 2018; Monaghan et al., 2014). Therefore, metamorphic axolotl has been used as complementary to neotenic counterpart by acting as a regeneration-deficient model (Sibai et al., 2020).
However, based on our observations, the SCI-induced paralyzed metamorphic axolotls were capable of regenerating the missing spinal cord, and hence, regained the lost motor activities. This unexpected finding raises concerns to utilize metamorphic axolotl as a regeneration deficient counterpart of neotenic stage unless the full regenerative capacity at this developmental stage is fully explored.
To understand the molecular mechanisms triggering the plasticity to facilitate the renewal process, gene expression comparison between sham and experimental groups was performed, and DEGs were thoroughly analyzed for 1 and 7 dpi samples. Obtaining a remarkably greater list of DEGs at 7 dpi compared to 1 dpi might pinpoint the diverse molecular profile between sham and experimental groups at 7 dpi. In axolotl, the first 6–72 h after an injury, wound healing program is characterized by infiltration of immune cells to the damaged zone, secretion of proregenerative signals to wound bed, and dedifferentiation of specialized cells into stem/progenitor cells (Godwin and Brockes, 2006; Godwin et al., 2013; Seifert et al., 2012).
Since the followed operative procedure in both sham and experimental groups was identical except for the last step, the generated wound healing-related response to the induced damage in both groups was most probably very much alike and, therefore, a limited number of genes were found as DEGs at 1 dpi. This finding was further supported by getting the limited number of enriched gene ontology (GO) terms and KEGG pathways between sham and experimental groups.
Among the observed KEGG pathways, considering the molecular processes at early phase of regeneration, “Apoptosis” was found noteworthy to have a closer look at. Because TUBA1A, TUBA1B, and TUBA1C genes are highly expressed in nervous system following an injury (Fournier and McKerracher, 1997; Liu et al., 2017; Yi et al., 2018), upregulation of these genes in SCI-induced group may suggest the remodeling of microtubules in damaged neurons to enhance axon regrowth. NF-KB pathway is essential for both activation of innate immunity and cell survival mechanisms under stress conditions (Piva et al., 2006), and in our dataset NFKBIA gene, an inhibitor of NF-KB signaling was found downregulated, which may lead to activation of survival mechanisms in SCI-induced animals.
Seven days postinjury is one of the commonly used timepoints to investigate the molecular composition of blastema, a regeneration specific tissue enriched by dedifferentiated stem/progenitor cells (Demircan et al., 2017; Monaghan et al., 2006; Voss et al., 2015). Upregulated genes in SCI-induced animals were enriched in “proteasome” and “beta-alanine metabolism” pathways.
In accordance with our data, previous studies have shown that proteasome pathway components are upregulated during both mammalian and axolotl regeneration to enhance the survival after an injury along with removal of the damaged proteins (Kitajima et al., 2018; Pasten et al., 2012; Rao et al., 2009). Together with its roles in carnosine biosynthesis (Sale et al., 2010), it has been postulated that beta-alanine functions in elevated levels of the cellular oxygen consumption (Schnuck et al., 2016). Although the putative function of this pathway in regeneration is an open question, it may contribute to regenerative capacity by improving the oxidative metabolism.
HCM was one of the enriched pathways by the downregulated DEGs in SCI-induced animals at 7 dpi. Based on the literature, after an injury or amputation in axolotl, degradation of muscle tissues is detected at around 1–2 dpi (Dwaraka et al., 2019), and acceleration of histolysis of muscle-tissues accompanied by dedifferentiation of specialized muscle cells into stem/progenitor cells is observed at around 7 dpi (Voss et al., 2015). Thus, as presented in our data, enrichment of muscle tissue-related pathways by downregulated DEGs at blastema stage is both expected and essential.
Enrichment of downregulated DEGs (HSP90AA1, USP25, MMP3, IL-6, and PTGS2) in “IL-17 signaling pathway” for injury group indicated a notable decrease in immune system activity at 7 dpi. This is required for successful regeneration as shown in previous studies (Sibai et al., 2020; Tica and Didangelos, 2018; Tsai et al., 2019).
The above notion is further confirmed by enrichment of downregulated IL-12B, SAR1B, TLR5, and IL-6 genes in “Legionellosis pathway” in injury group. Legionella species, whose main targets are epithelial cells and macrophages for intracellular replication, are causative agents for a severe infectious legionellosis disease (Isaac and Isberg, 2014; Liu and Shin, 2019). Activation of this pathway by membrane bounded Toll-like receptors (TLRs), such as TLR5, results in increased proinflammatory response and chemoattractant signal secretion (Liu and Shin, 2019). It has been postulated that, downregulation of proinflammatory pathways following an adequate immune response is imperative to limit the activity of immune system for successful regeneration (Godwin and Brockes, 2006; Mescher and Neff, 2005; Mescher et al., 2017).
Therefore, our findings highlighting the less active immune system in SCI-induced animals than the sham group at 7 dpi by downregulation of “IL-17 signaling pathway” and “Legionellosis pathway” may provide new evidences on the putative link between regeneration and immune system. In a previous study (Kumamaru et al., 2012), it has been demonstrated that lower amount of proinflammatory cytokines such as IL-6 and TNF-α were secreted in young animals compared with adult mice following the SCI. The results from this study underline the increased proinflammatory signals in aged animals that contribute to the amplification of a strong inflammatory response to SCI.
Furthermore, from a comparative model perspective, Tsujioka and Yamashita (2019) revealed the enrichment of genes signifying the activity of an immature immune system in neonatal mice compared to adult mice groups, and the selective enhanced activity of inflammatory-related pathways in adult mice postpyramidotomy that may suggest their low sprouting ability. Due to higher developmental stage, adult or aged mice exhibit a stronger immune response in comparison to young or neonatal mice, which may be linked to diminished regenerative capacity. Interestingly, however, regeneration capable axolotls and young mice negatively regulate the immune system-related pathways to promote the regenerative/healing process. Those observations may also suggest the possible feasibility to safely modulate the immune system in higher mammals, including humans, for interrogating their regenerative capacity.
STRING analysis at 7 dpi further supported the enrichment analyses from a different angle. In one of the main clusters of this interaction network, downregulated IL-6 and PTGS2 genes were positioned at the center of the cluster and these genes were found as interacting with immune system and extracellular matrix (ECM)-related genes such as CXCR2, COL1A1, ANXA5, and VEGF3. Recently published studies have demonstrated that immune system is actively functional in early period of regeneration, whereas at blastema stage the genes taking roles in immune system-related pathways were found downregulated together with ECM components (Sibai et al., 2020; Tsai et al., 2019). Higher abundance of SMAD4 and MMP14 in injury group compared to sham 1 at 7 dpi may explain the ECM remodeling activities at blastema stage before further steps of regeneration.
Interestingly, several RNA editing genes, such as FBL, TARSL2, NAT10, and NOL11, were found upregulated in SCI-induced animals that may be an indication of the increased RNA modifications during regeneration. Follow-up investigations on modified RNAs during in the timeline of regeneration might be remarkable considering the expanding field of epitranscriptome (Jantsch et al., 2018).
Conclusions
In this study, we aimed to untangle the regenerative capability of metamorphic axolotl after an acute SCI. Surprisingly, despite deficiency in limb regeneration, metamorphic axolotls were successful to functionally restore the removed spinal cord. To explore the molecular signature of spinal cord regeneration at metamorphic stage, we compared the gene expression profile between injury-induced and sham-operated metamorphic axolotls at 1 and 7 dpi. Based on the obtained results, a greater number of DEGs was found at 7 dpi than 1 dpi indicating the regeneration-specific processes such as downregulation of immune system-related processes in injury-induced animals. We believe that functional studies of identified DEGs during spinal cord regeneration would expand the molecular basis of regeneration and provide opportunities to develop new therapeutic strategies for SCI.
Footnotes
Acknowledgments
We thank the Scientific and Technological Research Council of Turkey (TÜBİTAK; Project No. 315S027) for support of the present study.
Author Disclosure Statement
The authors declare they have no competing financial interests.
Funding Information
Scientific and Technological Research Council of Turkey, grant number: 315S027.
Abbreviations Used
References
Supplementary Material
Please find the following supplemental material available below.
For Open Access articles published under a Creative Commons License, all supplemental material carries the same license as the article it is associated with.
For non-Open Access articles published, all supplemental material carries a non-exclusive license, and permission requests for re-use of supplemental material or any part of supplemental material shall be sent directly to the copyright owner as specified in the copyright notice associated with the article.
