Abstract

   Slider turtles are the only known amniotes with self-repair mechanisms
   of the spinal cord that lead to substantial functional recovery. Their
   strategic phylogenetic position makes them a relevant model to
   investigate the peculiar genetic programs that allow anatomical
   reconnection in some vertebrate groups but are absent in others. Here,
   we analyze the gene expression profile of the response to spinal cord
   injury (SCI) in the turtle Trachemys scripta elegans. We found that
   this response comprises more than 1000 genes affecting diverse
   functions: reaction to ischemic insult, extracellular matrix
   re-organization, cell proliferation and death, immune response, and
   inflammation. Genes related to synapses and cholesterol biosynthesis
   are down-regulated. The analysis of the evolutionary distribution of
   these genes shows that almost all are present in most vertebrates.
   Additionally, we failed to find genes that were exclusive of
   regenerating taxa. The comparison of expression patterns among species
   shows that the response to SCI in the turtle is more similar to that of
   mice and non-regenerative Xenopus than to Xenopus during its
   regenerative stage. This observation, along with the lack of conserved
   “regeneration genes” and the current accepted phylogenetic placement of
   turtles (sister group of crocodilians and birds), indicates that the
   ability of spinal cord self-repair of turtles does not represent the
   retention of an ancestral vertebrate character. Instead, our results
   suggest that turtles developed this capability from a non-regenerative
   ancestor (i.e., a lineage specific innovation) that was achieved by
   re-organizing gene expression patterns on an essentially
   non-regenerative genetic background. Among the genes activated by SCI
   exclusively in turtles, those related to anoxia tolerance,
   extracellular matrix remodeling, and axonal regrowth are good
   candidates to underlie functional recovery.

   Keywords: spinal cord injuries, genomics, RNAseq analysis, turtles,
   amniotes, animal models, axon regeneration, glial scarring

Introduction

   Spinal cord injury (SCI) is the cause of one of the most devastating
   conditions for humans because leads to permanent loss of nervous
   functions below the site of injury. Considerable research efforts have
   been made to find new therapies aimed to reduce (or revert) damage or
   to improve the quality of life of patients living with SCI (Silva et
   al., [33]2014). Some groups explored the use of mesenchymal (Park et
   al., [34]2004; Kanno, [35]2013; Li and Lepski, [36]2013) and neural
   stem cells (Sandner et al., [37]2012), whereas others focused on gene
   therapy. The latter approach appears as a promising alternative
   oriented to stimulate intrinsic repairing mechanisms of the injured
   cord (Weidner et al., [38]2012). In this regard, special attention has
   been paid to the modulation of diverse cell growth factors. Among them,
   those from the neurotrophin gene family (particularly BDNF that
   promotes axonal growth) yielded the most promising results in relation
   to therapeutic possibilities. However, currently available gene therapy
   treatments do not achieve a significant functional recovery (Neirinckx
   et al., [39]2014). Moreover, complications from potential treatments,
   such as an increased risk of tumorigenesis (Lee et al., [40]2013), have
   hampered possible clinical applications.

   The inability to regenerate the spinal cord is not restricted to adult
   mammals, but is also observed in birds, most reptiles, and even in
   non-amniotes like frogs (Diaz Quiroz and Echeverri, [41]2013). In
   non-amniote vertebrates, regeneration of the spinal cord is more
   common. Among the latter, urodele amphibians (salamanders) appear as
   one of the most notable examples, as they regenerate several parts of
   their body after amputation, including the spinal cord (Tanaka and
   Ferretti, [42]2009; Diaz Quiroz and Echeverri, [43]2013). In fact,
   2.5–5 months after spinalization, salamanders recover the ability to
   walk and swim almost as uninjured animals (Chevallier et al., [44]2004;
   Cabelguen et al., [45]2010). Some Osteichthyes (bony fishes) can also
   regenerate their spinal cords. Such is the case of the model species
   Danio rerio (Zebrafish), yet its regeneration ability is more limited
   than in salamanders (Becker et al., [46]1997).

   So far, the only known amniotes with spinal cord regenerating capacity
   are the slider turtles, a group of American freshwater turtles commonly
   used as pets (Rehermann et al., [47]2009, [48]2011). The phylogenetic
   location of this group implies that they are the closest relative to
   mammals exhibiting this feature. This evolutionary aspect is an
   essential consideration that makes these turtles a unique model to
   investigate the peculiar cellular and molecular mechanisms that allow
   regeneration or reconnection in some vertebrate groups but not in
   others.

   The cellular and functional processes that take place after traumatic
   spinal injury in turtles were characterized by Rehermann et al.
   ([49]2009) in the species Trachemys dorbignyi. Briefly, the first
   histological indications of healing appear about 5 days after injury,
   consisting in the formation of a cellular bridge permissive for growing
   axons, which connects the cephalic and caudal portions of the injured
   cord. Although, full anatomical restoration and functional recovery was
   never observed, given a proper recovery time, turtles are able to walk
   again using their four legs though stepping is slower than in normal
   animals. Moreover, their hindlimbs are not employed for swimming
   (Rehermann et al., [50]2009). Therefore, this group of turtles
   represents an intermediate stage between non-regenerating vertebrates,
   like mammals, and fully regenerating animals like most anamniote
   vertebrates. In summary, the crucial phylogenetic location of turtles
   (relative proximity to mammals) and their incomplete regeneration
   capacity are strong reasons to consider this group a unique model to
   study spinal cord regeneration.

   Most previous molecular studies on the response to SCI were focused on
   limited sets of genes, which were known or suspected to have a direct
   correlation with the process (van Horssen et al., [51]2006; García et
   al., [52]2012; Demircan et al., [53]2013; Verslegers et al., [54]2013).
   However, considering the complexity of SCI response, both in
   regenerative and non-regenerative vertebrates, it would be also
   important to tackle its study from a genome wide perspective rather
   than focusing only on some potential candidate genes.

   At present, there is no genome-scale study characterizing the response
   to SCI in turtles, and to the best of our knowledge, neither there is
   any global comparative analysis among other vertebrate groups (having
   different or similar regenerating abilities). Besides, genome scale
   studies on mammals are scarce with just some assessments in rodents
   (Velardo et al., [55]2004; Chen K. et al., [56]2013; Torres-Espín et
   al., [57]2013). In this work, we explore the gene expression patterns
   (in a genome wide scale) associated with the response to SCI in
   turtles, aiming to understand the molecular mechanisms that allow
   regeneration in these animals but are absent or inhibited in mammals.

Materials and methods

Animals and surgical procedure

   Fresh-water turtles of the species Trachemys scripta elegans (carapace
   lengths 6 cm) were used (n = 12) in this study. The detailed surgical
   procedures and post-operative care have been described previously
   (Rehermann et al., [58]2009). As control we used “sham injured” animals
   that were obtained using the same surgical procedure as in actual
   injured animals, but leaving the spinal cord intact. This control
   allowed us to “subtract” the side effects of surgery on gene expression
   from those caused by SCI itself. All animal handling and experimental
   procedures were performed according to the animal welfare guidelines of
   CNEA (National Commission for Animal Experimentation). The protocols
   were evaluated and approved (license # 005-09-2015) by the CEUA-IIBCE
   (Comisión de Ética en el Uso de Animales), the Committee for Animal
   Care and Use of Instituto de Investigaciones Biológicas Clemente
   Estable.

Tissue preparation

   For RNASeq differential expression analysis we studied a total of 9
   turtles that were killed 4 days after surgery. In injured turtles, a 4
   mm long cord segment of the spinal cord centered on the lesion
   epicenter was quickly dissected out and placed in RNAlater
   (Sigma-Aldrich, St. Louis, MO, USA). Spinal cord segments of sham
   injured animals were processed in the same way as those with a spinal
   cord transection.

RNA extraction, library preparation, and sequencing

   The tissue was mechanically homogenized using a cordless motor and
   pellet pestle (Sigma–Aldrich). Total RNA was extracted and purified
   using TRIZOL reagent (Invitrogen, Carlsbad, CA, USA) following the
   manufacturer's protocol. The concentration of RNA samples were
   quantified using the Nanodrop (Nano-Drop Inc., Rockland, DE, USA) and
   RNA integrity was determined using Bioanalyzer (Agilent Technologies,
   Santa Clara, CA, USA). Only samples with a RNA Integrity Number (RIN)
   higher than 7 were used to make sequencing libraries. RNA samples were
   polyA selected and sequencing libraries were constructed using
   ScriptSeq RNA Sample Prep Kit as described in the ScriptSeq RNA Sample
   Preparation V2 Guide (Illumina). ScriptSeq Index PCR primers were used
   to allow sequencing several libraries in a single lane. In total, six
   samples were sequenced, two derived from sham animals, three from
   injured ones, and one sample which consisted on a pool of different
   tissues. The latter was used only for transcriptome assembly (see
   Section De novo Transcriptome Assembly). Two replicates, one derived
   from sham and the second from injured animals consisted in a pool of 3
   specimens each and were sequenced using 2 × 100 paired-end reads in an
   Illumina GAIIX platform. Three additional biological replicates (one
   sham and two spinalized animals) were derived from one specimen each,
   and were sequenced using 2 × 75 paired end reads in Illumina MiSeq
   system. RNASeq data used for this study was deposited in the SRA
   repository under the accession number [59]SRP082501.

Gene expression analysis

   The genome of a close relative species, Chrysemys picta, was combined
   with de novo transcriptome assembly (from T. scripta itself) to obtain
   a mapping reference. The assembled transcriptome allowed us to analyze
   genes which are too divergent between the two species (hence preventing
   read mapping on C. picta) or not present in the C. picta genome
   assembly.

Improving the annotation of C. picta (reference) genome using RNAseq data

   Reads were mapped to the C. picta genome (NCBI cpic v 3.0.1) with
   Tophat (v 2.0.8) using default options. Reads that failed to map into
   this genome were isolated and analyzed separately (see Section
   Differential Expression Analysis).

   Mapped reads were assembled using RBTA (Reference Based Transcriptome
   Assembler) with Cufflinks (Trapnell et al., [60]2012; v 2.1.1).
   Assemblies for each condition were combined with the C. picta reference
   genome (NCBI C. picta reference v3.0.1) using Cuffmerge (Trapnell et
   al., [61]2012; v 1.0.0) to produce a final enriched consensus genomic
   reference file (gtf file) that was used for downstream analyses (for
   details see Figure [62]1 and Figure [63]SF4).

Figure 1.

   [64]Figure 1
   [65]Open in a new tab

   Schematic depiction of the experimental pipeline used to generate an
   “enriched” C. picta mapping reference and subsequent differential gene
   expression analysis. Further details in section under the title “de
   novo transcriptome assembly” and Figures [66]S1–S6 in supplementary
   file.

De novo transcriptome assembly

   Only Illumina GAIIX reads (2 × 100 pe) were used for assembling. These
   were quality filtered: undetermined and low quality bases were trimmed
   out, and the resulting read fragments shorter than 70 nt were discarded
   (results are summarized in Figure [67]SF1 and Table [68]S1).

   De novo assembly was performed using Trinity RNA-seq assembler (v
   r2013-05-25) which in the case of our dataset, outperformed other
   assembling algorithms tested (Abyss, Oases, and a two-step strategy
   combining Abyss and Cap3 as described in Supplementary File, in section
   under the title “de novo transcriptome Assembly” and Figure [69]SF2).

   The experimental pipeline of the analysis is depicted in Figure [70]1.

Differential expression analysis

   TopHat2 (Trapnell et al., [71]2009) was used to map reads on the
   enriched C. picta mapping reference obtained with Cuffmerge as
   described in Section Improving the Annotation of C. picta (Reference)
   Genome Using RNAseq Data (i.e., C. picta genome, gtf enriched
   annotation file). When assessing the expression levels of divergent
   genes (or genes not present in the C. picta genome) using the “de novo”
   mapping reference, read mapping was done with Bowtie2 (Langmead and
   Salzberg, [72]2012).

   Read counts were obtained with HTSeq-count (v 0.6.0) and differential
   expression analysis was conducted using DEseq2 (Anders and Huber,
   [73]2010; Love et al., [74]2014). Genes were considered as
   differentially expressed (DE) when the three following conditions were
   met: a value of FDR lower than 0.1 (FDR is the False Discovery Rate, a
   correction of the p-value to account for multiple simultaneous tests);
   a fold change in transcript abundance of at least two (in either
   direction) and a minimum number of mapped reads in at least one of the
   conditions not lower than 50.

Functional characterization of differentially expressed genes: gene ontology
and metabolic pathways analyses

   Because the annotation of C. picta genome publicly available has no
   Gene Ontology (GO) terms associated to its genes, we generated a GO
   annotation using Blast2go (v 2.7.0). Genes containing several isoforms
   may pose biases because these isoforms share GO terms, therefore
   leading to their overrepresentation. To avoid possible artifacts,
   isoforms were collapsed into their corresponding genes using a homemade
   pipeline. Gene Ontology enrichment analysis was done using the Fisher's
   exact test applying FDR correction for multiple tests. Terms with FDR <
   0.05 were denoted as enriched and used in downstream analysis. Enriched
   terms were filtered to extract the most specific ones (i.e., those that
   have higher information content; defined as relative frequency of
   occurrence of the term in the protein annotation dataset under
   consideration, see Lord et al., [75]2003), using the filtering tool
   available in blast2go suite. Analysis of metabolic pathway enrichment
   was conducted using Panther Classification System
   ([76]http://pantherdb.org/).

Protein-protein interaction analysis

   Interactions among the proteins encoded by differentially expressed
   genes were assessed using STRING ([77]http://string-db.org). In order
   to simplify the network and highlight the most biologically relevant
   interactions, we excluded those interactions with a confidence score
   below 0.7. Clustering of the network was obtained using the Louvain
   algorithm (Blondel et al., [78]2008), implemented in the software Pajek
   (v 4.09). The resulting network was exported to Gephi (graph
   visualization and manipulation software; Bastian et al., [79]2009) to
   conduct additional analysis and obtain graphical representations.

Results

Mapping references: reference genome and divergent gene set