BACKGROUND
echinoderms constitute a phylum of marine invertebrates closely related to chordates.
their phylogenetic position as a non-chordate deuterostome phylum combined with the ability to regenerate various body parts makes them a valuable group, which can give unique insights into fundamental issues of regenerative biology and provide new clues to finding better treatment of human conditions.
recent analyses of cellular events underlying post-traumatic regeneration in echinoderms  identified interesting parallels with corresponding processes in regeneration-competent vertebrates, pointing to reparative mechanisms that might have been evolutionary preserved throughout deuterostomia and could potentially be re-activated in poorly regenerating vertebrates.
however, the paucity of genomic and transcriptomic information has precluded major progress in understanding key regulators of regeneration at the molecular level  <cit> .
the body of an echinoderm has radial organization and is composed of multiple  units, which are arranged as sectors or rays around the central axis connecting the mouth and the anus.
each of these radial units is supplied with a set of major organs, which in sea cucumbers  include a radial nerve cord , a canal of the water-vascular system, and a longitudinal muscle band .
holothuria glaberrima
selenka,  <dig> .
our previous research  <cit>  showed that following a transverse cut, the injured organs of the radial organ complex on either side of the wound start growing across the wound gap and eventually reconnect to restore the anatomical continuity .
the newly regenerated structures then completely re-acquire their normal tissue architecture and resume their functions.
radial nerve cord regeneration involves extensive dedifferentiation of radial glial cells in the vicinity of the injury.
these dedifferentiated glial cells play the key role in subsequent regeneration through extensive proliferation, ecm invasion, and differentiation into new neurons and glial cells.
the newly produced neurons are thought to be functionally integrated into the cns circuitry, as they survive for extended periods of time and form typical synaptic connections  <cit> .
the main goal of the present study is to help understand the molecular basis underlying the extensive regenerative capacity of the central nervous system in the sea cucumber holothuria glaberrima by providing an outline of the transcriptomic landscape and thus identifying possible directions of future research.
to this end, we used deep rna sequencing on both the  <dig> and illumina platforms to analyze changes in the transcriptome that occurred on day  <dig>  day  <dig>  and day  <dig> after injury.
these time points were chosen based on our previous studies of cellular events in the regenerating radial nerve cord of h. glaberrima <cit> .
day  <dig> post-injury  is the early post-injury phase of extensive dedifferentiation of radial glia, axonal degeneration, and programmed cell death in the injured radial nerve.
day  <dig> post injury  corresponds to a period of active growth across the wound gap; dividing cells are most abundant at this stage.
day  <dig>  post injury is a late regeneration phase, when the two growing regenerates have restored their anatomical continuity and started to resume their typical histological architecture  <cit> .
the present study provides first insights into gene expression changes that underlie these previously described cellular events.
RESULTS
sequencing and assembly: technical information about the assembly
the present study was performed on the brown rock sea cucumber holothuria glaberrima selenka,  <dig> , an established model organism in echinoderm regenerative biology.
we injured the radial nerve cord and the surrounding tissues, including the body wall connective tissue, water-vascular canal, and the longitudinal muscle band , by performing a single transverse cut through these organs in the mid-ventral ambulacrum at about the mid-body level.
we then harvested the regenerating tissues on day  <dig>  day  <dig>  day  <dig>  and day  <dig> post-injury, as well as non-injured radial nerve cord and used the samples for sequencing library preparation.
the major steps involved in our sequencing and data analysis pipeline are outlined in figure  <dig>
the high-throughput transcriptome sequencing yielded a total of  <dig> , <dig> roche  <dig> reads with the average/modal length of 595/ <dig> bases and  <dig> , <dig> single-end illumina reads  .
we first combined the reads from all libraries to assemble a reference transcriptome for annotation purposes.
prior to the de-novo assembly step, we subjected the raw reads from both sequencing platforms to a rigorous cleaning/clipping procedure  leaving  <dig> , <dig> roche  <dig> reads  and  <dig> , <dig> illumina reads  for further analysis .
all cleaned  <dig> reads were assembled with mira  into  <dig>  contigs with the mean length of  <dig> bp and n <dig>  of  <dig> bp .
illumina reads were assembled with velvet and oases .
de bruijn graph assemblers, such as velvet/oases, are known to be sensitive to the value of parameter k, which is under the control of the user.
it has been shown that the use of a single k-mer length may result in suboptimal de novo transcriptome assemblies  <cit> .
for a given assembly, the best k-mer value depends on the coverage depth , the sequencing error rate, and the complexity of the transcriptome.
changing the k-mer value in either direction  is associated with both pros and cons.
using higher k-mer values increases specificity and leads to a better assembly of highly expressed transcripts, whereas lower k-mer values improves sensitivity and allows better reconstruction of weakly expressed transcripts  <cit> .
therefore, seven separate assembly runs were performed with different k-mer length .
they yielded between  <dig>  and  <dig>  contigs , with the mean length of  <dig> –  <dig> bp, the maximum length of  <dig>  –  <dig> , and n <dig> ranging between  <dig>  and  <dig>  .
then all mira and velvet contigs were pooled together and used as an input for the final assembly run with the cap <dig> program , which produced  <dig>  contigs with dramatically improved mean length  and n <dig>  and the maximum length of  <dig>  bp .
we refer to this final set of  <dig>  contigs  as a reference library, which represents the transcriptional diversity in the normal and regenerating radial organ complex of the sea cucumber h. glaberrima, and which we used for functional annotation and read mapping .
previously,  <dig> of these contigs were found to correspond to  <dig> long terminal repeat  transposons.
homology analysis and expression pattern of these mobile genetic elements were characterized in detail elsewhere  <cit>  and, therefore, will not be covered in this study.
to assess the accuracy of our sequencing and assembly pipeline, we re-sequenced  <dig> randomly selected contigs from the reference library using sanger technology .
in total,  <dig>  nt of  <dig>  nt  were analyzed by re-sequencing.
the overall error rate  was  <dig> %, with insertions, deletions, and substitutions affecting  <dig> nt,  <dig> nt, and  <dig> nt, respectively.
eleven of the resequenced transcripts had a complete open reading frame , in four contigs both the 5’ and 3’ ends of the orf were incomplete, in one contig, only the 5’ end was missing, and one contig was likely a non-coding rna.
one of the  <dig> orf-containing contigs had a frameshift error.
in order to determine the proportion of the contigs in the reference library, which corresponded to known proteins, we performed a blastx search against publicly available protein databases.
in the first round of search, the sequences were matched against the swissprot database.
the contigs, which did not produce hits passing the significance threshold corresponding to ane-value < 1e- <dig> were subsequently subjected to a second round of blastx search against the non-redundant  ncbi protein database with the same threshold.
overall,  <dig>   contigs had significant blast hits.
the results of the blast analysis are listed in additional file  <dig>
we also annotated the sea cucumber transcriptome by performing reciprocal best blast hit analysis  versus the ncbi’s collection of the sea urchin strongylocentrotus purpuratus predicted protein sequences  <cit> , the echinoderm species whose genome has been most thoroughly characterized so far.
there were  <dig>  one-way blastx matches between the h. glaberrima contigs and the s. purpuratus proteins and  <dig>  one-way tblastn matches between the sea urchin proteins and the contigs of our transcriptome.
the number of reciprocal blast matches  was  <dig> .
the results of this analysis are listed in additional file  <dig>
we, obviously, do not expect all our  <dig>  contigs to represent individual sea cucumber genes.
it is impossible to give an exact answer on the number of unique genes represented by de novo assembled contigs without having complete genomic information.
thus, we can only provide an educated guess here.
the genome of the sea urchin s. purpuratus <cit>  encodes ∼ <dig>  genes, of which ∼ <dig>  are presumed orthologs of mammalian genes.
in the course of our analysis, a similar number  of unique mouse proteins showed significant similarity to sequences contained in the assembled transcriptome of h. glaberrima.
if we assume the same level of genomic similarity between mammals and each of the two echinoderm species, the similar number of “mammalian genes” observed in the sea urchin genome and the sea cucumber transcriptome would imply that  sea cucumbers have roughly the same number of genes as sea urchins , and  most of these genes are represented in our reference library constructed from mrna samples from non-injured and regenerating animals.
if the above reasoning is correct, there is a ∼3× redundancy in our assembly.
a similar redundancy ratio  can be obtained, if we consider that, when blasted against the reference mouse proteome,  <dig>  of  <dig>  contigs of the sea cucumber reference library matched  <dig>  unique mouse proteins.
some of this redundancy is part of the “natural” variation that is expected from differences in mrna processing .
in fact, when manually inspecting assembled contigs, we saw polymorphic transcripts, which were especially common among retroelements.
additional variation can be due to the limitations of the sequencing techniques and/or de novo assembly programs.
differential gene expression in radial organ complex regeneration
to characterize differentially regulated genes involved in regeneration, we mapped the reads from each of the eight illumina libraries  to the contigs of the reference library .
gene expression values in each pair of replicates showed high correlation  .
we then used the deseq package  <cit>  to identify sets of significantly up- and down-regulated genes in the regenerating tissues as compared with the intact animals .
upon imposing a threshold adjusted p value  of less than  <dig>  and a threshold fold change of  <dig>  we identified a total of  <dig>  and  <dig>  transcripts, which were differentially up-regulated and down-regulated, respectively.
the distribution of unique and overlapping sets of the differentially expressed genes at different time points is summarized in figure 4b.
note that the early post-injury phase  is characterized by the highest values of both the total number of differentially expressed transcripts and the number of unique up-regulated and down-regulated genes.
on the other hand, in line with the general histological similarity between the late regenerates and uninjured radial organs   <cit> , the samples analyzed on day  <dig> post-injury showed the lowest number of differentially regulated transcripts.
at all three time points after the injury, the majority  of the significant differentially expressed contigs were up- or down-regulated between 2- and 4-fold relative to the normal animals , suggesting that relatively small changes in transcript abundance of most genes are sufficient for regeneration to occur.
among the most extreme outliers in our differential expression analyses were sequences that were identified as retrotransposon-derived transcripts.
some of them showed an over 50-fold change in expression during regeneration.
this unexpected finding prompted us to undertake a separate study of these mobile genetic elements that has already been published elsewhere  <cit> .
to verify the validity of our approach and the accuracy of the large-scale digital gene expression assay, we selected  <dig> genes for analysis by quantitative real-time rt-pcr .
this dataset also includes previously determined expression values for three ltr retroelements  characterized earlier elsewhere  <cit> .
our goal was to validate the sensitivity of the assay at different levels of transcript abundances.
therefore, we chose not only highly over-expressed or down-regulated genes, but also genes with moderate fold change and genes, whose transcript abundance remained stable throughout the experiment.
we found highly significant strong positive correlation between the rna-seq and qrt-pcr data  .
the corresponding fold change ratios  for both quantification techniques are listed in additional file  <dig>  this correlation analysis suggests that our rna-seq data reliably represent relative changes of mrna levels in the regenerating tissues.
functional annotation of differentially expressed genes at different time points of regeneration
from the biological perspective, it is important to know which processes predominate or become suppressed at different time points during regeneration.
for functional annotation, we employed an approach, which was conceptually similar to that described by stewart et al.
<cit>  for the analysis of the limb blastema transcriptome in axolotl.
in order to gain insight into the biological meaning behind the large numbers of differentially regulated transcripts, we analyzed the lists of genes that showed significant changes in expression levels  with david v <dig>   <cit>  to determine enriched functional categories.
to be able to use this tool, we had to match the contigs resulting from our de novo assembly of the sea cucumber transciptome to genes of some well-annotated model organism.
we thus annotated the assembled transcriptome of h. glaberrima against the non-redundant reference proteome of the mouse.
our assembled contigs matched to  <dig>  unique mouse genes at the blastx cut-off value of 1e- <dig>  this entire list of annotated contigs representing both uninjured and regenerating animals was submitted to david along with the lists of differentially expressed genes and used as a custom reference background.
functional annotation clustering was performed using gene ontology  and pathway  annotations.
the functional clustering algorithm in david groups together statistically enriched annotation categories if they share significant degrees of similarity in terms of associated gene members.
the ease score  threshold was set at  <dig> .
functional annotation clusters with enrichment scores  >  <dig>   were considered significant.
the following major findings emerged from this analysis .
<dig>  among the most enriched functional categories associated with both up-regulated and down-regulated genes at all time points were those related to synthesis and organization of the extracellular matrix  components, ecm remodeling, and interaction between cells and the extracellular matrix .
<dig>  functional annotation terms associated with normal physiology, differentiation, and development of the nervous tissue were over-represented among the down-regulated genes on days  <dig> and  <dig> post-injury, and some of those annotation categories were also enriched among negatively regulated genes even as late as after  <dig> days post-injury.
this observation correlates well with previously published morphological studies showing extensive dedifferentiation in the regenerating radial nerve at these time points  <cit> .
interestingly, annotation categories related to glycolysis were also enriched among the down-regulated genes on day  <dig> and  <dig>
since our tissue samples, in addition to the radial nerve cord per se, also contained some adjacent tissues, including the contractile epithelium of the water-vascular canal and the coelomic myoepithelium of the body wall , annotation terms associated with structural components of muscular tissue, as well as with muscle development and physiology, also appeared in the analysis.
these terms were enriched in sets of down-regulated genes at all three time points of regeneration corroborating earlier observations of muscular dedifferentiation during body wall regeneration in sea cucumbers  <cit> .
<dig>  on day  <dig> post-injury, the up-regulated genes were associated with over-represented terms related to initiation of dna replication and protein translation.
on day  <dig>  positively regulated genes were characterized by continued over-representation of terms associated with dna synthesis and cell cycle.
these data corroborate our cell proliferation assay, which suggested that the peak of cell division in the regenerating radial nerve cord occurs on days  <dig> thru  <dig> after injury  <cit> .
on day  <dig>  up-regulated genes are also enriched in annotation categories associated with activation of the innate immune response.
the latter group of functional terms remains over-represented among up-regulated genes on day  <dig>
<dig>  all three regeneration time points are characterized by over-representation of annotation categories associated with developmental morphogenesis, indicating involvement of ontogenic processes in regeneration.
in order to get further insight into mechanisms controlling gene expression in post-traumatic regeneration of the radial nerve cord, we performed clustering of differentially expressed genes using the unsupervised clustering algorithm implemented in the autosome program  <cit>  to identify sets of co-regulated genes with similar patterns of expression throughout the time course of regeneration.
this gene co-expression analysis yielded eight distinct clusters with more than  <dig> genes in each .
in two of these clusters , the median gene expression in regenerating animals was consistently higher at all three regeneration stages when compared with uninjured animals.
these genes were enriched in functional annotation terms related to the extracellular matrix, initiation of dna replication and translation, positive regulation of immune response, wound healing, and cell migration.
three other clusters  contained genes, whose median expression at all three time points was lower that in normal animals.
associated enriched functional annotation terms were related to extracellular matrix, transmission of nerve impulse, regulation of neurotransmitter activity, muscle cell and neuron development, cytoskeleton, mapk signaling pathway, negative regulation of cell death, gated channel activity, and carbohydrate binding.
clusters  <dig> and  <dig> contained genes, whose expression initially either remained unchanged or somewhat decreased at the early post-injury stage , but then was elevated during the growth  and late regeneration  phases.
these genes were associated with positive regulation of developmental processes, extracellular matrix, cell adhesion, serine peptidase activity, endocytosis.
cluster  <dig> contained genes, which were markedly down-regulated at the early post-injury stage , but whose expression than returned to the normal levels as regeneration progressed.
enriched functional annotation terms associated with this cluster included extracellular matrix, neuromuscular junction development, and synaptic transmission.
putative regulation of groups of co-expressed genes by transcription factors
we then sought to get insight into which upstream transcription factors  may regulate coordinated expression of genes within each of the eight clusters identified above .
over-representation of transcription factor binding sites associated with co-expressed genes of each cluster was predicted using opossum software  <cit> .
we considered only those of the predicted tfs whose significant  changes in expression level correlated in time  with changes in the median gene expression within the respective cluster.
this analysis yielded  <dig> tfs, which putatively regulate regeneration-related genes  and therefore constitute good candidates for further functional analysis.
it is worth noting that there was some overlap between the clusters.
thus, ctcf and nfkb <dig> were associated with four and three clusters, respectively, and srf, fli <dig>  and plag <dig> were associated with two different clusters each.
some of the tfs  were identified as functioning as both transcriptional activators and transcriptional repressors in different clusters.
among these predicted  <dig> tfs, there are genes that have been previously implicated in post-traumatic processes and developmental processes in various organisms.
for example, nfkb <dig> is known to affect expression of broad range of downstream genes involved in various biological processes including immunity, differentiation, and programmed cell death.
for example, elevated nfkb signaling has been previously shown to activate the wnt signaling pathway and thus induce dedifferentiation of nonstem cells  <cit> .
another tf, serum response factor , was shown to be one of the key genes involved in post-traumatic regeneration initiation in planarians  <cit> .
in gastric ulcer healing, srf promotes re-epithelialization and muscle regeneration through activation of cell migration and proliferation  <cit> .
moreover, srf is also implicated in neuronal cell migration and axonal guidance through regulation of components of the actin cytoskeleton  <cit> .
still another gene, ccaat/enhancer binding protein  was suggested to be involved in regulation of the innate immune response during tissue injury repair  <cit> .
interestingly, the candidate tfs identified in this study included not only previously known regeneration-related genes, but also factors not previously known to act as regulators of post-traumatic tissue regrowth, such as, for example, liver x receptor, fli <dig>  plag <dig>  ebf <dig>  esrrb.
the potential role of these genes in regeneration deserves further attention in future research.
discussion of selected functional gene groups
the preceding section provided an overview of unbiased functional characterization of the regeneration-associated genes at the global, transcriptome-wide level.
below, we further zoom in on certain groups of differentially expressed genes, which we picked from our database based on the results of the above analysis, prior biological knowledge, or both.
cancer-related pathways
as revealed by david analysis, the set of differentially expressed genes at the early post-injury stage  includes many known cancer-related genes .
for example, pathway mapping revealed differential regulation of wnt receptors  and ligands .
another notable observation is down-regulation of survivin  .
this gene codes for a multifunctional protein, highly expressed in most human cancers, and implicated in both suppression of programmed cell death  <cit>  and regulation of cell division  <cit> .
we have previously shown that elevated expression of survivin in regenerating sea cucumber intestinal tissues correlated with low levels of apoptosis  <cit> .
reduced expression of this gene in the regenerating radial nerve cord may thus be at least in part explained as being associated with extensive programmed cells death in the vicinity of the injury  <cit> .
differential expression of oncogenes early in regeneration has been also reported for other model systems, including early limb blastema in axolotl  <cit>  and can be explained by the fact that both regeneration and cancer progression are developmental processes that share a number of key mechanisms including cell division, programmed cell death, and differentiation.
ecm-related pathways
as mentioned above, among the most significantly over-represented pathways associated with differentially expressed genes at all three time points of the radial nerve cord regeneration were those related to the extracellular matrix  components and ecm-cell interactions .
one example is significantly reduced expression of genes coding for basal lamina proteins on day  <dig> , but up-regulation of the fibrillar collagens col5a <dig> and col11a <dig>  these data correlate well with previous observations of breakdown of the basal lamina during the early post-injury phase  <cit> .
of particular interest is concurrent up-regulation of both matrix metalloproteases  and their inhibitors, timps, in the regenerating radial organ complex.
by breaking down components of the connective tissue matrix at the injury site and thus affecting cell migration and proliferation, epithelialization, differentiation, and apoptosis, mmps are known to facilitate wound healing and regeneration  <cit> .
however, up-regulation of timps is as important for precise regulation of the regenerative processes, as they protect newly synthesized ecm matrix from degradation  <cit> .
proper interactions between cells and the ecm are essential both for correct tissue organization and signal transduction.
among the main cell adhesion molecules mediating these interactions are integrins  <cit> .
in regeneration of the radial nerve in the sea cucumber, integrin beta and integrin alpha- <dig> are significantly up-regulated .
integrins affect cytoskeletal organization through cdc <dig>  a member of the rho family of small gtpases  <cit> .
expression levels of the cdc <dig> transcript were consistently highly elevated in the sea cucumber radial nerve regeneration, suggesting a possible involvement of the integrin-cdc <dig> pathway in control of ecm-cell interactions in this animal model.
besides integrins, other cell adhesion molecules are also up-regulated during the early post-injury phase.
these include selectins .
interestingly, expression of selectins is known to be induced by nf-kb  <cit> , which was putatively identified as one of the putative key transcription factors driving differential gene expression during the early post-injury stage of neural regeneration in our model .
pluripotency factors
as has been documented by morphological studies  <cit> , regeneration of the radial organs in sea cucumbers involves extensive dedifferentiation of specialized cells of the adult tissues in the vicinity of the injury.
the dedifferentiation step endows these cells with increased proliferative potential, as well as with the ability to migrate and subsequently re-acquire their initial morphology or, more interestingly, transdifferentiate to other cell types.
as this dedifferentiation phase is an indispensable component of the regenerative response, we asked if it requires elevated expression of factors known to be involved in acquisition and/or control of pluripotency.
we, therefore, searched our assembled transcriptome for homologs of those genes and analyzed their expression values based on our rna-seq read count data .
out of  <dig> factors analyzed,  <dig> were found to be expressed in both normal and regenerating tissues.
of those  <dig> genes, only the myc homolog was expressed at significantly higher levels in regeneration.
myc genes have a wide variety of functions affecting the developmental potential of cells, including driving cells into the cell cycle, promoting the open state of chromatin, and dna replication.
c-myc is one of the yamanaka factors involved in the production of induced pluripotent cells  <cit>  from specialized mammalian cells.
during this reprogramming procedure, myc represses differentiation-associated genes  <cit> .
the functional role of myc in echinoderm regeneration is a subject of ongoing research in our group.
to our surprise, we also found that expression of the bmi- <dig> homolog was significantly reduced in the regenerating tissues on days  <dig> and  <dig> post-injury.
the remaining  <dig> genes were expressed at the same level both in the non-injured and regenerating tissues.
a possible explanation that can be proposed to account for the observed data is that although the analyzed pluripotency factors do not show any large scale over-expression in regeneration, most of them are still expressed at a certain level both in the normal and regenerating tissues.
therefore, it may be hypothesized that this basal level constitutes sufficient environment for dedifferentiation to be triggered by myc upregulation alone.
interestingly, a similar situation has been reported in non-injured and regenerating tissues of lower vertebrates, such as danio rerio and xenopus, where pluripotency markers are never shut off completely under normal conditions.
even though the expression level of these genes remained largely at the same level after the injury, it was hypothesized that this basal expression is neverhteless sufficient to facilitate regeneration upon damage  <cit> .
genes associated with neurogenesis
post-traumatic neurogenesis in the lesioned radial nerve cord involves extensive proliferation of radial glial cells, which then act as progenitors generating both new glial cells and neurons  <cit> .
surprisingly, our functional analysis of differentially expressed genes did not yield statistically significant enrichment of neurogenesis-associated functional categories.
the only exception was the enrichment of the category “neuron development” in the set of down-regulated genes during the growth phase on day  <dig> post-injury .
we then manually searched our databases for homologs of genes known to be implicated in maintenance of neural stem cells and in specification of the neuronal and glial lineages in other animal models  .
the results of this analysis are summarized in figure  <dig>  the most probable explanation why neurogenesis was not a dominant theme in the functional analysis of differentially expressed genes is that most of the neurogenesis-related genes, whose homologs were identified in the transcriptome of h. glaberrima, were already expressed in the non-injured radial nerve cord.
their expression level showed some fluctuations in the regenerating tissues, but not above the cut-off threshold , with the exception of up-regulated sox <dig>  bmp <dig> and down-regulated prom <dig>  isl <dig>  and lhx <dig>
to our surprise, homologs of some of the key markers of vertebrate neural stem cells, including nestin and vimentin, were absent from the sea cucumber transcriptome.
likewise, transcripts of some of the important pro-neural factors  <cit> , such as mash <dig> and neurog <dig> were also not detected in either the normal or regenerating animals.
this may be due to the fact that the program of post-traumatic neurogenesis in sea cucumbers is not known and, obviously, may not entirely consist of the mechanisms that regulate neurogenesis in vertebrates, whose genes were used as the reference.
prevention of excitotoxic neuronal cell death
one of the adaptations that may facilitate neural regeneration in the sea cucumber h. glaberrima is inhibition of excitotoxicity .
in vertebrates, nervous system trauma results in two types of degenerative processes: the primary damage as the direct result of the injury, which affects the cells in the immediate vicinity of the trauma, and the secondary degeneration, which significantly increases the amount of the affected tissue and is mediated by endogenous mechanisms triggered by the initial injury  <cit> .
excitotoxicity is one of these secondary pathological processes.
it involves neuronal death via excessive stimulation by glutamate, whose extracellular concentration increases after cns injury.
our results suggest that transection of the radial organ complex in the sea cucumber h. glaberrima results in significant down-regulation of glutamate carboxypeptidase ii  in the vicinity of the injury .
this enzyme hydrolyzes naag  to glutamate and n-acetylaspartate.
it can be hypothesized that down-regulation of gcpii will eventually result in inhibition of naag hydrolysis and accumulation of intact naag, which is known to further decrease the amount of glutamate released from presynaptic terminals and therefore prevent excitotoxic neuronal cell death  <cit> .
moreover, naag has been reported to stimulate production of important regulators of neuron survival, such as tgf- β <cit> .
a sequence strongly related to tgf- β is present in the sea cucumber transcriptome.
it shows a trend toward up-regulation on day  <dig> and is significantly up-regulated above the cut-off threshold on days  <dig> and  <dig> .
on the other hand, the transcripts coding for another enzyme involved in glutamate production  and for the glutamate receptors  are also significantly down-regulated.
suppression of pathways leading to induction of excitotoxic neuronal death and concomitant induction of factors involved in neuronal survival can represent one of the adaptations aimed at minimizing the traumatic impact of the neuronal injury and sparing the maximum possible number of neurons, so that less of them will need to be replenished in the course of subsequent regeneration.
a promising further line of research would be to experimentally test if changes in the severity of excitotoxicity would actually affect the extent of tissue damage and the rate of re-connection in the injured radial nerve cord.
CONCLUSIONS
this study was conceived as a first stage in exploring molecular mechanisms behind the observed cellular processes in echinoderm cns regeneration.
in general, functional annotations of the differentially expressed genes corroborate well our previous morphological data, but also open up potential new avenues for future research.
as mentioned above, our results point out to the important role of the ecm remodeling in regeneration of the radial complex in the sea cucumber.
so far, reorganization of the connective tissue in regenerating echinoderms has received little attention.
there have been just two experimental studies addressing this issue directly and in both cases they were focused on regeneration of the digestive tube only  <cit> .
another promising line of future research would be to experimentally test the predictions of this study suggesting the existence of mechanisms suppressing excitotoxicity in the injured cns.
if these mechanisms actually exist, their understanding will be valuable for devising new therapies to prevent excitocytotic neuronal death following, for example, brain and spinal cord injury.
one of the most important outcomes from this study is a predicted list of putative transcription factors, which presumably control differential expression of large groups of downstream genes and thus occupy key positions in regulatory networks controlling regeneration.
these genes represent promising candidates for future functional analysis.
among other interesting findings is that post-traumatic regeneration of the radial nerve cord did not involve large-scale over-expression of pluripotency factors.
many of these genes were already expressed in intact tissues, and only myc showed up-regulation after injury.
we are well aware of the limitations in our study.
for example, by design, our functional annotation was dependent on matching the sea cucumber contigs against a well studied proteome from another organism .
this approach could have led to many potentially relevant sea cucumber-specific sequences being excluded from the analysis.
this issue cannot be resolved without carrying out a separate study aimed at characterization of the ‘new’ or ‘unknown’ sequences, which do not have significant homologs in current databases.
the present paper, like most of other high-throughput studies of gene expression, uses mrna abundance levels to get insight into how changes in gene expression might affect the phenotype of tissues and cells, although it is largely the quantity of the protein that directly determines the phenotype.
however, for the time being, transcriptomic approaches are justified by the fact that the current methodologies for direct quantification of protein expression are either less reliable or more laborious and expensive than mrna-based studies.
notwithstanding the limitations, the results and predictions reported above are valuable, because they provide a number of clearly defined testable hypotheses, whereas the associated pitfalls and limitations are well known to many researchers working with ‘non-model’ organisms and are not unresolvable in the future.

 echinoderms constitute a phylum of marine invertebrates closely related to chordates. their phylogenetic position as a non-chordate deuterostome phylum combined with the ability to regenerate various body parts makes them a valuable group, which can give unique insights into fundamental issues of regenerative biology and provide new clues to finding better treatment of human conditions. recent analyses of cellular events underlying post-traumatic regeneration in echinoderms identified interesting parallels with corresponding processes in regeneration-competent vertebrates, pointing to reparative mechanisms that might have been evolutionary preserved throughout deuterostomia and could potentially be re-activated in poorly regenerating vertebrates. however, the paucity of genomic and transcriptomic information has precluded major progress in understanding key regulators of regeneration at the molecular level. the body of an echinoderm has radial organization and is composed of multiple units, which are arranged as sectors or rays around the central axis connecting the mouth and the anus. each of these radial units is supplied with a set of major organs, which in sea cucumbers include a radial nerve cord, a canal of the water-vascular system, and a longitudinal muscle band. holothuria glaberrima selenka,. our previous research showed that following a transverse cut, the injured organs of the radial organ complex on either side of the wound start growing across the wound gap and eventually reconnect to restore the anatomical continuity. the newly regenerated structures then completely re-acquire their normal tissue architecture and resume their functions. radial nerve cord regeneration involves extensive dedifferentiation of radial glial cells in the vicinity of the injury. these dedifferentiated glial cells play the key role in subsequent regeneration through extensive proliferation, ecm invasion, and differentiation into new neurons and glial cells. the newly produced neurons are thought to be functionally integrated into the cns circuitry, as they survive for extended periods of time and form typical synaptic connections. the main goal of the present study is to help understand the molecular basis underlying the extensive regenerative capacity of the central nervous system in the sea cucumber holothuria glaberrima by providing an outline of the transcriptomic landscape and thus identifying possible directions of future research. to this end, we used deep rna sequencing on both the and illumina platforms to analyze changes in the transcriptome that occurred on day day and day after injury. these time points were chosen based on our previous studies of cellular events in the regenerating radial nerve cord of h. glaberrima. day post-injury is the early post-injury phase of extensive dedifferentiation of radial glia, axonal degeneration, and programmed cell death in the injured radial nerve. day post injury corresponds to a period of active growth across the wound gap; dividing cells are most abundant at this stage. day post injury is a late regeneration phase, when the two growing regenerates have restored their anatomical continuity and started to resume their typical histological architecture. the present study provides first insights into gene expression changes that underlie these previously described cellular events.  sequencing and assembly: technical information about the assembly the present study was performed on the brown rock sea cucumber holothuria glaberrima selenka,, an established model organism in echinoderm regenerative biology. we injured the radial nerve cord and the surrounding tissues, including the body wall connective tissue, water-vascular canal, and the longitudinal muscle band, by performing a single transverse cut through these organs in the mid-ventral ambulacrum at about the mid-body level. we then harvested the regenerating tissues on day day day and day post-injury, as well as non-injured radial nerve cord and used the samples for sequencing library preparation. the major steps involved in our sequencing and data analysis pipeline are outlined in figure  the high-throughput transcriptome sequencing yielded a total of, roche reads with the average/modal length of 595/ bases and, single-end illumina reads. we first combined the reads from all libraries to assemble a reference transcriptome for annotation purposes. prior to the de-novo assembly step, we subjected the raw reads from both sequencing platforms to a rigorous cleaning/clipping procedure leaving, roche reads and, illumina reads for further analysis. all cleaned reads were assembled with mira into contigs with the mean length of bp and n of bp. illumina reads were assembled with velvet and oases. de bruijn graph assemblers, such as velvet/oases, are known to be sensitive to the value of parameter k, which is under the control of the user. it has been shown that the use of a single k-mer length may result in suboptimal de novo transcriptome assemblies. for a given assembly, the best k-mer value depends on the coverage depth, the sequencing error rate, and the complexity of the transcriptome. changing the k-mer value in either direction is associated with both pros and cons. using higher k-mer values increases specificity and leads to a better assembly of highly expressed transcripts, whereas lower k-mer values improves sensitivity and allows better reconstruction of weakly expressed transcripts. therefore, seven separate assembly runs were performed with different k-mer length. they yielded between and contigs, with the mean length of – bp, the maximum length of –, and n ranging between and. then all mira and velvet contigs were pooled together and used as an input for the final assembly run with the cap program, which produced contigs with dramatically improved mean length and n and the maximum length of bp. we refer to this final set of contigs as a reference library, which represents the transcriptional diversity in the normal and regenerating radial organ complex of the sea cucumber h. glaberrima, and which we used for functional annotation and read mapping. previously, of these contigs were found to correspond to long terminal repeat transposons. homology analysis and expression pattern of these mobile genetic elements were characterized in detail elsewhere and, therefore, will not be covered in this study. to assess the accuracy of our sequencing and assembly pipeline, we re-sequenced randomly selected contigs from the reference library using sanger technology. in total, nt of nt were analyzed by re-sequencing. the overall error rate was %, with insertions, deletions, and substitutions affecting nt, nt, and nt, respectively. eleven of the resequenced transcripts had a complete open reading frame, in four contigs both the 5’ and 3’ ends of the orf were incomplete, in one contig, only the 5’ end was missing, and one contig was likely a non-coding rna. one of the orf-containing contigs had a frameshift error. in order to determine the proportion of the contigs in the reference library, which corresponded to known proteins, we performed a blastx search against publicly available protein databases. in the first round of search, the sequences were matched against the swissprot database. the contigs, which did not produce hits passing the significance threshold corresponding to ane-value < 1e- were subsequently subjected to a second round of blastx search against the non-redundant ncbi protein database with the same threshold. overall, contigs had significant blast hits. the results of the blast analysis are listed in additional file we also annotated the sea cucumber transcriptome by performing reciprocal best blast hit analysis versus the ncbi’s collection of the sea urchin strongylocentrotus purpuratus predicted protein sequences, the echinoderm species whose genome has been most thoroughly characterized so far. there were one-way blastx matches between the h. glaberrima contigs and the s. purpuratus proteins and one-way tblastn matches between the sea urchin proteins and the contigs of our transcriptome. the number of reciprocal blast matches was. the results of this analysis are listed in additional file we, obviously, do not expect all our contigs to represent individual sea cucumber genes. it is impossible to give an exact answer on the number of unique genes represented by de novo assembled contigs without having complete genomic information. thus, we can only provide an educated guess here. the genome of the sea urchin s. purpuratus encodes ∼ genes, of which ∼ are presumed orthologs of mammalian genes. in the course of our analysis, a similar number of unique mouse proteins showed significant similarity to sequences contained in the assembled transcriptome of h. glaberrima. if we assume the same level of genomic similarity between mammals and each of the two echinoderm species, the similar number of “mammalian genes” observed in the sea urchin genome and the sea cucumber transcriptome would imply that sea cucumbers have roughly the same number of genes as sea urchins, and most of these genes are represented in our reference library constructed from mrna samples from non-injured and regenerating animals. if the above reasoning is correct, there is a ∼3× redundancy in our assembly. a similar redundancy ratio can be obtained, if we consider that, when blasted against the reference mouse proteome, of contigs of the sea cucumber reference library matched unique mouse proteins. some of this redundancy is part of the “natural” variation that is expected from differences in mrna processing. in fact, when manually inspecting assembled contigs, we saw polymorphic transcripts, which were especially common among retroelements. additional variation can be due to the limitations of the sequencing techniques and/or de novo assembly programs. differential gene expression in radial organ complex regeneration to characterize differentially regulated genes involved in regeneration, we mapped the reads from each of the eight illumina libraries to the contigs of the reference library. gene expression values in each pair of replicates showed high correlation. we then used the deseq package to identify sets of significantly up- and down-regulated genes in the regenerating tissues as compared with the intact animals. upon imposing a threshold adjusted p value of less than and a threshold fold change of we identified a total of and transcripts, which were differentially up-regulated and down-regulated, respectively. the distribution of unique and overlapping sets of the differentially expressed genes at different time points is summarized in figure 4b. note that the early post-injury phase is characterized by the highest values of both the total number of differentially expressed transcripts and the number of unique up-regulated and down-regulated genes. on the other hand, in line with the general histological similarity between the late regenerates and uninjured radial organs, the samples analyzed on day post-injury showed the lowest number of differentially regulated transcripts. at all three time points after the injury, the majority of the significant differentially expressed contigs were up- or down-regulated between 2- and 4-fold relative to the normal animals, suggesting that relatively small changes in transcript abundance of most genes are sufficient for regeneration to occur. among the most extreme outliers in our differential expression analyses were sequences that were identified as retrotransposon-derived transcripts. some of them showed an over 50-fold change in expression during regeneration. this unexpected finding prompted us to undertake a separate study of these mobile genetic elements that has already been published elsewhere. to verify the validity of our approach and the accuracy of the large-scale digital gene expression assay, we selected genes for analysis by quantitative real-time rt-pcr. this dataset also includes previously determined expression values for three ltr retroelements characterized earlier elsewhere. our goal was to validate the sensitivity of the assay at different levels of transcript abundances. therefore, we chose not only highly over-expressed or down-regulated genes, but also genes with moderate fold change and genes, whose transcript abundance remained stable throughout the experiment. we found highly significant strong positive correlation between the rna-seq and qrt-pcr data. the corresponding fold change ratios for both quantification techniques are listed in additional file this correlation analysis suggests that our rna-seq data reliably represent relative changes of mrna levels in the regenerating tissues. functional annotation of differentially expressed genes at different time points of regeneration from the biological perspective, it is important to know which processes predominate or become suppressed at different time points during regeneration. for functional annotation, we employed an approach, which was conceptually similar to that described by stewart et al. for the analysis of the limb blastema transcriptome in axolotl. in order to gain insight into the biological meaning behind the large numbers of differentially regulated transcripts, we analyzed the lists of genes that showed significant changes in expression levels with david v to determine enriched functional categories. to be able to use this tool, we had to match the contigs resulting from our de novo assembly of the sea cucumber transciptome to genes of some well-annotated model organism. we thus annotated the assembled transcriptome of h. glaberrima against the non-redundant reference proteome of the mouse. our assembled contigs matched to unique mouse genes at the blastx cut-off value of 1e- this entire list of annotated contigs representing both uninjured and regenerating animals was submitted to david along with the lists of differentially expressed genes and used as a custom reference background. functional annotation clustering was performed using gene ontology and pathway annotations. the functional clustering algorithm in david groups together statistically enriched annotation categories if they share significant degrees of similarity in terms of associated gene members. the ease score threshold was set at. functional annotation clusters with enrichment scores > were considered significant. the following major findings emerged from this analysis. among the most enriched functional categories associated with both up-regulated and down-regulated genes at all time points were those related to synthesis and organization of the extracellular matrix components, ecm remodeling, and interaction between cells and the extracellular matrix. functional annotation terms associated with normal physiology, differentiation, and development of the nervous tissue were over-represented among the down-regulated genes on days and post-injury, and some of those annotation categories were also enriched among negatively regulated genes even as late as after days post-injury. this observation correlates well with previously published morphological studies showing extensive dedifferentiation in the regenerating radial nerve at these time points. interestingly, annotation categories related to glycolysis were also enriched among the down-regulated genes on day and since our tissue samples, in addition to the radial nerve cord per se, also contained some adjacent tissues, including the contractile epithelium of the water-vascular canal and the coelomic myoepithelium of the body wall, annotation terms associated with structural components of muscular tissue, as well as with muscle development and physiology, also appeared in the analysis. these terms were enriched in sets of down-regulated genes at all three time points of regeneration corroborating earlier observations of muscular dedifferentiation during body wall regeneration in sea cucumbers. on day post-injury, the up-regulated genes were associated with over-represented terms related to initiation of dna replication and protein translation. on day positively regulated genes were characterized by continued over-representation of terms associated with dna synthesis and cell cycle. these data corroborate our cell proliferation assay, which suggested that the peak of cell division in the regenerating radial nerve cord occurs on days thru after injury. on day up-regulated genes are also enriched in annotation categories associated with activation of the innate immune response. the latter group of functional terms remains over-represented among up-regulated genes on day all three regeneration time points are characterized by over-representation of annotation categories associated with developmental morphogenesis, indicating involvement of ontogenic processes in regeneration. in order to get further insight into mechanisms controlling gene expression in post-traumatic regeneration of the radial nerve cord, we performed clustering of differentially expressed genes using the unsupervised clustering algorithm implemented in the autosome program to identify sets of co-regulated genes with similar patterns of expression throughout the time course of regeneration. this gene co-expression analysis yielded eight distinct clusters with more than genes in each. in two of these clusters, the median gene expression in regenerating animals was consistently higher at all three regeneration stages when compared with uninjured animals. these genes were enriched in functional annotation terms related to the extracellular matrix, initiation of dna replication and translation, positive regulation of immune response, wound healing, and cell migration. three other clusters contained genes, whose median expression at all three time points was lower that in normal animals. associated enriched functional annotation terms were related to extracellular matrix, transmission of nerve impulse, regulation of neurotransmitter activity, muscle cell and neuron development, cytoskeleton, mapk signaling pathway, negative regulation of cell death, gated channel activity, and carbohydrate binding. clusters and contained genes, whose expression initially either remained unchanged or somewhat decreased at the early post-injury stage, but then was elevated during the growth and late regeneration phases. these genes were associated with positive regulation of developmental processes, extracellular matrix, cell adhesion, serine peptidase activity, endocytosis. cluster contained genes, which were markedly down-regulated at the early post-injury stage, but whose expression than returned to the normal levels as regeneration progressed. enriched functional annotation terms associated with this cluster included extracellular matrix, neuromuscular junction development, and synaptic transmission. putative regulation of groups of co-expressed genes by transcription factors we then sought to get insight into which upstream transcription factors may regulate coordinated expression of genes within each of the eight clusters identified above. over-representation of transcription factor binding sites associated with co-expressed genes of each cluster was predicted using opossum software. we considered only those of the predicted tfs whose significant changes in expression level correlated in time with changes in the median gene expression within the respective cluster. this analysis yielded tfs, which putatively regulate regeneration-related genes and therefore constitute good candidates for further functional analysis. it is worth noting that there was some overlap between the clusters. thus, ctcf and nfkb were associated with four and three clusters, respectively, and srf, fli and plag were associated with two different clusters each. some of the tfs were identified as functioning as both transcriptional activators and transcriptional repressors in different clusters. among these predicted tfs, there are genes that have been previously implicated in post-traumatic processes and developmental processes in various organisms. for example, nfkb is known to affect expression of broad range of downstream genes involved in various biological processes including immunity, differentiation, and programmed cell death. for example, elevated nfkb signaling has been previously shown to activate the wnt signaling pathway and thus induce dedifferentiation of nonstem cells. another tf, serum response factor, was shown to be one of the key genes involved in post-traumatic regeneration initiation in planarians. in gastric ulcer healing, srf promotes re-epithelialization and muscle regeneration through activation of cell migration and proliferation. moreover, srf is also implicated in neuronal cell migration and axonal guidance through regulation of components of the actin cytoskeleton. still another gene, ccaat/enhancer binding protein was suggested to be involved in regulation of the innate immune response during tissue injury repair. interestingly, the candidate tfs identified in this study included not only previously known regeneration-related genes, but also factors not previously known to act as regulators of post-traumatic tissue regrowth, such as, for example, liver x receptor, fli plag ebf esrrb. the potential role of these genes in regeneration deserves further attention in future research. discussion of selected functional gene groups the preceding section provided an overview of unbiased functional characterization of the regeneration-associated genes at the global, transcriptome-wide level. below, we further zoom in on certain groups of differentially expressed genes, which we picked from our database based on the results of the above analysis, prior biological knowledge, or both. cancer-related pathways as revealed by david analysis, the set of differentially expressed genes at the early post-injury stage includes many known cancer-related genes. for example, pathway mapping revealed differential regulation of wnt receptors and ligands. another notable observation is down-regulation of survivin. this gene codes for a multifunctional protein, highly expressed in most human cancers, and implicated in both suppression of programmed cell death and regulation of cell division. we have previously shown that elevated expression of survivin in regenerating sea cucumber intestinal tissues correlated with low levels of apoptosis. reduced expression of this gene in the regenerating radial nerve cord may thus be at least in part explained as being associated with extensive programmed cells death in the vicinity of the injury. differential expression of oncogenes early in regeneration has been also reported for other model systems, including early limb blastema in axolotl and can be explained by the fact that both regeneration and cancer progression are developmental processes that share a number of key mechanisms including cell division, programmed cell death, and differentiation. ecm-related pathways as mentioned above, among the most significantly over-represented pathways associated with differentially expressed genes at all three time points of the radial nerve cord regeneration were those related to the extracellular matrix components and ecm-cell interactions. one example is significantly reduced expression of genes coding for basal lamina proteins on day, but up-regulation of the fibrillar collagens col5a and col11a these data correlate well with previous observations of breakdown of the basal lamina during the early post-injury phase. of particular interest is concurrent up-regulation of both matrix metalloproteases and their inhibitors, timps, in the regenerating radial organ complex. by breaking down components of the connective tissue matrix at the injury site and thus affecting cell migration and proliferation, epithelialization, differentiation, and apoptosis, mmps are known to facilitate wound healing and regeneration. however, up-regulation of timps is as important for precise regulation of the regenerative processes, as they protect newly synthesized ecm matrix from degradation. proper interactions between cells and the ecm are essential both for correct tissue organization and signal transduction. among the main cell adhesion molecules mediating these interactions are integrins. in regeneration of the radial nerve in the sea cucumber, integrin beta and integrin alpha- are significantly up-regulated. integrins affect cytoskeletal organization through cdc a member of the rho family of small gtpases. expression levels of the cdc transcript were consistently highly elevated in the sea cucumber radial nerve regeneration, suggesting a possible involvement of the integrin-cdc pathway in control of ecm-cell interactions in this animal model. besides integrins, other cell adhesion molecules are also up-regulated during the early post-injury phase. these include selectins. interestingly, expression of selectins is known to be induced by nf-kb, which was putatively identified as one of the putative key transcription factors driving differential gene expression during the early post-injury stage of neural regeneration in our model. pluripotency factors as has been documented by morphological studies, regeneration of the radial organs in sea cucumbers involves extensive dedifferentiation of specialized cells of the adult tissues in the vicinity of the injury. the dedifferentiation step endows these cells with increased proliferative potential, as well as with the ability to migrate and subsequently re-acquire their initial morphology or, more interestingly, transdifferentiate to other cell types. as this dedifferentiation phase is an indispensable component of the regenerative response, we asked if it requires elevated expression of factors known to be involved in acquisition and/or control of pluripotency. we, therefore, searched our assembled transcriptome for homologs of those genes and analyzed their expression values based on our rna-seq read count data. out of factors analyzed, were found to be expressed in both normal and regenerating tissues. of those genes, only the myc homolog was expressed at significantly higher levels in regeneration. myc genes have a wide variety of functions affecting the developmental potential of cells, including driving cells into the cell cycle, promoting the open state of chromatin, and dna replication. c-myc is one of the yamanaka factors involved in the production of induced pluripotent cells from specialized mammalian cells. during this reprogramming procedure, myc represses differentiation-associated genes. the functional role of myc in echinoderm regeneration is a subject of ongoing research in our group. to our surprise, we also found that expression of the bmi- homolog was significantly reduced in the regenerating tissues on days and post-injury. the remaining genes were expressed at the same level both in the non-injured and regenerating tissues. a possible explanation that can be proposed to account for the observed data is that although the analyzed pluripotency factors do not show any large scale over-expression in regeneration, most of them are still expressed at a certain level both in the normal and regenerating tissues. therefore, it may be hypothesized that this basal level constitutes sufficient environment for dedifferentiation to be triggered by myc upregulation alone. interestingly, a similar situation has been reported in non-injured and regenerating tissues of lower vertebrates, such as danio rerio and xenopus, where pluripotency markers are never shut off completely under normal conditions. even though the expression level of these genes remained largely at the same level after the injury, it was hypothesized that this basal expression is neverhteless sufficient to facilitate regeneration upon damage. genes associated with neurogenesis post-traumatic neurogenesis in the lesioned radial nerve cord involves extensive proliferation of radial glial cells, which then act as progenitors generating both new glial cells and neurons. surprisingly, our functional analysis of differentially expressed genes did not yield statistically significant enrichment of neurogenesis-associated functional categories. the only exception was the enrichment of the category “neuron development” in the set of down-regulated genes during the growth phase on day post-injury. we then manually searched our databases for homologs of genes known to be implicated in maintenance of neural stem cells and in specification of the neuronal and glial lineages in other animal models. the results of this analysis are summarized in figure  the most probable explanation why neurogenesis was not a dominant theme in the functional analysis of differentially expressed genes is that most of the neurogenesis-related genes, whose homologs were identified in the transcriptome of h. glaberrima, were already expressed in the non-injured radial nerve cord. their expression level showed some fluctuations in the regenerating tissues, but not above the cut-off threshold, with the exception of up-regulated sox bmp and down-regulated prom isl and lhx to our surprise, homologs of some of the key markers of vertebrate neural stem cells, including nestin and vimentin, were absent from the sea cucumber transcriptome. likewise, transcripts of some of the important pro-neural factors, such as mash and neurog were also not detected in either the normal or regenerating animals. this may be due to the fact that the program of post-traumatic neurogenesis in sea cucumbers is not known and, obviously, may not entirely consist of the mechanisms that regulate neurogenesis in vertebrates, whose genes were used as the reference. prevention of excitotoxic neuronal cell death one of the adaptations that may facilitate neural regeneration in the sea cucumber h. glaberrima is inhibition of excitotoxicity. in vertebrates, nervous system trauma results in two types of degenerative processes: the primary damage as the direct result of the injury, which affects the cells in the immediate vicinity of the trauma, and the secondary degeneration, which significantly increases the amount of the affected tissue and is mediated by endogenous mechanisms triggered by the initial injury. excitotoxicity is one of these secondary pathological processes. it involves neuronal death via excessive stimulation by glutamate, whose extracellular concentration increases after cns injury. our results suggest that transection of the radial organ complex in the sea cucumber h. glaberrima results in significant down-regulation of glutamate carboxypeptidase ii in the vicinity of the injury. this enzyme hydrolyzes naag to glutamate and n-acetylaspartate. it can be hypothesized that down-regulation of gcpii will eventually result in inhibition of naag hydrolysis and accumulation of intact naag, which is known to further decrease the amount of glutamate released from presynaptic terminals and therefore prevent excitotoxic neuronal cell death. moreover, naag has been reported to stimulate production of important regulators of neuron survival, such as tgf- β. a sequence strongly related to tgf- β is present in the sea cucumber transcriptome. it shows a trend toward up-regulation on day and is significantly up-regulated above the cut-off threshold on days and. on the other hand, the transcripts coding for another enzyme involved in glutamate production and for the glutamate receptors are also significantly down-regulated. suppression of pathways leading to induction of excitotoxic neuronal death and concomitant induction of factors involved in neuronal survival can represent one of the adaptations aimed at minimizing the traumatic impact of the neuronal injury and sparing the maximum possible number of neurons, so that less of them will need to be replenished in the course of subsequent regeneration. a promising further line of research would be to experimentally test if changes in the severity of excitotoxicity would actually affect the extent of tissue damage and the rate of re-connection in the injured radial nerve cord.  this study was conceived as a first stage in exploring molecular mechanisms behind the observed cellular processes in echinoderm cns regeneration. in general, functional annotations of the differentially expressed genes corroborate well our previous morphological data, but also open up potential new avenues for future research. as mentioned above, our results point out to the important role of the ecm remodeling in regeneration of the radial complex in the sea cucumber. so far, reorganization of the connective tissue in regenerating echinoderms has received little attention. there have been just two experimental studies addressing this issue directly and in both cases they were focused on regeneration of the digestive tube only. another promising line of future research would be to experimentally test the predictions of this study suggesting the existence of mechanisms suppressing excitotoxicity in the injured cns. if these mechanisms actually exist, their understanding will be valuable for devising new therapies to prevent excitocytotic neuronal death following, for example, brain and spinal cord injury. one of the most important outcomes from this study is a predicted list of putative transcription factors, which presumably control differential expression of large groups of downstream genes and thus occupy key positions in regulatory networks controlling regeneration. these genes represent promising candidates for future functional analysis. among other interesting findings is that post-traumatic regeneration of the radial nerve cord did not involve large-scale over-expression of pluripotency factors. many of these genes were already expressed in intact tissues, and only myc showed up-regulation after injury. we are well aware of the limitations in our study. for example, by design, our functional annotation was dependent on matching the sea cucumber contigs against a well studied proteome from another organism. this approach could have led to many potentially relevant sea cucumber-specific sequences being excluded from the analysis. this issue cannot be resolved without carrying out a separate study aimed at characterization of the ‘new’ or ‘unknown’ sequences, which do not have significant homologs in current databases. the present paper, like most of other high-throughput studies of gene expression, uses mrna abundance levels to get insight into how changes in gene expression might affect the phenotype of tissues and cells, although it is largely the quantity of the protein that directly determines the phenotype. however, for the time being, transcriptomic approaches are justified by the fact that the current methodologies for direct quantification of protein expression are either less reliable or more laborious and expensive than mrna-based studies. notwithstanding the limitations, the results and predictions reported above are valuable, because they provide a number of clearly defined testable hypotheses, whereas the associated pitfalls and limitations are well known to many researchers working with ‘non-model’ organisms and are not unresolvable in the future. 
we also searched the sea cucumber transcriptome for homologs of factors known to be involved in acquisition and/or control of pluripotency. the identified candidate transcription factors included not only known regeneration-related genes, but also factors not previously implicated as regulators of post-traumatic tissue regrowth. our analysis yielded eleven putative transcription factors, which constitute good candidates for further functional studies. our de novo assembly yielded contigs, of which showed significant similarity to known protein-coding sequences. functional analysis of sets of differentially expressed genes suggested that among the most extensively over-represented pathways were those involved in the extracellular matrix remodeling and ecm-cell interactions, indicating a key role of the ecm in regeneration. more importantly, however, it also illuminates new aspects of echinoderm regeneration, which have been scarcely studied or overlooked altogether. in this study, we employ high-throughput sequencing to analyze the transcriptome of the normal and regenerating radial nerve cord, in the sea cucumber holothuria glaberrima. functional annotation also suggested that one of the possible adaptations contributing to fast and efficient neural regeneration in echinoderms may be related to suppression of excitotoxicity. echinoderms are emerging as important models in regenerative biology. significant amount of data are available on cellular mechanisms of post-traumatic repair in these animals, whereas studies of gene expression are rare. expression profiling revealed large-scale changes in gene expression associated with regeneration. we also sought to get insight into which transcription factors may operate at the top of the regulatory hierarchy to control gene expression in regeneration. the most significant outcome of the present work is that it lays out a roadmap for future studies of regulatory mechanisms by providing a list of key candidate genes for functional analysis. our transcriptomic analysis corroborates existing data on cellular mechanisms implicated in regeneration in sea cucumbers.
 echinoderms are emerging as important models in regenerative biology. significant amount of data are available on cellular mechanisms of post-traumatic repair in these animals, whereas studies of gene expression are rare. in this study, we employ high-throughput sequencing to analyze the transcriptome of the normal and regenerating radial nerve cord, in the sea cucumber holothuria glaberrima.  our de novo assembly yielded contigs, of which showed significant similarity to known protein-coding sequences. expression profiling revealed large-scale changes in gene expression associated with regeneration. functional analysis of sets of differentially expressed genes suggested that among the most extensively over-represented pathways were those involved in the extracellular matrix remodeling and ecm-cell interactions, indicating a key role of the ecm in regeneration. we also searched the sea cucumber transcriptome for homologs of factors known to be involved in acquisition and/or control of pluripotency. we identified eleven genes that were expressed both in the normal and regenerating tissues. of these, only myc was present at significantly higher levels in regeneration, whereas the expression of bmi- was significantly reduced. we also sought to get insight into which transcription factors may operate at the top of the regulatory hierarchy to control gene expression in regeneration. our analysis yielded eleven putative transcription factors, which constitute good candidates for further functional studies. the identified candidate transcription factors included not only known regeneration-related genes, but also factors not previously implicated as regulators of post-traumatic tissue regrowth. functional annotation also suggested that one of the possible adaptations contributing to fast and efficient neural regeneration in echinoderms may be related to suppression of excitotoxicity.  our transcriptomic analysis corroborates existing data on cellular mechanisms implicated in regeneration in sea cucumbers. more importantly, however, it also illuminates new aspects of echinoderm regeneration, which have been scarcely studied or overlooked altogether. the most significant outcome of the present work is that it lays out a roadmap for future studies of regulatory mechanisms by providing a list of key candidate genes for functional analysis. transcriptomerna-seqgene expressionregenerationechinodermnervous systemtranscription factorsinjuryextracellular matrix 