- Research article
- Open Access
Using transcriptomics to enable a plethodontid salamander (Bolitoglossa ramosi) for limb regeneration research
BMC Genomicsvolume 19, Article number: 704 (2018)
Tissue regeneration is widely distributed across the tree of life. Among vertebrates, salamanders possess an exceptional ability to regenerate amputated limbs and other complex structures. Thus far, molecular insights about limb regeneration have come from a relatively limited number of species from two closely related salamander families. To gain a broader perspective on the molecular basis of limb regeneration and enhance the molecular toolkit of an emerging plethodontid salamander (Bolitoglossa ramosi), we used RNA-Seq to generate a de novo reference transcriptome and identify differentially expressed genes during limb regeneration.
Using paired-end Illumina sequencing technology and Trinity assembly, a total of 433,809 transcripts were recovered and we obtained functional annotation for 142,926 non-redundant transcripts of the B. ramosi de novo reference transcriptome. Among the annotated transcripts, 602 genes were identified as differentially expressed during limb regeneration. This list was further processed to identify a core set of genes that exhibit conserved expression changes between B. ramosi and the Mexican axolotl (Ambystoma mexicanum), and presumably their common ancestor from approximately 180 million years ago.
We identified genes from B. ramosi that are differentially expressed during limb regeneration, including multiple conserved protein-coding genes and possible putative species-specific genes. Comparative analyses reveal a subset of genes that show similar patterns of expression with ambystomatid species, which highlights the importance of developing comparative gene expression data for studies of limb regeneration among salamanders.
Amphibians belonging to the order Caudata (Urodela) have been studied for more than 100 years in developmental biology and more specifically in the area of tissue regeneration . Studies of relatively few species suggest that salamanders, in general, have a broad capacity to regenerate different tissues and organs, including the heart, brain, jaws, tail and the complex structures of complete limbs . However, recent studies have revealed a surprising degree of interspecific variation in regenerative capacity and in the cellular mechanisms by which regeneration is accomplished. For example, satellite cells serve as the progenitors for reforming muscle during regeneration in the Mexican axolotl (Ambystoma mexicanum; Family Ambystomatidae) and larvae of the eastern red-spotted newt (Notophthalmus viridescens; Family Salamandridae), whereas adult eastern red-spotted newts regenerate using progenitors that are derived from dedifferentiation of functional muscle fibers [3, 4]. These and other examples [5, 6], which are based on relatively few comparative data, suggest that at least some mechanisms of regeneration have diverged during salamander evolution. Clearly, there is need to study additional species to resolve conserved/common mechanisms and taxon-specific differences that have accrued during salamander evolution.
RNA-seq provides an efficient methodology to generate fundamental molecular information for enabling analyses that leverage new model species. To date, transcriptional analyses of limb regeneration have been limited to just a few ambystomatid and salamandrid species (Table 1), with the most transcript data generated for A. mexicanum [7, 8] and N. viridescens . Transcriptomic data were recently generated for the Chinese giant salamander (Andrias davidianus; Family Cryptobranchidae)  and Chinese salamander (Hynobius chinensis; Family Hynobiidae) , but these data were not collected from regenerating tissues. To gain a broader perspective on limb regeneration and salamander evolution in general, we generated transcriptomic data for Bollitoglossa ramosi (Fig. 1a), a South American species from the family Plethodontidae. Plethodontid salamanders diverged from all other salamander families approximately 180 million years ago (162.2–199.0 MYA)  and exhibit several traits not observed in other salamanders, including enucleated red blood cells, projectile tongues, absence of lungs , tail autotomy , nasolabial grooves, and postaxial development of the digits . Also, some plethodontids (including B. ramosi) undergo direct development, wherein individuals hatch from eggs in the adult form and lack a free-living larval phase . While limb regeneration has been investigated in plethodontid salamanders , no study to date has used a transcriptomic approach to globally characterize gene expression.
The transcript data that we report for B. ramosi follow a recent study that described limb regeneration in this species at anatomical and histological levels . Limb regeneration in B. ramosi exhibits morphological similarities to axolotl and newts but also some notable differences, including a more protracted period of blastema growth preceding digit formation, hyperpigmentation of the wound epithelium and collagen deposition in the mesenchymal tissue at 40 days post amputation (dpa). While these differences indicate that regeneration processes are not entirely conserved across divergent salamander families, salamanders, in general, share a blastema-based, developmental strategy to reform amputated limbs. Thus, it seems likely that many core regeneration processes are likely to be shared among salamanders, leading us to perform transcriptome analysis of limb regeneration in B. ramosi. We identified genes from B. ramosi that are differentially expressed during limb regeneration, including multiple conserved protein-coding genes, noncoding RNAs, and putative species-specific genes. Comparative analyses reveal a subset of genes that show similar patterns of expression in ambystomatid species, identifying these as key probes for identifying conserved mechanisms of limb regeneration.
Animals and surgical procedures
All animals used in this work were caught in the wild in non-private owned land, under the Contract on Access to Genetic Resources number 118–2015, which was delivered by the Ministerio del Medio Ambiente (Ministry of Environment) of Colombia to the Principal Investigator and all experimental procedures were approved by the Institutional Bioethics and Animal Care and Use Committee of the University of Antioquia (Medellín, Colombia). Wild adult salamanders (7–10 cm snout to tail length) of the species Bolitoglossa ramosi were collected from their type locality by the nocturnal visual encounter method  in the Andes region of Antioquia, Colombia. Specimens were kept in the laboratory under established protocols .
A total of 22 adults of B. ramosi were administered limb amputations and tissues from the limb, gut, and skin were collected at the time of amputation, and at 20, 40, 60 and 70 days post-amputation (dpa) (Fig. 1b). Briefly, animals were anesthetized by immersion in 1% Tricaine (SIGMA, USA) before bilateral-proximal amputations of the forelimb were performed at mid-humerus level using microscissors and forceps (Fine Precision Tools, USA) . The protruding bone and muscle were trimmed to produce a flat wound surface. Under anesthesia, approximately 2 mm of tissue was collected from the distal tip of amputated limbs. Animals were then euthanized in 2% Tricaine and gut and skin tissues were collected. Overall, we did pools of different animals to get the tissues of controls (D0) (three biological replicates, n = 9 animals), gut (one biological replicate, n = 1 animal), skin (one biological replicate, n = 1 animal), blastema 20 dpa (one biological, replicate n = 1 animal), blastema 40 dpa (one biological replicate, n = 2 animals), blastema 60 dpa (two biological replicates, n = 5 animals) and 70 dpa (one biological replicate, n = 3 animals). Samples were collected and stored in TRIzol® reagent until total RNA was extracted using the protocol from the reagent manufacturer (Life Technologies).
The quality of RNA samples was assessed using an Agilent Technologies 2100 Bioanalyzer. Only samples with an RNA Integrity Number > 8 were used for preparing sequencing libraries. RNA-Seq libraries were prepared from total RNA using TruSeq RNA Sample Prep Kit (Illumina) and the resulting libraries were paired-end (PE) sequenced (2 × 100 bp) using an Illumina Hiseq-2000. The average depth of sequencing for each sample was ~50 million reads (Additional file 1).
De novo transcriptome assembly
The quality of the raw data was assessed using FastQC . Transcripts were assembled using the Trinity (V 2.0.6) software pipeline  and default Trimmomatic parameters to remove sequence adapters and low quality reads (Phred score < 5). In silico reads were further normalized according to the depth of sequence coverage using default settings for Kmer coverage (k = 25). The reference transcriptome was assembled from all RNA samples and all contigs with length ≥ 200 nucleotides were extracted to generate an initial reference transcriptome. Transcriptome assembly quality was assessed based on the calculated E90N50 contig length and BUSCO annotation (Benchmarking Universal Single-Copy Orthologs) .
Gene annotation and GO analysis
Using Reciprocal Best Hits of translation Blast searches (RBH-Blast)  with BLASTx and tBLASTn, B. ramosi contigs were searched against bacterial, viral, single-celled eukaryote, fungal, salamander ribosomal, and salamander mitochondrial sequence databases compiled from NCBI in order to identify potential contaminants. Sequences with sequence identity ≥50% and bit scores values ≥50 were removed from the reference transcriptome (N = 4217). We further predicted long open reading frames (ORF) using TransDecoder (version 3.0.0) software . For gene annotation, B. ramosi contigs were reciprocally Blast searched against vertebrate sequence databases to identify high identity alignments that could be used to infer homology relationships. Blast searches were performed using translated nucleotide sequences of A. mexicanum [25,26,27] and N. viridescens , and protein-coding sequences from seven vertebrate taxa (Anolis carolinensis: GCA_000090745.1, Danio rerio: GCA_000002035.3, Gallus gallus: PRJNA10808, Homo sapiens: PRJNA168, Mus musculus: GCA_000001635.7, Xenopus tropicalis: PRJNA205740, Latimeria chalumnae: GCA_000225785.1; available through Ensembl or Refseq NCBI). Gene names were assigned to B. ramosi contigs if an alignment showed ≥50% identity and returned a bit scores value ≥50.
Transcript abundance (RSEM) and expression level analyses (EBSeq)
Sequence reads generated from limb tissue samples were aligned to the reference transcriptome using Bowtie2  and RSEM (RNA-Seq by Expectation Maximization) was used to obtain estimates of transcript abundance for all transcripts . Expression levels were calculated as transcripts per million (TPM). For cases where two or more significant transcripts mapped to the same gene identifier, only the transcript with the highest expression estimate was evaluated statistically using EBseq. Transcripts were considered differentially expressed between control and post-amputation samples when TPM was ≥0.95 (4 to 800,000 expected counts) for at least a single time point and fold change (log2FC) was ≤ − 2 and ≥ 2 with an FDR p < 0.05 (Fold discovery rate). Differentially expressed genes were further analyzed using PANTHER gene expression tools (Version 11.1) to identify corresponding Gene Ontology (GO) terms .
Three replicate tissue samples from collected from unamputated and regenerating (40 and 60 days post amputation) forelimbs and total RNA was prepared using TRIzol® reagent (Life Technologies). The total RNA was reverse-transcribed to single-stranded cDNA with reverse transcriptase (Thermo) in the presence of random hexamer primers, oligoDT primers, and dNTPs for 60 min at 42 °C. Expression levels of specific mRNAs were determined by qPCR using gene-specific primer pairs (Additional file 2) (two technical replicates). Each reaction was performed at a total volume of 10 μL containing 50 ng first-strand cDNA, 5 μL Syber greenMix (Biorad), and 0.1 μM of each primer pair, and cycled on a Biorad Real-Time PCR system. Real-time data were analyzed using Biorad software version 2.1. Relative mRNA expression was calculated using the 2 –ΔΔCT method with GAPDH as a cross-sample reference. Significance was determined using a two-tailed Unpaired T-test and Welch’s correction (P < 0.05). Pearson correlation (P < 0.05) was performed to analyze the correlation between the log2 Fold change of 40dpa and 60dpa (relative to the control unamputated sample) obtained in the in silico analyses by RNA-Seq and with the validation by RT-qPCR.
A de novo reference transcriptome of B. ramosi was generated from RNAs that were isolated from normal limb, regenerating limb, skin, and gut. The total number of high quality assembled PE reads recovered was 654,673,506. Using Trinity software, we obtained 433,809 contigs with an average GC content of 43.5%, an average length of 569 bp, and a maximum assembled contig length of 20,709 bp (Table 2). On the basis of read coverage, the E90N50 statistic was ~3 Kb (Additional file 3) and the reference transcriptome contained 85.1% of the conserved core eukaryotic genes using BUSCO annotation (Additional file 4). Blast searches of B. ramosi contigs identified 36 mitochondrial, zero rDNA, and 4181 putative microorganism hits. Notably, multiple hits (n = 66) matched sequences from the pathogenic chytrid fungus Batrachochytrium dendrobatides. This fungus is a primary cause of worldwide declines in amphibian populations  and was first described in Colombia in 2013 .
The contigs of B. ramosi were annotated by performing BLAST searches against A. mexicanum transcriptome databases (Ambystoma.org and Axolotl-omics.org) and N. viridescens, and by generating a cross-referenced dataset of orthologous genes between A. mexicanum and B. ramosi (n = 13,065) (Additional file 5). Large divergences were expected between salamanders species (180 MYA), the frequencies of the % identity between B. ramosi vs N. viridescens and B. ramosi vs A. mexicanum during the RHB-Blast showing than the major of the matches were > 80% (Additional file 6). We also searched Ensembl and NCBI protein sequences from seven vertebrate taxa, which identified gene annotations for 26,183 B. ramosi contigs (Additional file 7). Finally, we obtained homology information for contigs that did not return a sequence match by searching Treefam , UniRef90 , PFAM  (Additional file 8) and ncRNA (miRBase and RFam) [34, 35] (Additional file 9) databases. These searchers (Additional file 10) returned non-redundant annotations for 140,974 of B. ramosi contigs. Across all search strategies (Fig. 2), we obtained gene and RNA annotations for 142,926 non-redundant transcripts of the B. ramosi reference transcriptome.
Differential gene expression analysis
Due to limitations in obtaining source materials from wild B. ramosi, our study included three replicate samples for the D0 and two replicates for 60 dpa time points, but only single replicate samples for 20, 40, and 70 dpa time points. Although a previous transcriptional study of salamander limb regeneration reported genes as significant from contrasts of single replicates , we conservatively only report genes that exhibited statistically significant changes in expression between the control and two or more post-amputation time points, and where the magnitude of the expression difference was > 2-fold (Fig. 3). We reasoned that this approach reduced the number of false positives and differentially expressed genes associated with individual-specific differences (e.g. age, disease history, physiological state) that are unrelated to limb regeneration. Overall, 602 non-redundant differentially expressed genes were identified, of which 556 were assigned an official gene symbol or locus identifier from a model organism database (Additional file 11). Hierarchical clustering identified two primary clusters, one where genes were expressed more highly in regenerating limbs than controls across a majority of time points (n = 310) and another where genes were expressed more highly in controls (n = 292) (Fig. 3a). The cluster of highly expressed genes in regenerating limbs included extracellular matrix (different collagen isoforms, emilin1, adamts4, spon2, adamts17, mfap2, clec11a, p3h1) and developmental process (lef1, msx1, sox12, sall4) genes. The cluster of lowly expressed genes in regenerating samples encode transcription factors (nfkbiz, fhl1, znf385a, tsc22d1, nfil3, sfmbt2, fhl2, pdlim7, nr1d1, foxk2, per1, tfcp2, elf5, znf618, en1, ubp1, foxq1), and extracellular matrix components (ogn, mgp, mfap5, slitrk6, crim1, ap1m1, dpt).
Genes were primarily down-regulated during the earliest post-amputation time point (20 dpa) (Fig. 3b). These genes are associated with collagen biosynthesis (Col1a1, Col6a2, Col11a1, Col12a1, Col6a1, Serpinh1, Col4a5, Col6a3), axon guidance (Epha4, Nrp2, St8sia2, Kif4a, Ncam1, Shb, Dnm1, Robo1) and developmental regulation (Hoxa1, Ezh2). Genes that were up-regulated at 20 dpa showed decreased expression at subsequent sample times, and this list included genes that function in actin crosslinking (GO:0051764) (Rac2, Lcp1), regulation of the p38MAPK cascade (GO:1900744) (Per1), and hemidesmosome assembly (GO:0031581) (Krt14).
We note that 109 of the 602 genes were differentially expressed at all post-amputation time points in comparison to the non-amputated condition (Additional file 12), and many of these encode proteins for development processes (45.8%), including prrx1, sall4, snai2, sox12, pak1, lef1, lefty2, and smarca4.
Some of the genes exhibited complex patterns of gene expression. For example, some genes were expressed differently in regenerating limbs relative to controls at 20 dpa and 40 dpa; these included collagen (col11a1, col1a1, col6a1) and homeobox (prrx1, hoxa1) genes that were expressed lower than controls at 20 dpa but higher at 40 dpa (Fig. 3b). Also, some genes were more highly expressed in regenerating samples than controls specifically at 40 dpa (Table 3, Fig. 3c), including transcription factors (msx1, lmx1a), cell adhesion molecules (tgfbi), extracellular matrix components (emilin1, adamts4, spon2, adamts17, mfap2, clec11a, p3h1), cytoskeletal molecules (tubb, kifc1, tubb2b, lmx1a) and developmental process genes (fgd1, adamts4, spon2, rbp7, lef1, clec11a, alpi, msx1, lmx1a, dact2).
The annotated gene list (N = 556) was subjected to a statistical over-representation test using Panther Gene List Analysis tools and default settings  (Additional file 13). We identified enriched biological process gene ontologies associated with extracellular matrix components, cell differentiation, migration, proliferation, morphogenesis, and development of multiple tissues (epithelium, neurons, vasculature, cartilage, and bone). Many of these genes encode proteins that function in cell signaling pathways associated with tissue development, including Wnt (wnt11, wnt5a, tcf7l1, wnt4, fzd1, lef1, en1, tp53, smarca4), TGFb/BMP (fosl1, smad7, bambi, bmp4, junb), Hedgehog (gli3), Hox (hoxa1), and Notch (pofut1, hes1, hes4) (Fig. 3b, c).
The significant B. ramosi DEGs were further compared against 3053 significant genes identified by Voss et al. , the most comprehensive and statistically powered study of A. mexicanum limb regeneration to date. In that study, an Affymetrix microarray was used and thus estimates were obtained for a finite number of genes. We determined that the  study generated A. mexicanum transcript abundance estimates for 395 of the 602 B. ramosi DEG genes. The majority (71%; N = 282) of these A. mexicanum genes were significantly differentially expressed in the  study and showed a similar temporal pattern of expression relative to orthologous B. ramosi DEGs (Fig. 4). The DEGs shared between B. ramosi and A. mexicanum were enriched for GO terms associated with protein binding (GO:0005515) (sall4, sox12, fhl1, pdlim7), extracellular matrix organization (GO:0030198) (vcan, col6a2, col11a1, col6a1, tnc, col6a3, mmp1) and developmental growth (GO:0048589) (wnt5a, aurka, sema3f) (Additional file 14). Also, we found that 104 genes had an ortholog in the axolotl-omics.org database but not in the Voss et al. study (2015) and some genes (n = 103) had no ortholog in any of the axolotl databases (Additional file 15).
Finally, we recovered a group of transcripts (n = 85) that had good transcriptional support (TPM ≥ 0.95) but did not match sequences from any of the reference databases that were searched. From this group of anonymous transcripts, we identified 11 potentially novel, taxon-specific genes with complete ORFs that were differentially expressed during regeneration. We also identified 14 hits with homology to ncRNAs (Additional file 16). These genes were found to have good transcriptional evidence (TPM ≥ 0.95) and were differentially expressed in some of the conditions evaluated. These candidates included matches to ncRNAs such as eca-mir-9182 (MI0028412), signal recognition particle RNA (RF00017), and 7SK (RF00100) (Fig. 3d).
Validation by RT-qPCR
We made use of RT-qPCR to validate the expression estimates obtained by RNA-Seq for selected transcripts (sall4, fn1, myot, coll11a1, col1a1, col6a1) (Fig. 5, Additional file 17) that are known to be modulated during limb regeneration [25, 36, 37]. Also, we validated the expression of two B. ramosi taxon-specific transcripts. For these validations, biological replicate RNA samples were prepared for control limbs, and for regenerating forelimb tissues collected 40 and 60 dpa. We observed that sall4, fn1, col11a1, col1a1, col6a1 (higher in regenerating limbs) and myot (lower in regenerating limbs) showed the same pattern of expression in qPCR as was found using RNA-Seq. Also, the transcripts that we defined as taxon-specific showed the same patterns using both methodologies (Additional file 18).
Here we present the first de novo reference transcriptome of limb regeneration for B. ramosi, a salamander belonging to the largest family of salamanders. Plethodontids of the genus Bolitoglossa exhibit some notable biological differences in comparison to model salamanders (axolotls and newts) from other taxonomic families, including their restrictive terrestrial habits and direct development without larval stages [14, 16]. We recently reported anatomical and morphological changes during limb regeneration in B. ramosi . Many of the changes that we observed, including processes associated with wound healing, blastemal formation, and re-patterning, also occur during limb regeneration in other salamanders. These features are likely conserved among all salamanders and indeed we identified DEGs in B. ramosi that are known to be expressed during limb regeneration in A. mexicanum and N. virdescens [36,37,38]. These include genes and proteins that associate with the wound epidermis, extracellular matrix, basement membrane, blastema and differentiating chondrogenic precursor cells. Additionally, we identified 109 genes that were differentially expressed throughout regeneration in B. ramosi, of which 77 were previously identified as differentially expressed during limb regeneration in A. mexicanum [25, 39, 40,41,42]. However, we found genes that were DEG in B. ramosi but have not been reported yet as DEG in A. mexicanum and N. virdescens. This is the case for TGF-beta signaling pathway genes (Lefty2 and Bmp4). The gene Lefty2 was highly expressed during all time points of regeneration, including at 40 dpa, when the blastema shows increased pigmentation. This gene is implicated in left-right axis determination during development, however, it will be important to further analyze the Lefty2 function during limb regeneration. Also, Bmp4 function is important during limb development in mice  and direct development of frogs, where it provides an important biomarker of skeletogenic cell differentiation [43, 44]. This gene was up-regulated in B. ramosi at 70 dpa, which correlates to the palette stage of limb regeneration when limb skeletal patterning occurs. Thus it is possible that this gene has an important role during regeneration of salamanders with direct development. Collectively, these results suggest an important role for that the TGF-beta signaling pathway during limb regeneration in B. ramosi.
Through additional gene expression profile comparisons between B. ramosi and A. mexicanum, we observed conservation of temporal gene expression patterns (Fig. 4). Many genes that were up-regulated or down-regulated throughout limb regeneration in B. ramosi showed the same expression pattern in A. mexicanum, and this correlation became even greater when ignoring changes in gene expression at early post-amputation time points in A. mexicanum that were not sampled for B. ramosi. We note though that it is difficult to reconcile gene expression similarity for some genes, even when considering differences in temporal sampling. For example, col6a1, col6a2, and col6a3 show different patterns of expression between A. mexicanum and B. ramosi, however, their correlated expression within species strongly suggests conservation of mechanisms that regulate transcription during limb regeneration.
Although our results generally support the idea that many aspects of the limb regeneration process are highly conserved among salamander limb regeneration programs, regeneration differences [4, 6, 45, 46] are likely to evolve among lineages that present different modes of development  and life history [48, 49]. For example, ECM (Extracellular matrix) proteins are important modulators of cell proliferation, migration, and differentiation, and they also provide structural support to regenerating limbs that is important for blastemal cell survival [50, 51]. The ECM is predicted to be more rigid in regenerating limbs of terrestrial salamanders that support body weight during locomotion. Less rigidity is predicted for aquatic salamanders that use buoyancy to reduce load-bearing on limbs during regeneration, and swimming for locomotion. Consistent with these predictions, we observed high expression of collagen genes during B. ramosi limb regeneration and this correlated with a high abundance of collagen protein in the blastema . Generation of a more rigid ECM, perhaps to minimize compressive forces on blastemal cells, may be characteristic of terrestrial salamanders that continue to walk during their protracted, limb regeneration programs. Also, we identified DEGs at 20 dpa that further support this idea; the proteins encoded by Rac2, Lcp1, and Krt14 might confer mechanical support to the epithelium and protect the early blastema.
We also note the possibility that some features unique to B. ramosi may be associated with taxon-specific genes, a feature that is also known for other salamanders [52,53,54]. Our analyses identified 12 presumptive taxon-specific genes, two of which were validated using RT-qPCR for control limbs, and 40dpa and 60dpa time points (Fig. 5). All of these transcripts had ORFs > 1000 bp, and 3D prediction with I-TASSER  identified peptide sequences consistent with protein folding (primarily alpha helices). It is possible that these represent novel, taxon-specific proteins that could play an important role in B. ramosi limb regeneration. Future work will endeavor to determine the function of these genes in vivo.
We would further like to acknowledge that some limitations were encountered in the use of a wild-caught species as a research model. Mandatory environmental licensing was required for the collection of wild specimens, and while these regulations are in the best interests of preserving biodiversity, specimen sampling may be reduced below thresholds needed to perform a high-powered statistical analysis. Indeed, such conditions may present an unfortunate hindrance of biodiversity research in developing countries . The use of wild-caught animals is additionally complicated where standardization of optimal captive conditions may be unknown or difficult to replicate, and where suboptimal conditions place the animals under sustained stress leading to disease and decreased survivorship. Further study of the reproductive biology of B. ramosi is needed to establish a captive breeding program and in turn to grant experimental access to more animals of different stages to further test hypotheses generated during this study.
In this study, we report the first de novo reference transcriptome during limb regeneration in a plethodontid salamander, which allowed the identification of differentially expressed genes. Our study shows conservation of transcriptional regulation across plethodontids and ambystomatids which diverged some 180 Ma ago . The genes that we report should be especially good biomarkers for future comparative studies of limb regeneration among salamander taxa.
Benchmarking universal single-copy orthologs
Day zero, control limb
Differential expressed genes
Days post amputation
Fold discovery rate
Million years ago
Open reading frames
Reciprocal best hits of translation blast searches
RNA-seq by expectation maximization
Transcripts per million
Stocum DL, Cameron JA. Looking proximally and distally: 100 years of limb regeneration and beyond. Dev Dyn. 2011;240:943–68.
Zhao A, Qin H, Fu X. What determines the regenerative capacity in animals? Bioscience. 2016;66:735–46.
Sandoval-Guzmán T, Wang H, Khattak S, Schuez M, Roensch K, Nacu E, et al. Fundamental differences in dedifferentiation and stem cell recruitment during skeletal muscle regeneration in two salamander species. Cell Stem Cell. 2014;14:174–87.
Sugiura T, Wang H, Barsacchi R, Simon A, Tanaka EM. MARCKS-like protein is an initiating molecule in axolotl appendage regeneration. Nature. 2016;531:237–40.
Sousounis K, Bhavsar R, Looso M, Krüger M, Beebe J, Braun T, et al. Molecular signatures that correlate with induction of lens regeneration in newts: lessons from proteomic analysis. Hum Genomics. 2014;8:22.
Eguchi G, Eguchi Y, Nakamura K, Yadav MC, Millán JL, P a T. Regenerative capacity in newts is not altered by repeated regeneration and ageing. Nat Commun. 2011;2:384.
Bryant DM, Johnson K, Ditommaso T, Regev A, Haas BJ, Whited JL. A tissue-mapped axolotl De novo transcriptome enables identification of limb regeneration factors. Cell Rep. 2017;18:762–76.
Baddar NWAH, Woodcock MR, Khatri S, Kump DK, Voss SR. Sal-site: research resources for the Mexican axolotl. Methods Mol Biol. 2015;1290:321–36.
Looso M, Braun T. Data mining in newt-omics, the repository for omics data from the newt. Methods Mol Biol. 2015;1290:337–51.
Jiang X, Wang Y, Zhang X. Data set for transcriptome analysis of the Chinese giant salamander (Andrias davidianus). Data Br. 2016;6:12–4.
Che R, Sun Y, Wang R, Xu T. Transcriptomic analysis of endangered Chinese salamander: identification of immune, sex and reproduction-related genes and genetic markers. PLoS One. 2014;9:e87940.
Kumar S, Stecher G, Suleski M, Hedges SB. TimeTree: a resource for timelines, Timetrees, and divergence times. Mol Biol Evol. 2017;34:1812–9.
Wake DB. What salamanders have taught us about evolution. Annu Rev Ecol Evol Syst. 2009;40:333–52.
Mueller RL, Macey JR, Jaekel M, Wake DB, Boore JL. Morphological homoplasy, life history evolution, and historical biogeography of plethodontid salamanders inferred from complete mitochondrial genomes. Proc Natl Acad Sci U S A. 2004;101:13820–5.
Wake DB, Hanken J. Direct development in the lungless salamanders: what are the consequences for developmental biology, evolution and phylogenesis? Int J Dev Biol. 1996;40:859–69.
Chippindale PT, Bonett RM, Baldwin AS, Wiens JJ. Phylogenetic evidence for a major reversal of life history evolution in plethodontid salamanders. Evolution. 2004;58:2809–22.
Scadding SR. Limb regeneration in adult amphibia. Can J Zool. 1981;59:34–46.
Arenas Gómez CM, Gomez Molina A, Zapata JD, Delgado JP. Limb regeneration in a direct-developing terrestrial salamander, Bolitoglossa ramosi (Caudata: Plethodontidae). Regeneration. 2017;4:227–35.
Grover MC. Comparative effectiveness of nighttime visual encounter surveys and cover object searches in detecting salamanders. Herpetol Conserv Biol. 2006;1:93–9.
Arenas CM, Gómez-Molina A, Delgado JP. Maintaining plethodontid salamanders in the laboratory for regeneration studies. Methods Mol Biol. 2015;1290:71–8.
Simon Andrews. FastQC: a quality control tool for high throughput sequence data. Http://www.bioinformatics.babraham.ac.uk/projects/fastqc 2010.
Haas BJ, Papanicolaou A, Yassour M, Grabherr M, Blood PD, Bowden J, et al. De novo transcript sequence reconstruction from RNA-Seq using the trinity platform for reference generation and analysis. Nat Protoc. 2013;8:1494–512.
Simão FA, Waterhouse RM, Ioannidis P, Kriventseva EV, Zdobnov EM. BUSCO: assessing genome assembly and annotation completeness with single-copy orthologs. Bioinformatics. 2015;31:3210–2.
Moreno-Hagelsieb G, Latimer K. Choosing BLAST options for better detection of orthologs as reciprocal best hits. Bioinformatics. 2008;24:319–24.
Caballero-Pérez J, Espinal-Centeno A, Falcon F, García-Ortega LF, Curiel-Quesada E, Cruz-Hernández A, et al. Transcriptional landscapes of Axolotl (Ambystoma mexicanum). Dev Biol. 2018;433:227–39.
Voss SR, Palumbo A, Nagarajan R, Gardiner DM, Muneoka K, Stromberg AJ, et al. Gene expression during the first 28 days of axolotl limb regeneration I: experimental design and global analysis of gene expression. Regeneration. 2015;2:120–36.
Axolotl-omics website. www.axolotl-omics.org. Accesed 25 Mar 2018.
Langmead B, Salzberg SL. Fast gapped-read alignment with bowtie 2. Nat Methods. 2012;9:357–9.
Li B, Dewey CN. RSEM: accurate transcript quantification from RNA-Seq data with or without a reference genome. BMC Bioinformatics. 2011;12:323.
Mi H, Huang X, Muruganujan A, Tang H, Mills C, Kang D, et al. PANTHER version 11: expanded annotation data from gene ontology and Reactome pathways, and data analysis tool enhancements. Nucleic Acids Res. 2017;45:D183–9.
Olson DH, Aanensen DM, Ronnenberg KL, Powell CI, Walker SF, Bielby J, et al. Mapping the global emergence of Batrachochytrium dendrobatidis, the amphibian Chytrid fungus. PLoS One. 2013;8:e56802.
Flechas SV, Medina EM, Crawford AJ, Sarmiento C, Cárdenas ME, Amézquita A, et al. Characterization of the first Batrachochytrium dendrobatidis isolate from the Colombian Andes, an amphibian biodiversity hotspot. Ecohealth. 2013;10:72–6.
Li H, Coghlan A, Ruan J, Coin LJ, Hériché J-K, Osmotherly L, et al. TreeFam: a curated database of phylogenetic trees of animal gene families. Nucleic Acids Res. 2006;34:D572–80.
Suzek BE, Wang Y, Huang H, McGarvey PB, Wu CH, UniProt Consortium the U. UniRef clusters: a comprehensive and scalable alternative for improving sequence similarity searches. Bioinformatics. 2015;31:926–32.
Punta M, Coggill PC, Eberhardt RY, Mistry J, Tate J, Boursnell C, et al. The Pfam protein families database. Nucleic Acids Res. 2012;40:D290–301.
Stewart R, Rascón CA, Tian S, Nie J, Barry C, Chu L-F, et al. Comparative RNA-seq analysis in the Unsequenced axolotl: the oncogene burst highlights early gene expression in the Blastema. PLoS Comput Biol. 2013;9:e1002936.
Kozomara A, Griffiths-Jones S. miRBase: annotating high confidence microRNAs using deep sequencing data. Nucleic Acids Res. 2014;42:D68–73.
Mercer SE, Cheng C-H, Atkinson DL, Krcmery J, Guzman CE, Kent DT, et al. Multi-tissue microarray analysis identifies a molecular signature of regeneration. PLoS One. 2012;7:e52375.
Burge SW, Daub J, Eberhardt R, Tate J, Barquist L, Nawrocki EP, et al. Rfam 11.0: 10 years of RNA families. Nucleic Acids Res. 2013;41:D226–32.
Rao N, Jhamb D, Milner DJ, Li B, Song F, Wang M, et al. Proteomic analysis of blastema formation in regenerating axolotl limbs. BMC Biol. 2009;7:83.
Santosh N, Windsor LJ, Mahmoudi BS, Li B, Zhang W, Chernoff EA, et al. Matrix metalloproteinase expression during blastema formation in regeneration-competent versus regeneration-deficient amphibian limbs. Dev Dyn. 2011;240:1127–41.
Monaghan JR, Athippozhy A, Seifert AW, Putta S, Stromberg AJ, Maden M, et al. Gene expression patterns specific to the regenerating limb of the Mexican axolotl. Biol Open. 2012;1:937–48.
Bandyopadhyay A, Tsuji K, Cox K, Harfe BD, Rosen V, Tabin CJ. Genetic analysis of the roles of BMP2, BMP4, and BMP7 in limb patterning and Skeletogenesis. PLoS Genet. 2006;2:e216.
Kerney R, Gross JB, Hanken J. Early cranial patterning in the direct-developing frog Eleutherodactylus coqui revealed through gene expression. Evol Dev. 2010;12:373–82.
Monaghan JR, Epp LG, Putta S, Page RB, Walker JA, Beachy CK, et al. Microarray and cDNA sequence analysis of transcription during nerve-dependent limb regeneration. BMC Biol. 2009;7:1.
Wu C-H, Tsai M-H, Ho C-C, Chen C-Y, Lee H-S. De novo transcriptome sequencing of axolotl blastema for identification of differentially expressed genes during limb regeneration. BMC Genomics. 2013;14:434.
Sousounis K, Athippozhy AT, Voss SR, Tsonis PA. Plasticity for axolotl lens regeneration is associated with age-related changes in gene expression. Regeneration. 2014;1:47–57
Frobisch NB, Bickelmann C, Witzmann F. Early evolution of limb regeneration in tetrapods: evidence from a 300-million-year-old amphibian. Proc R Soc B Biol Sci. 2014;281:20141550–20141550.
Fröbisch NB, Shubin NH. Salamander limb development: integrating genes, morphology, and fossils. Dev Dyn. 2011;240:1087–99.
Godwin J, Kuraitis D, Rosenthal N. Extracellular matrix considerations for scar-free repair and regeneration: insights from regenerative diversity among vertebrates. Int J Biochem Cell Biol. 2014;56C:47–55.
Meredith JE, Fazeli B, Schwartz MA. The extracellular matrix as a cell survival factor. Mol Biol Cell. 1993;4:953–61.
da Silva SM, Gates PB, Brockes JP. The newt Ortholog of CD59 is implicated in Proximodistal identity during amphibian limb regeneration. Dev Cell. 2002;3:547–55.
Kumar A, Godwin JW, Gates PB, Garza-Garcia a A, Brockes JP. Molecular basis for the nerve dependence of limb regeneration in an adult vertebrate. Science. 2007;318:772–7.
Smith JJ, Putta S, Zhu W, Pao GM, Verma IM, Hunter T, et al. Genic regions of a large salamander genome contain long introns and novel genes. BMC Genomics. 2009;10:19.
Roy A, Kucukural A, Zhang Y. I-TASSER: a unified platform for automated protein structure and function prediction. Nat Protoc. 2010;5:725–38.
Prathapan KD, Pethiyagoda R, Bawa KS, Raven PH, Rajan PD. When the cure kills—CBD limits biodiversity research. Science. 2018;360:1405–6.
Smith JJ, Putta S, Walker JA, Kump DK, Samuels AK, Monaghan JR, et al. Sal-site: integrating new and existing ambystomatid salamander research and informational resources. BMC Genomics. 2005;6:181.
Monaghan JR, Walker JA, Page RB, Putta S, Beachy CK, Voss SR. Early gene expression during natural spinal cord regeneration in the salamander Ambystoma mexicanum. J Neurochem. 2007;101:27–40.
Makarev E, Call M, Grogg M, Atkinson D. Gene expression signatures in the newt irises during lens regeneration. FEBS Lett. 2007;581:1865–70.
Maki N, Martinson J, Nishimura O, Tarui H, Meller J, Tsonis PA, et al. Expression profiles during dedifferentiation in newt lens regeneration revealed by expressed sequence tags. Mol Vis. 2010;16:72–8.
Campbell L. Gene expression profile of the regeneration epithelium during axolotl limb regeneration. Dev Dyn. 2011;240:1826–40.
Holman EC, Campbell LJ, Hines J, Crews CM. Microarray analysis of microRNA expression during axolotl limb regeneration. PLoS One. 2012;7:e41804.
Sousounis K, Michel CS, Bruckskotten M, Maki N, Borchardt T, Braun T, et al. A microarray analysis of gene expression patterns during early phases of newt lens regeneration. Mol Vis. 2013;19:135–45.
Looso M, Preussner J, Sousounis K, Bruckskotten M, Michel CS, Lignelli E, et al. A de novo assembly of the newt transcriptome combined with proteomic validation identifies new protein families expressed during tissue regeneration. Genome Biol. 2013;14:R16.
Abdullayev I, Kirkham M, Björklund Å, Simon A, Sandberg R. A reference transcriptome and inferred proteome for the salamander Notophthalmus viridescens. Exp Cell Res. 2013;319:1187–97.
Nakamura K, Islam MR, Takayanagi M, Yasumuro H, Inami W, Kunahong A, et al. A transcriptome for the study of early processes of retinal regeneration in the adult newt, Cynops pyrrhogaster. PLoS One. 2014;9:e109831.
Elewa A, Wang H, Talavera-López C, Joven A, Brito G, Kumar A, et al. Reading and editing the Pleurodeles waltl genome reveals novel features of tetrapod regeneration. Nat Commun. 2017;8:1–9.
Nowoshilow S, Schloissnig S, Fei J-F, Dahl A, Pang AWC, Pippel M, et al. The axolotl genome and the evolution of key tissue formation regulators. Nature. 2018;554:50–5.
We thank the lab of Juan Fernando Alzate from the University of Antioquia for their help in the bioinformatic methodological approach in the beginning of this study, to Diego Fernanado Uribe and Carlos Muskus in the RT-qPCR analysis. We thank Andrea Gómez and Melisa Hincapie for their help in animal collection and husbandry.
This work was funded by grants from the University of Antioquia (CODI) and COLCIENCIAS (569,027-2013) to JPD and CMAG (COLCIENCIAS 567), the grants allowed the design of the study and collection, analysis, interpretation of data and writing of the manuscript. Grants from the National Institutes of Health (R24OD010435) and Army Research Office (W911NF1410165) supported computational resources, bioinformatics training, data analysis, data interpretation, and manuscript writing efforts of MRW, SRV, and JJS at the University of Kentucky.
Availability of data and materials
The raw sequence reads have been deposited in the Sequence Read Archive under the accession number SRP120553. The output from RSEM, EBSeq and assembled transcriptome are deposited in The Gene Expression Omnibus (GEO) under the accession number GSE105232.
The Andes Community Environmental Regulations, and also the Andes bill (Decision Adina 391 de 1996), stipulate that the biodiversity, animals and all molecular resources are property of the Colombian Government. The salamanders used in this work were caught in the wild in non-private owned land, under the Contract on Access to Genetic Resources number 118–2015, which was delivered by the Ministerio del Medio Ambiente (Ministry of Enviroment) of Colombia to the Principal Investigator. All experimental procedures were approved by the Institutional Bioethics and Animal Care and Use Committee of the University of Antioquia (Medellín, Colombia).
Consent for publication
The authors declare that they have no competing interests.
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Transcript sequence depth for samples used in this study. (DOCX 11 kb)
RT-qPCR primers used for validation of gene expression during regeneration in Bolitoglossa ramosi. (DOCX 12 kb)
E90N50 statistic for de novo reference transcriptome of Bolitoglossa ramosi . The reference transcriptome has an E90N50 of ~3kb (red arrow). (TIF 221 kb)
BUSCO Analysis of Transcriptome Completeness. (DOCX 11 kb)
Output of BLASTn alignments between the Bolitoglossa ramosi reference transcriptome, Ambystoma mexicanum and Notophthalmus viridescens databases. (XLSX 16976 kb)
Histogram that shows the frequencies of the % identity between B. ramosi vs N. viridescens and B. ramosi vs A. mexicanum during the RHB-Blast. (TIF 3750 kb)
Output of the RBH-BLAST alignments of the Bolitoglossa ramosi reference transcriptome using seven vertebrate databases. (XLSX 36198 kb)
Output from the BLASTp alignments of the Bolitoglossa ramosi predicted ORFs against TreeFam, PFAM and UniRef90 databases. (XLSX 73915 kb)
Output from the BLASTn alignments of the Bolitoglossa ramosi reference transcriptome against the ncRNA database. (XLSX 950 kb)
Homology assignments recovered in de novo reference transcriptome assembly of Bolitoglossa ramosi . The B. ramosi transcriptome was surveyed by Reciprocal Best Hits of translated BLAST searches (RBH-BLAST) to protein or translated databases from different vertebrates. Additional gene family homologs were assigned to B. ramosi using protein BLAST against the UniRef90, TreeFam and PFAM domain databases, as well as BLASTN against ncRNA databases. (TIF 48 kb)
DEGs identified from Bolitoglossa ramosi with ≥2-fold (Log2) difference relative to the unamputated control at two or more post-amputation time points, and with TPM support ≥0.95. (XLSX 227 kb)
List of 109 genes that were differently expressed at all post-amputation time points in comparison to the non-amputated condition. (XLSX 31 kb)
Statistical over-representation test of biological processes using Panther Gene List Analysis tools of the annotated gene list. (XLSX 13 kb)
Enriched GO associated terms for DEGs shared between Bolitoglossa ramosi and A. mexicanum. (XLSX 18 kb)
List of DEG of B. ramosi not reported yet as differential expressed genes during limb regeneration in A. mexicanum. (XLSX 8 kb)
Transcripts annotated as possible ncRNA which were differentially expressed during limb regeneration in Bolitoglossa ramosi. (XLSX 8 kb)
RT-qPCR output and comparison of RT-qPCR and RNA-Seq estimates of gene expression. (XLSX 26 kb)
Correlation between RNA-Seq and quantitative real time PCR (RT-qPCR). The expression patterns are similar between the RNA-Seq and RT-qPCR data, and statistically significant Pearson correlation is shown for the expression levels of eight genes measured with both methodologies during 40 dpa and 60 dpa of limb regeneration. The X and Y axis show log2 fold change for RT-qPCR and RNA-Seq values, respectively. (TIF 133 kb)