Skip to main content

A transcriptional time-course analysis of oral vs. aboral whole-body regeneration in the Sea anemone Nematostella vectensis

Abstract

Background

The ability of regeneration is essential for the homeostasis of all animals as it allows the repair and renewal of tissues and body parts upon normal turnover or injury. The extent of this ability varies greatly in different animals with the sea anemone Nematostella vectensis, a basal cnidarian model animal, displaying remarkable whole-body regeneration competence.

Results

In order to study this process in Nematostella we performed an RNA-Seq screen wherein we analyzed and compared the transcriptional response to bisection in the wound-proximal body parts undergoing oral (head) or aboral (tail) regeneration at several time points up to the initial restoration of the basic body shape. The transcriptional profiles of regeneration responsive genes were analyzed so as to define the temporal pattern of differential gene expression associated with the tissue-specific oral and aboral regeneration. The identified genes were characterized according to their GO (gene ontology) assignations revealing groups that were enriched in the regeneration process with particular attention to their affiliation to the major developmental signaling pathways. While some of the genes and gene groups thus analyzed were previously known to be active in regeneration, we have also revealed novel and surprising candidate genes such as cilia-associated genes that likely participate in this important developmental program.

Conclusions

This work highlighted the main groups of genes which showed polarization upon regeneration, notably the proteinases, multiple transcription factors and the Wnt pathway genes that were highly represented, all displaying an intricate temporal balance between the two sides. In addition, the evolutionary comparison performed between regeneration in different animal model systems may reveal the basic mechanisms playing a role in this fascinating process.

Background

The regeneration phenomenon in animals

How some animals readily manage to reform lost or damaged body parts has always intrigued man and represents a major classical biological riddle. This phenomenon was first subjected to scientific examination by early researchers/naturalists such as Abraham Trembley in the mid-18th century, who conducted elaborate experiments on a cnidarian he termed Hydra [1], and later, at the turn of the 20th century, by Thomas Morgan, who explored many plants and animals including Hydra and the flatworm Planaria [2, 3].

In its broader sense, regeneration results in the compensatory growth of any injury or detraction from the animal body, such as that which occurs during wound healing. However, researchers in the field define this term as a major restoration of body parts in a process that resembles some aspects of embryonic development (reviewed in [48]). The ability to restore organs and whole body parts is widespread in animals and also quite common in vertebrates, albeit limited in its scope to body parts like limbs and tail. The underlying mechanism responsible for this extraordinary regeneration ability is thought to be the deployment of pluripotent stem cells (in Planaria known as neoblasts), which can produce all cell types in response to injury or dissections [6, 811].

Cnidarians have remarkable regenerative abilities, and Hydra has served as a model organism for this purpose for almost three centuries since Tremblay’s time [12]. Hydra was found to be the most regenerative animal explored to date, because it can regenerate from small parts, and can even form new animals from dissociated cell masses, a phenomenon not observed in other more complex animals (reviewed in [1215]). This great regeneration potential is also thought to be mediated by stem cells, which give rise to all of the Hydra cell types [8, 16, 17]. Many of the key molecular players involved in the regeneration of Hydra have been identified, and include Wnt3 as a major head inducer, six additional Wnts, and members of the canonical Wnt pathway, including the major effector β-catenin and its DNA binding partner Tcf [1822]. The Wnt pathway is most likely to be conserved in all of the cnidarian regeneration systems and was shown to also be active in the vertebrate limb [23] and in Planarians [2429].

Polarized regeneration

The open question of how the amputated animals retain their proper polarity while regenerating still represents a mystery. In Hydra, upon mid-gastric bisection, asymmetric molecular events occur, wherein interstitial lineage cells in the head (oral)-facing part of the regenerate, undergo cell death, which leads to Wnt3 secretion and head organizer induction, resulting in subsequent cell proliferation and head formation [3032]. In the foot (aboral)-facing part of the regenerate, little cell death is observed and thus no head organizer forms, subsequently leading to foot formation. The cell death was found to be induced by the injury-driven asymmetrical activation of the ERK/MAPK pathway, which was detected only in head-facing wounds [33, 34]. However, in Hydra undergoing regeneration after decapitation, there is no asymmetrical cell death and no involvement of interstitial stem cells, which are not present in this region, but rather, the epithelial cells respond to injury signals and release Wnt3 [34].

In Planaria, Wnt1, which is necessary for executing the posterior (tail) regeneration program, does not exhibit polarized expression upon amputation and becomes enriched at the posterior-facing wounds only at later stages of the process [28]. The gene that does exhibit asymmetrical expression upon dissection and appears only at anterior wounds is Notum, which inhibits the activity of Wnt factors and also serves as a target gene of the pathway [29]. Recently, Wurtzel and colleagues reported that Notum was the first gene to show polarized expression in the time course of Planarian regeneration [35]. The execution of the head-to-tail polarity in Planaria is thus thought to be controlled by gradients of several Wnt factors that are expressed at different times after initiation of regeneration, and in different anterior-posterior (AP) domains, with feedback inhibition from Notum and other anterior inhibitors such as sFRPs [26, 28, 29]. The Wnt pathway is key to the polarized regeneration program because it establishes the A-P axis upon embryonic development in most bilaterians and in the basal metazoans at large [26]. Several other developmental pathways (e.g., Hedgehog, FGF, BMP) are also instrumental for polarized regeneration as well as many transcriptional regulators, which are asymmetrically expressed and most likely play crucial roles in this complex developmental feat [3639].

Regeneration in Nematostella

The starlet sea anemone, Nematostella vectensis (Nv), belongs to the basal cnidarian class Anthozoa and is now a well-established model organism, owing to its ease of culturing under simple laboratory conditions [4042]. It can reproduce both sexually, undergoing embryonic and planula larva stages, and asexually [40, 42, 43]. The availability of the Nematostella genome revealed its unexpected complexity in terms of its rich genetic repertoire [44] and molecular genetic studies demonstrated that the main axis of the animal (aboral to oral) bears a resemblance to the bilateral AP axis as it exhibits some differential expression of Hox genes and contains anterior markers in the aboral region and typical posterior markers (like Wnts) in the oral part [42, 4551].

Nematostella undergoes regeneration after transverse dissection; thus, it generates two complete and viable animals and also sometimes responds to injuries by forming extra heads or physa [43]. Like some of the sea anemones, the mode of its asexual reproduction is by transverse fission called physal pinching and occasionally it also performs polarity reversal; both of these events resemble regeneration. The main morphological stages that can be observed in Nematostella undergoing head regeneration from aboral (physa) fragments [52] and most recently after sub-pharynx amputation [42, 53] have been described and specific morphological stages were assigned.

The cell dynamics in Nematostella is expected to be very different from that of Hydra, since dividing cells in sea anemones, e.g. Aiptasia, are not limited to the main body column as is in Hydra, but instead are distributed in an oral-aboral gradient with the highest level of proliferating cells present in the oral/tentacle region [54]. In Nematostella undergoing oral regeneration, proliferation of cells started about a day after dissection [55]. Inhibition of cell division aborted the regeneration, demonstrating that unlike Hydra, which exhibits “morphallactic” regeneration, the mode of regeneration in Nematostella is “epimorphic-like” and relies on cell proliferation. Amiel et al. reported that cell proliferation started at 12 h after head amputation, was followed by later periods at 24–48 h, and that proliferation is necessary for the later stages of regeneration but not for the initial wound healing stage [53]).

Regarding molecular studies, like in the Hydra, the Wnt pathway has been shown to be involved in oral regeneration, as demonstrated by treating regenerates with alsterpaullone, a Wnt pathway activator. This led to ectopic oral structures growing from polyps undergoing aboral regeneration and from wounds in its torso region [56]. Recently, a genome-wide microarray-based transcriptional profiling screen was conducted to obtain the expression patterns of genes during the early wound healing process of regeneration. The healing process was observed during a few hours after puncture wounds were performed in the aboral region of juvenile polyps [57]. A group of wound-healing response genes, some of which, like metalloproteinases and several transcription and growth factors, were found to be induced also in many other species [37, 58], whereas other genes were specific for Nematostella. When pharmacological inhibitors were used, the ERK/MAPK pathway was shown to be essential for wound healing and regeneration in Nemtostella, similar to the previously reported Hydra and other wound-healing systems such as in Drosophila imaginal discs [32, 58].

To date, the oral vs. aboral regeneration response in Nematostella, which is currently the most basal metazoan model animal for the study of regeneration, has not been addressed. While Hydra and Planaria exhibit the highest regeneration capacities in terms of the extent of body rebuilding ability [5], likely due to specialized stem cells [6, 8], Nematostella may present the most primordial regeneration prototype and is therefore important for understanding the basis of this phenomenon.

Our goal in this study was to obtain insights into the transcriptional regulation programs of oral/aboral regeneration and compare the molecular time-course patterns of the two sides. We therefore used an RNA-Seq screen covering the 72 h of the regeneration process during which wound healing takes place and the animals complete the morphogenesis/patterning stage to assume their overall body shape. This screen was verified and extended by a comprehensive qPCR study of gene expression in the wound-proximal regions. We report the general trends of gene expression along this process and characterize the genes and gene pathways which responded differentially in oral vs. aboral regeneration. We elaborate on the kinetics of genes and gene pathways with polarized expression that are known from previous studies in Nematostella and other regeneration models, and reveal newly discovered genes that are likely to be important for the process. A comparison of gene expression patterns is made with the well-studied Planarian anterior-posterior regeneration paradigm, which unmasks evolutionary commonalities and differences in this major re-development event.

Results and discussion

General trends of the oral and aboral regeneration transcriptional programs in Nematostella

In order to study the transcriptional basis of oral versus aboral regeneration we performed a transcriptional profile screen in which we examined the time course of polarized regeneration. To this end, we dissected adult polyps (ca 3 months old) into two parts in their mid-body, which corresponds to mid-gastric dissection in the Hydra. The wound-proximal “blastema”–like regions from 20 oral (‘head’) or aboral (‘tail’ or physa) fragments were taken immediately upon amputation (hour 0) and at three later time points (8, 24, and 72 h post dissection) (see Fig. 1a). These time points represent early regeneration at the wound-healing stage (up to ~8 h), the active growth and cell proliferation stage (24 h), and the later stage of remodeling and morphogenesis at 72 h [42, 52, 53, 55, 57]. The morphological changes in the appearance of the bisected polyps were followed. We observed that at 8 h wound closure was complete due to wound healing processes, as reported by DuBuc et al. [57] for puncture and bisection assays performed on very young polyps [57] and by Amiel et al. [53] for adult and young polyps. At this time we did not find any visible signs of new tissue formation and the wounded mesenteries are seen to be in close proximity to the site of the cut, while at the first hours they seem to recede from the wound plane (Fig. 1b). At 24 h we observed that the mesenteries seemed to be attached to the center of the cut area in both the oral and aboral fragments as reported before for oral regeneration [53, 57]. At 72 h tentacle buds can begin to be detected in the regenerating oral side, and in the regenerating physa the mesenteries appeared disconnected from the previously wounded area and the physal end is more rounded and clear. We chose 72 h as our last time point since the general re-establishment of the body plan after dissection and the major morphogenesis processes seem to have been determined by this time. While the regenerating polyps still grow extensively for several more days until their body parts assume normal proportions, and several more weeks for returning to their original size, this is secondary to the reorganization of the body plan.

Fig. 1
figure 1

Experimental design of a time-course analysis of regeneration in Nematostella vectensis. a. A schematic depiction of the experimental procedure. Aboral “tail or physa” regeneration appears in the top row and oral “head” regeneration appears in the bottom row. The dissection planes are denoted by dashed red lines and the wound proximal parts were used for RNA extraction and the RNA-Seq screen. b. Photographs of regeneration process corresponding to the time course of the experiment. Left panels show a polyp before and after bisection at a plane marked by the yellow dashed line. The right panel shows pictures of the oral and physal-regenerating parts of the bisected fragments along the time course

We selected those regions close to the wound sites in order to focus on the parts that undergo the most conspicuous changes in their oral-aboral characteristics and in which we expected to observe a high level of response. RNA was extracted from the regenerating oral and physal sides. For each time series of regeneration two independent pools consisting of 20 animals each were created for each of the four time points (Fig. 1a); thus, we prepared 16 Illumina Hi-seq libraries. Sequencing these libraries returned an average of ~24.5 M 50 bp single end reads per sample with a standard deviation (SD) of ~9.2 M reads. An average of 73 % of the reads per sample (SD 3 %) uniquely mapped to the predicted genes in the draft Nematostella genome [44]. The number of expressed genes that were detected through read alignment was ~23000, whereas the total number of predicted Nematostella genes is ~24000.

To identify genes with time-dependent expression, a multi-factorial design analysis was carried out using generalized linear models (GLM) fitting and inference (see Methods). In this analysis we tested for the time effect on gene expression, taking into account that gene expression may also be affected by the regeneration-side. A total of 4205 genes, which corresponds to ~17 % of all Nematostella genes, displayed a significant time dependent change (BH adjusted p-value < 0.05). In Additional file 1: Table S1 these differentially expressed, time-course responsive genes (“DE genes”) are listed with their normalized number of reads together with their annotations obtained from UniProt (when available) or with the use of the Blast2GO tool. This table was the basis for the different expression patterns analyses that we present in this study.

We first examined the global nature of the gene response pattern in the oral and aboral regenerating parts. A hierarchical clustering analysis of all the DE genes indicated that more genes (about 60 %) were downregulated following dissection compared to the upregulated genes (Fig. 2a). Within the two groups many blocks of variation in one or more time point were evident, indicating temporal-, and tissue-specific differential gene expression. Similarly, sample-wise hierarchical clustering according to the temporal and spatial variables (Fig. 2b) indicates that gene expression patterns of the two tissues were largely shared at each time point, and that the different time points exhibited more distinct expression patterns from each other. As expected, the expression profiles of the two cut ends at time 0 were the most similar to each other and differed from all the ensuing regeneration time points. It was also evident that the 24 and 72 h expression profiles are more similar to each other than they are to 8 h (Fig. 2b), which is in line with previous studies defining the first 6–8 h as a distinct wound-healing stage [53, 57]. The similarity in the overall expression patterns of the two regions of regeneration, at each individual time point, suggests that a large part of the regeneration response is of a general nature which may govern general wound healing and tissue growth at large, and that side-specific differential expression patterns may be more limited. Interestingly though, the oral and physal 72 h samples were more diverse than in the other time points as reflected in their branch lengths (Fig. 2b), which may indicate that at 72 h there is an enriched expression of genes responsible for oral or aboral specific tissue traits compared to a more general earlier wound healing and initial regenerative response.

Fig. 2
figure 2

Clustering analyses of the oral and physal time course. a. Gene-wise hierarchical clustering of the differentially expressed (DE) genes. The yellow rectangle highlights genes being largely upregulated during regeneration. The blue rectangle highlights genes largely being downregulated during regeneration. The color scale indicates z-scored expression values b. A dendrogram of sample-wise hierarchical clustering of the samples from the different regeneration-sides and time points. c. The six largest clusters of gene expression in oral and physal regeneration. Clusters are ordered by the number of genes in each cluster which is indicated in the bottom left corner. The number in the top left corner indicates the STEM cluster number. The full list of clusters can be found in Additional file 2: Figure S1 and a list of the genes and which cluster they belong to can be found in Additional file 3: Table S2

In a study performed in Planaria, in which the transcriptional profiles of head and tail segments following bisection were compared by RNA-Seq analysis [38], the two parts were highly diverse at the early hours and converged to a similar pattern at the later time points. In this screen, similar to our screen there were also more downregulated than upregulated differential transcripts along most of the time course. Unlike our protocol, however, the Planaria study analyzed the entire regenerating fragments and not the wound-proximal regions as in our study, which explains the arrival to a similar state of gene expression towards the culmination of the regeneration process.

In order to identify and characterize the distinct expression profiles presented by each of the regeneration parts, we used the Short Time-series Expression Miner (STEM) software. This software is specifically tailored for high-throughput time series gene expression datasets having small number of time points with few replicates [59]. Using STEM clustering algorithm, we grouped the DE genes to clusters representing the major expression profiles in the oral and physal sides (Fig. 2c). Similar to the hierarchical clustering, the STEM clustering showed that the largest and most prominent cluster was of consistently downregulated genes. Most of the clusters in one of the sides have similar counterparts in the other side and indeed there are many shared genes in these clusters, however, the gene numbers in the clusters vary between the sides (Fig. 2c) due to side specific gene trends. The complete profile of all oral and physa clusters and their genes is presented in Additional file 2: Figure S1 and Additional file 3: Table S2.

Further analysis of the general time course trends revealed that about half of the DE genes could be defined as being always up or down-regulated throughout the sampled time points of regeneration, while the remaining genes displayed more complex expression patterns, with one or more changes in the direction of expression along the time course (Fig. 3a). We characterized a gene as being always upregulated if all time points of regeneration showed an expression level higher than that of hour 0, and always downregulated if expression in all time points were lower than in hour 0. When examining the expression modes in the overall always up and downregulated genes, the largest group consisted of genes that behaved similarly in both sides (~50 %), while the rest were either oral or physal specific (Fig. 3b). Of the genes that are regulated in one side, the oral side contained more of these and this was particularly evident among the downregulated genes. A list of the genes and their trend of expression can be found in Additional file 3: Table S2.

Fig. 3
figure 3

A general analysis and comparison of side and time course expression trends. a. Distribution of genes in either the oral or physal regenerating sides according to their time course response pattern: always upregulated, always downregulated or other, having a more complex trend. Genes were considered always up or downregulated if all time points of regeneration were at least 1.6 fold higher or lower than hour 0, respectively. b. Distribution of genes being either always upregulated or downregulated according to their oral or physal response in regeneration. c. The number of oral or physal genes being up or downregulated at each time point along regeneration. Genes were considered up or downregulated at each time point if their expression was at least 1.6 fold different compared to their expression level in the previous time point d. Examples of qPCR results of highly DE genes that show a similar expression pattern in both sides during early regeneration and a prominent divergence during later regeneration

For the comparison between time points we examined whether expression was being up or downregulated compared to the previous hour of regeneration. At the early regeneration stage (8 h), the majority of the DE genes (~75 %) displayed the same expression trend (up or down) in both the oral and physal sides (Fig. 3c). During the later stages of regeneration, a minority of genes exhibited such a parallel expression pattern and the majority showed different trends of expression in the two sides. Thus, the early response was more similar between the sides and also comprised a larger number of responsive genes, while at the later stages we observed a larger proportion of side specific genes. There was a cumulative increase in side specific upregulated genes up to 72 h whereas the number of side-specific downregulated genes decreased dramatically between 24 and 72 h, suggesting that gene upregulation characterizes tissue specific regeneration (Fig. 3c).

Similarly when we followed the time course of individual genes which displayed prominent differences in expression between the sides, we often detected similar expression at both sides in early regeneration followed by a divergence in the expression pattern. This fits the general trend we observed in our global analysis of diversification with time. An example of such an expression pattern, verified by qPCR at a higher time resolution, can be seen for the gene Wntless (Fig. 3d), which is an important regulator of Wnt factors and was implicated in polarized regeneration in Planaria [60]. This profile is typical for many of the polarized genes, both oral as shown for Wnt4 and Notum, and aboral as shown for OtxC (Fig. 3d). Thus, it seems that often the initial gene response in the two sides is similar and with time the appropriate pole-specific gene expression pattern is established. This suggests an interesting strategic mechanism of “repair missing parts first and implement polarization later”. However, the question remains whether this early response is due to a dual functionality of these genes with an early non side-specific wounding/regeneration function and a later side specific role, or whether the similar expression pattern in the immediate response is due to the animal creating a reserve of available strategic transcripts for use in the regeneration response, which will be utilized eventually only in the appropriate side with the establishment of polarity. This scenario greatly differs from the results reported for a head vs. tail RNA-Seq regeneration screen in Planaria, where the transcriptional profiles of many genes were very different at the start and later on converged to a similar pattern [38]. However, this difference is likely due to the experimental strategy used in the Planaria study of examining the entire fragments rather than the proximal regenerating regions. Yet, the initial similarity in expression between the sides in our screen does not mean that some key regulators are not differentially expressed from the start in Nematostella regeneration and indeed, we detected such early asymmetric genes, which may participate in determining polarity already at an early stage as will be reported ahead.

Characterization of side specific gene responses in regeneration

Next, we embarked on our main aim of elucidating the molecular basis of the oral vs. aboral (physal) regeneration programs, which is paramount to understanding the nature of the polarized regeneration logic. In order to characterize differentially expressed genes as tissue-specific we tested several methods of correlation analysis of the time course results. We found that the concordance method for correlation calculation (CCC) best revealed differences in both geometry and levels of gene expression. According to this analysis, genes with a low CCC are less similar and more side-specific in their expression pattern. A list of all time course responsive genes in the screen according to their CCC score, together with their oral or aboral nature and their response ratios in the time span of their largest change, is found in Additional file 3: Table S2. For a systematic identification of functional gene groups we considered genes to have a low correlation if they had a CCC of 0.6 or lower and analyzed them for enriched GO terms (as described in Methods). 1804 genes out of the 4205 time responsive genes passed this cutoff. In Table 1 the top most differentially expressed genes out of this list with the above data is shown. It can be clearly observed that the table contains a striking abundance of two types of gene groups: transcription factors and Wnt pathway members. The major GO groups that were detected in the list and their enrichment fold are depicted in Fig. 4 and a full list of the major GO groups with low CCC can be found in Additional file 4: Table S3. Some of the GO categories are redundant while some are a combination of other groups, for instance the large “extracellular region or matrix” groups are composed of the protease/peptidase groups, the Wnt pathway group and the chitin related group. The main functional groups that we found to be enriched among the low concordance genes are described:

Table 1 Differentially expressed oral vs. aboral genes
Fig. 4
figure 4

Major Gene Ontology (GO) terms with highly differential oral to aboral expression. Fold change of the most prominent GO terms that are enriched in genes with low concordance between the sides is shown. The numbers next to bars represent Benjamini-Hochberg (BH) adjusted p-values for the enrichment of the GO term

Peptidases and proteases

These gene families are known to take part in tissue remodeling in many regeneration models ([57] and references therein). These are necessary for the wound healing and injury repair program [61] taking place in both sides. This entails breakdown of the mesoglea, the basal lamina, at the amputation site, and its reconstruction after wound closure, as well as release of growth factors for the following reconstruction phase. Massive protease involvement in the early stages was observed in the previous study on wound healing in Nematostella [57], and indeed our screen shows a wealth of many types of such gene groups. Interestingly, but as may be expected, the “proteases”, “peptidases” and “extracellular regions” categories were enriched in the differentially regulated genes groups. This can indicate the widespread differential tissue and extracellular matrix (ECM) sculpting programs and growth factor shedding in the two sides (Fig. 4). Predominantly, the group of zinc metalloproteinases (MPs) of the ADAMTS (A Disintegrin And Metalloproteinase with Thrombospondin Motifs) type [62, 63] were differentially expressed in the two sides, with the majority highly expressed at the oral side at some or all time points, while several were expressed higher in the aboral side, as shown in Additional file 5: Table S4 and Additional file 6: Figure S3. Some other MPs showed polarized expression as well, including several members of the Astacin group, that were reported to be essential for both head and foot regeneration in the Hydra (reviewed in: [64]), and also the membrane type ADAM10, which is a sheddase of cell surface proteins such as growth factors. Nematostella matrix metalloproteinase (MMP) 17 was detected to be oral responsive. Its Hydra ortholog HMMP was also reported to be expressed in development and regeneration [64], and the Planarian ortholog was found to be enriched in the anterior regeneration program of the worm [38]. Regarding MP inhibitors, a tissue inhibitor of MP (TIMP) gene was reported before in Nematostella wound repair [57], to be expressed early in both the oral and aboral sides after amputation, and in our screen displayed an early higher transient aboral expression followed by an oral only expression at later times, suggesting a complex regulation of tissue consistency along the regeneration process. An additional inhibitor with an exclusive early oral expression pattern was also identified (Additional file 6: Figure S3). This complex and dynamic expression pattern of several proteinase types that are at interplay with their inhibitors is likely to be instrumental for the formation of oral and physal specific features.

The noticeable abundance of multiple ADAMTS species is an interesting avenue for further research in order to identify the substrates for these enzymes, which will allow a better understanding of their functions. Some of the orthologs we have detected which include multiple ADAMTS6-like but also 16 and 17 are orphan, while ADAMTS9 is a proteoglycanase, ADAMTS2 and 3 are procollagen peptidases and 7 and 12, are cartilage oligomeric matrix protein (COMP) proteases [57], with COMP also appearing as a down regulated gene in the screen. It will be necessary to verify the identity of the substrates and importantly – explore the patterns of ADAMTS expression and activities during the regeneration course. Many genes with a high correlation score were also found to be associated with the “proteases” and “peptidases” GO terms although these terms were not significantly enriched as for the side specific genes of these categories (data not shown). This is to be expected since there is a large general remodeling response of the tissues in both sides as was described above.

Transcription factors

The “DNA binding” and “transcription factor” groups were found to be prominent amongst the side specific genes in our screen both in terms of gene numbers and fold enrichment (Fig. 4). In Table 2, the top differentially expressed transcription factors (TFs) were marked according to their type and differential expression pattern in the oral or aboral sides and the full list is shown in Additional file 7: Table S5. Of these factors the homeobox genes were the most dominant group with multiple representatives (about 1 out of every 3 DE TFs), many of which displayed higher induction of expression. This may be expected of developmental “toolkit genes” that are known to be involved in developmental patterning and in evolutionary changes of the developmental process [65]. There are three Hox genes within the homeobox genes. However, the subgroups with the most members (4 each) are NK-homeobox type genes (4), which are considered to be the most ancient type of homeobox genes [66] and the Mox genes, which are the most highly induced TFs in our screen. Interestingly, the Mox family is represented by twice the number of genes in Nematostella as compared to human [50]. We can speculate that there may be a connection between this high copy number and strong induction of the Mox genes to a specific function in regeneration, which will be explored in future studies. The screen also included three aboral expressed Six genes including Six3/6, which is a master aboral pole regulator in Nematostella development [51], a key anterior development gene in bilaterians in general, and is necessary for forebrain formation in vertebrates [67].

Table 2 Differentially expressed homeobox genes and other transcription factors

Other types of TF genes we detected include several factors of the following types: Forkhead (7 members), Sox HMG box domain (7), SP zinc-fingers containing (4), T-box (4), Ets domain (3), Krueppel-like (2) and multiple TFs with one copy such as the bHLH factor Twist-related2, the Runt-domain Runx, p63 and HIF-1.

Most of the differentially expressed homeobox genes and other TFs are more highly induced in the oral side (67 and 61 % respectively), which can reflect the higher complexity of the structures formed at the oral side. When comparing oral side expression patterns to those of Hydra head regeneration, we can observe similar expression pattern along the process, such as the orthopedia-like (otp) and Sox14 genes which are late expressed oral genes both in Hydra as well as in this screen. On the other hand, some factors show a strikingly different pattern between these two cnidarian groups, such as the paired-type homeobox repressor Dmbx gene, which is downregulated in regenerating Hydra heads [68] but upregulated in Nematostella oral regeneration, and to a lesser degree also in the aboral part. It is more difficult to compare Hydra and Nematostella aboral regeneration since there are few reports regarding aboral gene responses in Hydra regeneration. The comparison of gene expression accompanying the regeneration process of these two systems would be interesting, as their aboral parts differ in shape and function; Hydra have an attachment disc whereas the physa of Nematostella is simplified and adapted to burrowing.

The Wnt pathway

This major developmental signaling cascade was found to represent a most prominent group of genes that was notably enriched in our screen and for which the expression pattern of many of its members was highly polarized. This pathway is known to be essential for establishing A-P polarity in the developmental stages of many diverse animals [26]. It is also known to play a key role in animal regeneration and specifically in cnidarians [13, 19, 20, 22, 57, 69]. An examination of the expression patterns of these genes along the time course of regeneration verified the expected oral predominance in a large majority of them.

Three distinct clusters of Wnt pathway genes could be identified based on clustering of their oral regeneration expression patterns (Fig. 5a). The largest cluster had an expression pattern in which its genes were upregulated at 8 h and therein either remained upregulated thereafter or showed a decrease in expression but always above their initial time 0 levels. We termed this cluster “early upregulated” (Fig. 5a1) and it contained Wnts4, 11 and 16, Frizzled1 and 10, intracellular Wnt pathway members like Tcf/Lef-1, Axin and naked-cuticle homolog1(Nkd1) as well as the Wnt secretion regulator Wntless and the Wnt inhibitors Notum and Dkk (Additional file 8: Table S6, Additional file 9: Figure S4 and Fig. 3d). Another cluster included genes that were initially downregulated and then upregulated in later time points with some returning to the baseline by hour 72 (Fig. 5a2). This cluster contained Wnt6 and 10, Frizzled2, two Wnt co-receptor LRP adapter protein genes and two protocadherin fat4 like genes, which are associated with planar cell polarity (PCP) activity and growth control via the Hippo pathway [70, 71]. The third cluster contained genes that were upregulated at a later stage (Fig. 5a3), beginning at hour 24, and was thus termed “late upregulated”, with only three Wnt factors: Wnts1, 2 and A, which are closely related in a Wnt family cladistics analysis as reported by Hensel et al. [72] but are more distantly related according to Kusserow et al. [73]. This clustering of different pathway components indicates a pattern of co-regulation, elucidation of which may be important for deciphering the different stages of the oral vs. aboral patterning process in regeneration in the future studies.

Fig. 5
figure 5

Analysis of the Wnt pathway and other developmental signaling groups. a. Clusters of expression of Wnt pathway genes of the GO group in the oral regeneration response. Three clusters were discerned including genes of the following expression patterns (with representative genes): 1. early upregulated (Nv-Wnt11, Nv-Wnt16, fzd-10 precursor). 2. downregulated and then upregulated (Nv-Wnt6 ,10). 3.late upregulated (Nv-Wnt1, Nv-WntA). Each gene is plotted (gray) in addition to the mean expression profile for that cluster (blue). A full list of the genes and their cluster can be found in Additional file 8: Table S6. b. A proposed hypothetical schematic of the canonical and PCP Wnt pathways interaction in aboral regeneration. The question marks depict signaling interactions which are less established and are thus more assumptive. The bent arrows represent a possible negative feedback loop between cilia and the early canonical Wnt complex that is highly speculative at this stage

Regarding the physal side, the majority of Wnt pathway genes showed lower and often only very little expression as compared to the oral side and thus we did not manage to cluster these genes’ expression patterns. A notable exception to the oral preponderance of the Wnt pathway members was the Frizzled5 gene which was strongly upregulated at 8 h in the aboral part, while its oral expression was downregulated. Nematostella Wnt5, reported to function as a non-canonical PCP/convergent extension regulator in experiments conducted on Xenopus embryos [74], displayed a unique higher aboral as compared to oral expression pattern along the regeneration time course, except at 72 h when its oral levels were higher. Thus Wnt5 may signal via its cognate receptor Frizzled5 in the aboral side to allow elongation of the physa after the wound repair stage (see hypothetical model at Fig. 5b).

In the oral side the late Wnt5 expression may likely induce evagination of tentacles as was reported in an elegant work on Hydra [75]. Wnt5 is expressed in Nematostella at the base of tentacles [73], however the PCP function here may also be carried out by the non-canonical Wnt11 [74], which is also expressed specifically in the tentacles at the juvenile polyp stage [76], and is an upregulated oral specific gene in our screen (Additional file 8: Table S6 and Additional file 9: Figure S4). It is not clear which Frizzled may mediate the presumptive Wnt/PCP program at the oral side but a possible candidate is Frizzled2.

Of interest, a Wnt/PCP pathway transmembrane co-receptor Ror2, exhibited a high aboral expression pattern similar to that of Frizzled5, which can indicate that this gene is also an evolutionarily conserved morphogenetic gene active in our speculative model of aboral PCP patterning at early regeneration (Fig. 5b). Regarding inhibitors of the pathway, besides the aforementioned Notum, we detected a Dickkopf (DKK) factor [76, 77] with an early oral expression and a steadily rising aboral expression reaching high levels at 72 h. We also found APCDD1 [78, 79] (Table 1), a little known membrane-residing inhibitor reportedly studied only in humans until now, which is downregulated upon regeneration and whose expression increases only at the aboral side at the late stage. This expression mode may imply regulation at both sides of the canonical Wnt pathway at the earlier stages and a late establishment of the steady state Wnt oral to aboral gradient.

The dynamic and intricate expression mode of the Wnt factors that is usually higher in the oral regenerating regions but is often accompanied with a lower and shorter aboral activity may also imply that these genes play non-redundant roles in specifying different regions and structures along the oral-aboral axis in regeneration, partly recapitulating the embryonic/larval developmental patterning of this major axis in cnidarians [19, 22]. In Hydra Wnt3 is the leading Wnt factor in head formation, budding, and head regeneration, and is expressed in the hypostome at the apex of the head [20, 22, 80]. In Nematostella Wnt3 is expressed in the pharynx of the polyp [76], and we did not detect a response of the Nematostella Wnt3 ortholog in our regeneration screen. Thus, the identity of the lead early Wnt factor in Nematostella regeneration remains to be determined. A possible candidate is Wnt4, which was upregulated very early on and had both the highest absolute reads and the largest fold induction level along the regeneration time course. This said, other early Wnt factors may induce the Wnt response and functional studies could resolve this issue in the future. Early onset of expression was detected for Wntless, Frizzled1, Lef-1 (Additional file 9: Figure S4 and Fig. 3d) and many other Wnt pathway components, demonstrating that its signaling activity is turned on early in the process at both sides. We did not detect induction of β-catenin response but this does not preclude its activation at the protein level along with other members of the pathway, which may or may not be also induced transcriptionally.

In summary, the above reported observations allow us to propose a speculative model, in which the less explored aboral regeneration plan involves an early transient canonical Wnt pathway activation which may specify the aboral domain as such and eventually elicit cell proliferation (Fig. 5b). This stage may then be followed by the activation of the non-canonical (PCP) pathway either directly by the canonical pathway or indirectly by aboral specific TFs. Activation of the Wnt-PCP pathway may induce elongation of the aboral part and may participate in the inhibition of the canonical pathway at later times together with the canonical inhibitors such as DKK (Fig. 5b). Thus an intricate interaction between Wnt pathways may lie at the basis of the regeneration process and other cues may also be involved such as ciliogenesis as described ahead. Future experiments in which the canonical and PCP Wnt pathways will be manipulated using pharmaceutical agents or genetic tools will allow establishing our proposed model.

Additional developmental signaling pathways

Although the GO analysis did not point to other enriched pathways showing differential expression, we manually searched for highly differentially expressed genes belonging to other families of signaling factors known to be involved in development and regeneration, such as the BMP, Hedgehog, FGF and Notch pathways (see Methods). Indeed we have located much fewer genes with a low CCC score (<0.6) for each of these pathways as compared to the Wnt pathway (Additional file 10: Table S7), explaining why they did not show up in the pathway GO analysis. Nevertheless, there were several individual members of these families that were considerably differentially expressed (Additional file 10: Table S7). These included the following: 1. two genes coding for ligands of the FGF pathway. 2. Hedgehog and its receptors patched domain-containing protein 3-like and Smoothened and the zinc finger Zic1-like. 3. the Notch pathway Jagged1 and Delta-like ligands, a HLH domain Hes4-like TF gene and a Hey, a Notch transcriptional repressor gene. Thus, these pathways may also contribute to the polarized regeneration process even though their participation seems to be less extensive than that of the Wnt pathway, which displays many DE components at all levels of the signaling cascade. Whether there are specific roles for the different pathways in regeneration and their possible interaction with the Wnt pathway are interesting avenues for future research.

With regard to the BMP pathway, most of its ligands and receptors were non-responsive or only weakly expressed, except for BMP2/4, which was downregulated during the course of regeneration (Additional file 10: Table S7). Of note was the Smad4 effector gene which did show induction of expression, an observation that is in line with recent work on Hydra head regeneration [68]. Interestingly, the Tolloid gene (Nve-Tld [81]), which encodes for an Astacin type metalloproteinase that has multiple substrates including the BMP2/4 antagonist Chordin [82], showed a shift in expression from high oral to aboral with time. However, when we examined BMP antagonist like genes, we found that while Chordin and Noggin orthologs had very low time course responses, kielin/chordin-like (KCP) [83], a gene whose product enhances BMP signaling, displayed high aboral expression. These findings may suggest that the BMP pathway functions differently in regeneration as compared to its activity during development in Nematostella.

Chitin-related genes

Surprisingly the group of genes that was most highly enriched in our analysis of differential expression was the chitin associated genes (Fig. 4 and Additional file 11: Figure S5). This group contained mostly chitinase genes, as well as von willebrand factor d and egf domain-containing genes that are also predicted to be secreted and to harbor a chitin binding domain found in animal chitinases. These enzymes have been reported to be expressed in the endoderm of the body column of three Hydra species [84] and were suggested to function in prey digestion. The first cnidarian chitinase that was cloned and studied in the hydroid Hydractinia, was rather found to be highly expressed in the ectodermal layers of the stolons of the colonies and of the basal (aboral) part of the polyps, and is proposed to act in immunological protection of the latter from chitinous pathogens [85]. However, some of the chitin associated genes that we located in the screen were chitin synthases and chitin deacetylases, which may strongly hint that Nematostella contains chitin, as do other anthozoans such as corals and a shell producing sea anemone [86].

Interestingly, all of the above chitin associated genes show a similar time course expression pattern upon regeneration; at time 0 the oral facing proximal region has a higher level of expression as compared to the aboral facing region, which is typical of genes exhibiting a sharp gradient in the mid-gastric zone. At later time points both of the oral and aboral expression levels decline, with some of the genes displaying late upregulation of aboral expression, which likely restores the initial aboral high/oral low gradient (Additional file 11: Figure S5). Thus, while the function of these chitin-related genes is yet unexplored, our study indicates that the activity of this gene assembly is quickly downregulated upon regeneration and is restored towards its later stages in order to re-establish the aboral specific steady-state function.

Tubulin and cilia associated genes

An additional GO group of interest was that of microtubule cytoskeleton, which stood out in the side-dissimilarity of the expression patterns of its members, although it was not one of the GO groups that were found to be enriched in the screen (Fig. 4). Further inspection of this group pointed to a number of genes, among them two tubulin-alpha3 genes, which are highly expressed in the aboral regenerating region at the early stages (Fig. 6 and Additional file 12: Table S8). This was surprising as these genes as well as many other microtubule-associated genes, are more highly expressed in the oral part in intact animals, compared to the aboral. Most interestingly, there was also high aboral expression of some cilia-related genes during regeneration (Fig. 6). These include IFTs (intraflagellar transport proteins), that are important for cilia formation and function, and several tubulin polyglutamylases, which modify tubulins into ciliary-enriched forms [87]. Another early aboral expressing gene we identified was the Bardet-Biedl syndrome protein 5 homolog (BBS5) gene, which is part of the basal body linked complex of proteins (BBSome) [88] that are mutated in genetic diseases linked to ciliary function, termed ciliopathies. All of this evidence suggests that an early ciliogenesis reaction occurs in the aboral regenerating region of Nematostella, which may imply that cilia are involved in aboral regeneration. This can be attributed to several possible mechanisms such as the directed flow of secreted factors in the gastrovascular lumen or cilia-based signaling.

Fig. 6
figure 6

Tubulin and cilia-associated genes in regeneration. A schematic depicting the possible involvement of tubulin and cilia enriched gene products in aboral regeneration. The axoneme of the cilia together with the basal body and its associated BBSome complex (see text) are illustrated. The graphs of representative early aboral enriched genes that are known to reside in cilia are shown: one of the two Tubulin3α whose products are part of the microtubules; IFT80, which is part of the dynein-driven IFT protein complex; and BBS5, which is a member of the BBSome coat complex that is linked to the basal body of the cilia. The inset in the Tubulin3α panel depicts the expression level in different body parts of intact polyps: T, tentacles; O, oral ring and pharynx; M, mesenteries and P, physa. IFT, intra-flagellar transport; BBSome, Bardet-Biedl Syndrome protein complex

Regarding the latter, it was established that components of the hedgehog pathway reside in cilia, and that cilia are essential for the function of this pathway in mammalian development [89]. We did not find, however, any conclusive concerted aboral specific induction of the major hedgehog pathway genes upon regeneration (Additional file 10: Table S7) but this cannot negate the possibility that there is an aboral enrichment at the protein level. Other pathways have been proposed to affect, and be affected by, cilia formation. Of these it is worthy to note that disruption in several ciliary genes elicited the canonical Wnt pathway in some developmental scenarios [90, 91]. Thus, ciliogenesis may be expected to have an effect on the inhibition of Wnt signaling. Accordingly, we can speculate that if indeed our tubulin and cilia-associated expression data indicates that cilia formation takes place in the aboral regenerating animals and this may explain the early down-regulation of the canonical Wnt response which we observed as part of the aboral reaction (Fig. 5b).

An evolutionary comparison between Nematostella and Planarians

Finally, we compared the expression patterns of the key genes we identified in our screening to those reported in the well-studied regeneration plan of the Planaria, where the head vs. tail regeneration program has been studied in depth [24, 27, 28, 38, 92]. To this end, we identified the Nematostella orthologs of the known Planaria anterior and posterior marker genes which are also induced during regeneration in a polarized fashion. When no bona fide orthologs could be found we added the closest homologous genes that also show polarized expression to complete the picture (Additional file 13: Table S9 and Fig. 7). Many of our oral-induced genes corresponded to the tail/posterior genes/markers of the Planaria and inversely, the aboral genes matched those of the head/anterior Planarian genes (Fig. 7). This relationship between the cnidarian oral/aboral axis and the anterior/posterior axis of bilaterians has been reported recently [51], where the aboral part of Nematostella shows a molecular resemblance to the anterior part of bilaterians, and inversely, the oral region has similarity to the posterior bilaterian side. Thus, the aforementioned master aboral territory gene Six3/6 [51] and the OtxC gene that are highly expressed in aboral regeneration (Fig. 3d) have orthologs in Planaria that are anteriorly expressed (Fig. 7). We also detected a higher aboral expression of FGF1a [93] (Fgf20 in [57]) and FGF-R1- like with a corresponding head expression of the FGF-receptor like ndk, ndl-3 and ndl-4 Planaria genes. Interestingly, we also found several oral enriched Nematostella FGF factors induced in regeneration (Fig. 7), however, the identity of their receptors as well as of the aboral FGF factors and their effect in each side awaits further studies. Interestingly, we found that the forkhead factors FoxA2 and FoxB in Nematostella, which were reported to be expressed in the pharynx during larval development and in young adults [94, 95], have a late high oral expression pattern. Accordingly, they most likely correspond to their Planarian counterpart FoxA, which was shown to be a pharynx marker [96] and a factor that is essential for the formation of this organ in Planaria and perhaps in all eumetazoans that contain pharynx/foregut [97].

Fig. 7
figure 7

Comparison of “head” vs. “tail” gene expression in the regeneration of Nematostella and Planaria. A cartoon of the two animals is presented with the oral (head) part of Nematostella aligned with the posterior (tail) part of Planaria according to the comparison in which these ends express an assembly of Wnts and other universal bilaterian posterior markers, while the aboral part resembles the Planarian anterior (head) part in the expression of universal anterior markers such as Six3/6 and Wnt inhibitors like sFRPs [51]. Some of the markers specific for the main axis ends are depicted with Wnts and Wnt pathway genes in red and orange, FGF pathway members in brown, homeobox factors in light blue, forkhead (Fox) factors in purple, and other TFs in green. An anterior Planaria regeneration marker and protease, proprotein convertase 2 (PC2) and its ortholog in Nematostella, the late aboral Pep8 gene are marked in grey. Highlighted in light green are bona fide orthologs that are expressed as expected according to this comparison, while yellow marks polarized genes that are expressed in an ‘inverted’ position, which can indicate evolutionary change in the regeneration program. aThe FoxD1 gene is induced at the oral regenerating region in Nematostella however it is aboral in embryonic/larval development and becomes orally expressed upon metamorphosis and tentacle budding [95]. b Dickkopf (DKK) is expressed higher orally at 8 h after the induction of regeneration and higher in the aboral part later than 24 h. c Runt-1 is induced in anterior sites upon eye and neural brain regeneration in Planaria but is also expressed in posterior wound sites [35, 37] and https://radiant.wi.mit.edu/app/

Notably, several Wnt genes are expressed in the posterior side as well as in the tail forming side of wounds in the Planaria and in the corresponding oral regenerating side of Nematostella, while the Wnt inhibitor sFRP is expressed in head facing wounds of Planaria and in the aboral regenerating side of Nematostella (Fig. 7). Although we cannot determine at this stage that the orthologs have a similar role in the regeneration process in these two phylogenetically remote animals, some similarities can be observed. For example, the involvement of several Wnts exhibiting early and late induction times such as Wnt4, which is initially expressed at both the regenerating sides (Fig. 3), can be likened to the key Planarian Wnt1 (WntP-1) that is initially expressed at both wound faces before becoming enriched at the posterior side [28], and induces later expressed Wnt genes. Additional examples are Frizzled1 and 10 which are orally expressed in Nematostella and may correspond to frizzled4 which is posterior in Planaria. We did not observe any polarized expression of the major effector β-catenin in Nematostella, while in Planaria it is a key posterior gene in regeneration [24]. Whether these homologous genes function similarly should be established in future perturbation experiments where the phenotypic and molecular effects can be compared.

However, some orthologs of the two species are inversely expressed with respect to the above rule. This could be seen for a group of side specific transcription factors including orthopedia (Otp), Zic, SoxB and FoxD which we found to be oral in Nematostella regeneration but are anterior in Planaria. Besides these factors, Wnt2 in Planaria is unusual as it is more anteriorly expressed in relation to all its other paralogs that are posterior oriented [28]. Most interestingly, the Wnt inhibitor Notum is orally expressed in Nematostella, while it is anterior in the Planaria [29]. This inversion of Notum expression in Planaria is of special interest as it was very recently found to represent the most early asymmetrically expressed gene in the Planaria [35], and thus orchestrates the subsequent polarization of Wnt factors [29], which is crucial for establishing the polarity of the A-P axis. Notum was found to be posteriorly induced in regeneration of the acoel worm Hofstenia [98], which resembles the Nematostella expression pattern with regard to the Wnt factors expression, and thus the Planarian regeneration scenario may be derived as compared to the basal state. It is interesting to note that we did not detect genes that are expressed in the aboral part of Nematostella which have their orthologs expressed in the posterior domain of the flat worm. A simple explanation is that besides the Wnt factors not many genes were defined in this part in Planaria, while the anterior region is more complex in term of structures and has also been more explored. However this finding could also suggest that the oral part of Nematostella has indeed similarities to both the anterior and the posterior parts of the Planaria and may attest to grand maneuvers which occurred in axis polarity along evolution.

Such changes in the polarity of gene expression along the major axis may reflect macro-evolutionary changes owing to differences in the body plan of bilaterians compared to that of radial-like symmetry in cnidarians, or in the regeneration mechanisms employed by these two evolutionary distant animals that also exhibit very different life styles. Understanding the basis for these evolutionary changes in the regeneration plan will also benefit from elucidating the gene networks involved and insight may be gained regarding the nature of changes that sculpt different body plans and the evolutionary process that shaped them. Eventually, the study of the mechanisms underlying regeneration in multiple animal models may be a most insightful approach to elucidate the molecular basis of this fascinating developmental process and may prove beneficial for future innovative injury treatments in man.

Conclusions

The re-formation of polarity upon the course of regeneration is a long time biological mystery but in general terms it is a specialized case of the formation of the polarized body of animals that occurs upon development. While the overall response in the two sides seems to be very similar with many genes being downregulated, there were nevertheless many genes which were expressed in a side-specific manner. Three main groups of genes were differentially expressed between the sides, namely proteases, transcription factors and Wnt pathway members, with the last especially demonstrating a complex dynamic expression pattern, suggesting an interplay between the canonical and non-canonical pathways at both sides. Among the other polarized groups, the microtubule and cilia associated genes stood out as novel potential players in the control of polarity induction. From an evolution wide comparison of the Nematostella regeneration plan to Planaria, we can discern a similarity of a conserved developmental pathway and other key genes along the presumed corresponding sides of the main axes of these animals, while multiple genes also show an inverse pattern of expression. Thus, like development in general, the whole body regeneration plan seems to harbor a common general trend with variations evolved to accommodate the specific body shape and structure of different animal groups.

Methods

Nematostella vectensis culture and regeneration

Nematostella polyps were raised and kept in 30 % artificial sea water (ASW) at 18° and fed freshly hatched Artemia three times a week. Adult 4 months old polyps of 1.5–2 cm length were used in the regeneration experiments, prior to which their water was changed and they were not fed for 3 days. Animals were relaxed in 7 % MgCl2 in 30 % ASW for 10 min prior to bisection halfway between the oral and aboral ends. For the hour 0 time point, a 1 mm wide section was collected from the area of the cut immediately following bisection. Regenerating Nematostella were kept in 30 % ASW at 18° and the regenerating ends were sampled at 8, 24 and 72 h following the bisection. At those specific time points a 1 mm wide section was taken from the areas undergoing regeneration and collected in Tri-reagent (T9424 Sigma), homogenized and frozen at −20° until further use. The animals assayed were cut twice: once upon the initial bisection and once again at the time of sample collection. 20 animals were used per sample with two biological replicates per time point, for a total of 16 samples (4 time points x 2 regeneration areas x 2 replicates). For quantitative PCR experiments the same procedure was performed but Nematostella were sampled at 0, 1, 2, 4, 6, 8, 12, 24, 72 and 96 h post amputation. Pictures of the area of regeneration (Fig. 1) were taken at times corresponding to those sampled for sequencing and used a Nikon DS-Fi1 with DS-U2 controller (resolution 2560X1920).

According to the ethics instructions of the Hebrew University of Jerusalem most invertebrates including cnidarians such as Nematostella are exempt from the ethics committee inspection.

RNA isolation and sequencing library preparation

RNA from the tissue samples was processed in Tri-reagent as described [99]. Libraries were prepared from RNA samples using Illumina TruSeq RNA Sample Preparation Kit v2 (Illumina # RS-122-2001) according to the manufacturers’ instructions. This protocol utilizes poly-A selection of mRNA from total RNA. Library concentrations were quantified with Agilent 2200 Tapestation, diluted to 10 nM and pooled with 8 samples in each pool. The samples were sequenced using 50-bp single-end reads on two Illumina HiSeq2000 lane and TruSeq v3 flow chambers at the The Technion Genome Center, Haifa, Israel.

Processing and RNA-Seq data analysis of Illumina sequencing reads

Nematostella vectensis reference genome was downloaded from Ensembl. Gene product annotations were retrieved from UniProt. However, proteins designated as ‘predicted protein’ were annotated using Blast2GO through blastp vs. refseq (considering up to top 20 hits, e-value cut-off 10−3).

The 50-bp reads were quality assessed using FASTQC and subsequently aligned to the N. vectensis genome using the STAR aligner [100]. Only those reads uniquely aligned to the reference genome were subsequently used for analyses. Raw read counts per gene were counted using HTSeq-count. Subsequent normalization and differential expression analyses were carried out using the DESeq package [59, 101]. As we used duplicates in this experiment, we decided not to perform pairwise comparisons among pairs of time points, but instead perform a test which considers the entire time course in both regeneration-sides. Using DESeq, we applied a multi-factorial design, with the time and regeneration-side as the biological factors influencing gene expression. For each gene, we checked whether the time in general had a significant effect on its expression, while accounting for the effect of the side of regeneration. In short, for each gene, two generalized linear models (GLMs) were fitted using the fitNbinomGLMs function: a full model regressed the genes’ expression on both the regeneration-side and the time effect; a reduced model regressed the gene expression only on the regeneration-side effect. The nbinomGLMTest was then used to compare the two models in order to infer whether the additional specification of the time (on top of the regeneration-side) improved the fit and hence, whether the time had a significant effect on gene expression. This analysis resulted in a p-value per gene, which was further adjusted for multiple testing using the BH method. Out of >23,000 expressed genes, 4205 genes had time effect adjusted p-value < 0.05. These genes were termed “differentially expressed” (DE) genes. These genes were then further analyzed throughout the rest of this work.

Clustering analyses

Hierarchical clustering of DE genes was performed in Partek® Genomics Suite® using Pearson’s dissimilarity and complete linkage. Prior to hierarchical clustering, the normalized counts were log2 transformed (after zeros were converted to 1) and standardized (z-scored) such that gene-wise means were shifted to zero and the standard deviation was scaled to 1. Partition clustering of DE genes was performed with STEM [102] using the STEM clustering method and a maximum of 20 model profiles.

To determine up- or down-regulation among certain time points in a given side, a difference of at least 1.6 fold between the time points was required, considering average of normalized counts per time point. We chose to use a rather permissive fold change cutoff since we were looking at the changes across a few time points together so that even a slightly lower but consistently changed expression was significant. A gene was considered to be always down or upregulated during regeneration if the expression level of all time points were at least 1.6 fold lower or higher than hour 0 respectively. DE genes that didn’t fall into these categories were counted as "other". For comparison between consecutive time points a similar process of 1.6 fold change was used, albeit between each time point and its previous one, to determine the overall change in transcriptome expression levels. For Wnt family genes, hierarchal clustering was performed using the heatmap.2 function in R, and the tree was separated to discrete clusters using the cuttree function.

Concordance calculation among regeneration sides followed by Gene Ontology (GO)

To assess the concordance or discordance of the time-dependent expression pattern between the oral and aboral sides of regeneration we calculated Lin’s concordance correlation coefficient (CCC) [103] for the 4 data points of oral and physal regeneration. This calculation was performed on the 4205 DE genes defined above. CCC ranges from 1 to −1 with different consensuses as to what constitutes a poor concordance. Based on observations of the data we found that genes with a CCC of 0.6 or less showed little to no similarity between sides and 0.6 was thus chosen as a cutoff for capturing genes with different expression patterns between regenerating sides. Genes with CCC < 0.6 were defined as having low concordance, and were subsequently analyzed for enriched GO terms using the DAVID server [104106]. A full list of enriched GO terms can be found in Additional file 4: Table S3. Many enriched GO terms contained similar genes and only one from each grouping was shown in Fig. 4.

To define Nematostella genes from several pathways of interest, gene lists were taken from KEGG (Kyoto Encyclopedia of Genes and Genomes). Since the KEGG Nematostella pathway annotation is not complete, pathway genes not available in the Nematostella annotation of KEGG were supplemented by using BLAST to find genes corresponding to prominent pathway orthologs.

Quantitative real-time PCR (qRT-PCR) analysis

To confirm gene expression levels detected by RNA-Seq, quantitative real-time PCR was performed using the Applied Biosystems StepOnePlus™ Real-Time PCR System using gene-specific primers (Additional file 14: Table S10). RNA was extracted using Tri-reagent [99] and RNA integrity was verified using agarose gel electrophoresis. RNA concentration and purity was measured using a NanoDrop spectrophotometer. cDNA was synthesized using High-Capacity cDNA Reverse Transcription Kit (Applied Biosystems). cDNA samples were run using Fast SYBR® Green Master Mix (Catalog Number 4385612/4). Relative quantity calculations were performed with the ΔΔCt method, normalized against the gene Nv-Hsp70 as described [56]. The error bars represent SD of three technical repeats and the graphs are representative of at least 2 biological repeats performed on different batches of animals. In general about 75 % of the genes analyzed by qRT-PCR showed overall similarity in their expression patterns to the results of the RNA-Seq.

Abbreviations

ADAMT:

A disintegrin and metalloproteinase with thrombospondin motifs

AP:

Anterior-posterior

ASW:

Artificial salt water

BH:

Benjamini & Hochberg

CCC:

Concordance correlation coefficieant

COMP:

Cartilage oligomeric matrix protein

DE:

Differentially expressed

DEG:

Diferentially expressed gene

DKK:

Dickkopf

ECM:

Extracellular matrix

FPKM:

Fragments Per Kilobase of transcript per Million fragments mapped

GLM:

Generalized linear models

GO:

Gene ontology

GRN:

Gene regulatory network

IFT:

Intraflagellar transport protein

MMP:

Matrix metalloproteinase

MP:

Metalloproteinase

PCP:

Planar cell polarity

qPCR:

quantitative PCR

SD:

Standard deviation

STEM:

Short time-series expression miner

TF:

Transcription factor

TIMP:

Tissue inhibitors of metalloproteinase

References

  1. Tremblay A. Memoires pour servir a l’histoire d’un genre de polypes d’eau douce. Leide:Jean et Herman verbeek; 1744.

  2. Morgan TH. Regeneration. New York: The Macmillan Company; 1901.

    Book  Google Scholar 

  3. Sunderland ME. Regeneration: Thomas hunt Morgan’s window into development. J Hist Biol. 2010;43:325–61.

    Article  PubMed  Google Scholar 

  4. Sanchez Alvarado A. Regeneration in the metazoans: why does it happen? Bioessays. 2000;22:578–90.

    Article  CAS  PubMed  Google Scholar 

  5. Sanchez Alvarado A, Tsonis PA. Bridging the regeneration gap: genetic insights from diverse animal models. Nat Rev Genet. 2006;7:873–84.

    Article  PubMed  CAS  Google Scholar 

  6. Bergmann A, Steller H. Apoptosis, stem cells, and tissue regeneration. Sci Signal. 2010;3:re8.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  7. Bonfanti L. From hydra regeneration to human brain structural plasticity: a long trip through narrowing roads. ScientificWorldJournal. 2011;11:1270–99.

    Article  PubMed  Google Scholar 

  8. Tanaka EM, Reddien PW. The cellular basis for animal regeneration. Dev Cell. 2011;21:172–85.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  9. Handberg-Thorsager M, Fernandez E, Salo E. Stem cells and regeneration in planarians. Front Biosci. 2008;13:6374–94.

    Article  CAS  PubMed  Google Scholar 

  10. Baguna J. The planarian neoblast: the rambling history of its origin and some current black boxes. Int J Dev Biol. 2012;56:19–37.

    Article  CAS  PubMed  Google Scholar 

  11. Scimone ML, Kravarik KM, Lapan SW, Reddien PW. Neoblast specialization in regeneration of the planarian Schmidtea mediterranea. Stem Cell Rep. 2014;3:339–52.

    Article  CAS  Google Scholar 

  12. Galliot B. Hydra, a fruitful model system for 270 years. Int J Dev Biol. 2012;56:411–23.

    Article  CAS  PubMed  Google Scholar 

  13. Holstein TW, Hobmayer E, Technau U. Cnidarians: an evolutionarily conserved model system for regeneration? Dev Dyn. 2003;226:257–67.

    Article  CAS  PubMed  Google Scholar 

  14. Bosch TC, Anton-Erxleben F, Hemmrich G, Khalturin K. The hydra polyp: nothing but an active stem cell community. Dev Growth Differ. 2010;52:15–25.

    Article  CAS  PubMed  Google Scholar 

  15. Technau U, Steele RE. Evolutionary crossroads in developmental biology: cnidaria. Development. 2011;138:1447–58.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  16. Bosch TC. Why polyps regenerate and we don’t: towards a cellular and molecular framework for hydra regeneration. Dev Biol. 2007;303:421–33.

    Article  CAS  PubMed  Google Scholar 

  17. Watanabe H, Hoang VT, Mattner R, Holstein TW. Immortality and the base of multicellular life: lessons from cnidarian stem cells. Semin Cell Dev Biol. 2009;20:1114–25.

    Article  CAS  PubMed  Google Scholar 

  18. Hobmayer E, Hatta M, Fischer R, Fujisawa T, Holstein TW, Sugiyama T. Identification of a Hydra homologue of the beta-catenin/plakoglobin/armadillo gene family. Gene. 1996;172:155–9.

    Article  CAS  PubMed  Google Scholar 

  19. Guder C, Philipp I, Lengfeld T, Watanabe H, Hobmayer B, Holstein TW. The Wnt code: cnidarians signal the way. Oncogene. 2006;25:7450–60.

    Article  CAS  PubMed  Google Scholar 

  20. Hobmayer B, Rentzsch F, Kuhn K, Happel CM, von Laue CC, Snyder P, Rothbacher U, Holstein TW. WNT signalling molecules act in axis formation in the diploblastic metazoan Hydra. Nature. 2000;407:186–9.

    Article  CAS  PubMed  Google Scholar 

  21. Holstein TW. Wnt signaling in cnidarians. Methods Mol Biol. 2008;469:47–54.

    Article  CAS  PubMed  Google Scholar 

  22. Lengfeld T, Watanabe H, Simakov O, Lindgens D, Gee L, Law L, Schmidt HA, Ozbek S, Bode H, Holstein TW. Multiple Wnts are involved in hydra organizer formation and regeneration. Dev Biol. 2009;330:186–99.

    Article  CAS  PubMed  Google Scholar 

  23. Kawakami Y, Rodriguez Esteban C, Raya M, Kawakami H, Marti M, Dubova I, Izpisua Belmonte JC. Wnt/beta-catenin signaling regulates vertebrate limb regeneration. Genes Dev. 2006;20:3232–7.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  24. Gurley KA, Rink JC, Sanchez Alvarado A. Beta-catenin defines head versus tail identity during planarian regeneration and homeostasis. Science. 2008;319:323–7.

    Article  CAS  PubMed  Google Scholar 

  25. Petersen CP, Reddien PW. Smed-betacatenin-1 is required for anteroposterior blastema polarity in planarian regeneration. Science. 2008;319:327–30.

    Article  CAS  PubMed  Google Scholar 

  26. Petersen CP, Reddien PW. Wnt signaling and the polarity of the primary body axis. Cell. 2009;139:1056–68.

    Article  CAS  PubMed  Google Scholar 

  27. Petersen CP, Reddien PW. A wound-induced Wnt expression program controls planarian regeneration polarity. Proc Natl Acad Sci U S A. 2009;106:17061–6.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  28. Gurley KA, Elliott SA, Simakov O, Schmidt HA, Holstein TW, Sanchez Alvarado A. Expression of secreted Wnt pathway components reveals unexpected complexity of the planarian amputation response. Dev Biol. 2010;347:24–39.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  29. Petersen CP, Reddien PW. Polarized notum activation at wounds inhibits Wnt function to promote planarian head regeneration. Science. 2011;332:852–5.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  30. Chera S, Ghila L, Dobretz K, Wenger Y, Bauer C, Buzgariu W, Martinou JC, Galliot B. Apoptotic cells provide an unexpected source of Wnt3 signaling to drive hydra head regeneration. Dev Cell. 2009;17:279–89.

    Article  CAS  PubMed  Google Scholar 

  31. Galliot B, Chera S. The Hydra model: disclosing an apoptosis-driven generator of Wnt-based regeneration. Trends Cell Biol. 2010;20:514–23.

    Article  CAS  PubMed  Google Scholar 

  32. Vriz S, Reiter S, Galliot B. Cell death: a program to regenerate. Curr Top Dev Biol. 2014;108:121–51.

    Article  CAS  PubMed  Google Scholar 

  33. Chera S, Ghila L, Wenger Y, Galliot B. Injury-induced activation of the MAPK/CREB pathway triggers apoptosis-induced compensatory proliferation in hydra head regeneration. Dev Growth Differ. 2011;53:186–201.

    Article  CAS  PubMed  Google Scholar 

  34. Galliot B. Injury-induced asymmetric cell death as a driving force for head regeneration in Hydra. Dev Genes Evol. 2013;223:39–52.

    Article  CAS  PubMed  Google Scholar 

  35. Wurtzel O, Cote LE, Poirier A, Satija R, Regev A, Reddien PW. A generic and cell-type-specific wound response precedes regeneration in planarians. Dev Cell. 2015;35:632–45.

    Article  CAS  PubMed  Google Scholar 

  36. Rink JC, Gurley KA, Elliott SA, Sanchez Alvarado A. Planarian Hh signaling regulates regeneration polarity and links Hh pathway evolution to cilia. Science. 2009;326:1406–10.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  37. Wenemoser D, Lapan SW, Wilkinson AW, Bell GW, Reddien PW. A molecular wound response program associated with regeneration initiation in planarians. Genes Dev. 2012;26:988–1002.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  38. Kao D, Felix D, Aboobaker A. The planarian regeneration transcriptome reveals a shared but temporally shifted regulatory program between opposing head and tail scenarios. BMC Genomics. 2013;14:797.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  39. Sandmann T, Vogg MC, Owlarn S, Boutros M, Bartscherer K. The head-regeneration transcriptome of the planarian Schmidtea mediterranea. Genome Biol. 2011;12:R76.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  40. Hand C, Uhlinger K. The culture, sexual and asexual reproduction and growth of the sea-anemone Nematostella vectensis. Biol Bull. 1992;182:169–76.

    Article  Google Scholar 

  41. Darling JA, Reitzel AR, Burton PM, Mazza ME, Ryan JF, Sullivan JC, Finnerty JR. Rising starlet: the starlet sea anemone, Nematostella vectensis. Bioessays. 2005;27:211–21.

    Article  CAS  PubMed  Google Scholar 

  42. Layden MJ, Rentzsch F, Rottinger E. The rise of the starlet sea anemone Nematostella vectensis as a model system to investigate development and regeneration. Wiley Interdiscip Rev Dev Biol. 2016;5(4):408–28.

    Article  PubMed  Google Scholar 

  43. Reitzel AM, Burton PM, Krone C, Finnerty JR. Comparison of developmental trajectories in the starlet sea anemone Nematostella vectensis: embryogenesis, regeneration, and two forms of asexual fission. Invertebr Biol. 2007;126:99–112.

    Article  Google Scholar 

  44. Putnam NH, Srivastava M, Hellsten U, Dirks B, Chapman J, Salamov A, Terry A, Shapiro H, Lindquist E, Kapitonov VV, et al. Sea anemone genome reveals ancestral eumetazoan gene repertoire and genomic organization. Science. 2007;317:86–94.

    Article  CAS  PubMed  Google Scholar 

  45. Finnerty JR, Martindale MQ. Homeoboxes in sea anemones (cnidaria:anthozoa): a PCR-based survey of nematostella vectensis and metridium senile. Biol Bull. 1997;193:62–76.

    Article  CAS  PubMed  Google Scholar 

  46. Finnerty JR, Martindale MQ. Ancient origins of axial patterning genes: Hox genes and ParaHox genes in the cnidaria. Evol Dev. 1999;1:16–23.

    Article  CAS  PubMed  Google Scholar 

  47. Finnerty JR, Paulson D, Burton P, Pang K, Martindale MQ. Early evolution of a homeobox gene: the parahox gene Gsx in the cnidaria and the bilateria. Evol Dev. 2003;5:331–45.

    Article  CAS  PubMed  Google Scholar 

  48. Chourrout D, Delsuc F, Chourrout P, Edvardsen RB, Rentzsch F, Renfer E, Jensen MF, Zhu B, de Jong P, Steele RE, Technau U. Minimal ProtoHox cluster inferred from bilaterian and cnidarian Hox complements. Nature. 2006;442:684–7.

    Article  CAS  PubMed  Google Scholar 

  49. Kamm K, Schierwater B. Ancient complexity of the non-Hox ANTP gene complement in the anthozoan nematostella vectensis: implications for the evolution of the ANTP superclass. J Exp Zoolog B Mol Dev Evol. 2006;306(6):589–96.

    Article  CAS  Google Scholar 

  50. Ryan JF, Mazza ME, Pang K, Matus DQ, Baxevanis AD, Martindale MQ, Finnerty JR. Pre-bilaterian origins of the Hox cluster and the Hox code: evidence from the sea anemone, Nematostella vectensis. PLoS One. 2007;2:e153.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  51. Sinigaglia C, Busengdal H, Leclere L, Technau U, Rentzsch F. The bilaterian head patterning gene six3/6 controls aboral domain development in a cnidarian. PLoS Biol. 2013;11:e1001488.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  52. Bossert PE, Dunn MP, Thomsen GH. A staging system for the regeneration of a polyp from the aboral physa of the anthozoan cnidarian nematostella vectensis. Dev Dyn. 2013;242:C1.

    Article  Google Scholar 

  53. Amiel AR, Johnston HT, Nedoncelle K, Warner JF, Ferreira S, Rottinger E. Characterization of morphological and cellular events underlying oral regeneration in the Sea anemone, Nematostella vectensis. Int J Mol Sci. 2015;16:28449–71.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  54. Singer II. Tentacular and oral-disc regeneration in the sea anemone, aiptasia diaphana. 3. Autoradiographic analysis of patterns of tritiated thymidine uptake. J Embryol Exp Morphol. 1971;26:253–70.

    CAS  PubMed  Google Scholar 

  55. Passamaneck YJ, Martindale MQ. Cell proliferation is necessary for the regeneration of oral structures in the anthozoan cnidarian Nematostella vectensis. BMC Dev Biol. 2012;12:34.

    Article  PubMed  PubMed Central  Google Scholar 

  56. Trevino M, Stefanik DJ, Rodriguez R, Harmon S, Burton PM. Induction of canonical Wnt signaling by alsterpaullone is sufficient for oral tissue fate during regeneration and embryogenesis in Nematostella vectensis. Dev Dyn. 2011;240:2673–9.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  57. DuBuc TQ, Traylor-Knowles N, Martindale MQ. Initiating a regenerative response; cellular and molecular features of wound healing in the cnidarian Nematostella vectensis. BMC Biol. 2014;12:24.

    Article  PubMed  PubMed Central  Google Scholar 

  58. Sun G, Irvine KD. Control of growth during regeneration. Curr Top Dev Biol. 2014;108:95–120.

    Article  CAS  PubMed  Google Scholar 

  59. Anders S, Huber W. Differential expression analysis for sequence count data. Genome Biol. 2010;11:R106.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  60. Adell T, Salo E, Boutros M, Bartscherer K. Smed-Evi/Wntless is required for beta-catenin-dependent and -independent processes during planarian regeneration. Development. 2009;136:905–10.

    Article  CAS  PubMed  Google Scholar 

  61. Gill SE, Parks WC. Metalloproteinases and their inhibitors: regulators of wound healing. Int J Biochem Cell Biol. 2008;40:1334–47.

    Article  CAS  PubMed  Google Scholar 

  62. Kelwick R, Desanlis I, Wheeler GN, Edwards DR. The ADAMTS (a disintegrin and metalloproteinase with thrombospondin motifs) family. Genome Biol. 2015;16:113.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  63. Porter S, Clark IM, Kevorkian L, Edwards DR. The ADAMTS metalloproteinases. Biochem J. 2005;386:15–27.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  64. Sarras Jr MP, Yan L, Leontovich A, Zhang JS. Structure, expression, and developmental function of early divergent forms of metalloproteinases in hydra. Cell Res. 2002;12:163–76.

  65. Carroll SB, Grenier JK, Weatherbee SD. From DNA to Diversity: Molecular Genetics and the Evolution of Animal Design. Malden: Wiley-Blackwell; 2004.

  66. Larroux C, Fahey B, Liubicich D, Hinman VF, Gauthier M, Gongora M, Green K, Worheide G, Leys SP, Degnan BM. Developmental expression of transcription factor genes in a demosponge: insights into the origin of metazoan multicellularity. Evol Dev. 2006;8:150–73.

    Article  CAS  PubMed  Google Scholar 

  67. Lagutin OV, Zhu CC, Kobayashi D, Topczewski J, Shimamura K, Puelles L, Russell HR, McKinnon PJ, Solnica-Krezel L, Oliver G. Six3 repression of Wnt signaling in the anterior neuroectoderm is essential for vertebrate forebrain development. Genes Dev. 2003;17:368–79.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  68. Petersen HO, Hoger SK, Looso M, Lengfeld T, Kuhn A, Warnken U, Nishimiya-Fujisawa C, Schnolzer M, Kruger M, Ozbek S, et al. A comprehensive transcriptomic and proteomic analysis of hydra head regeneration. Mol Biol Evol. 2015;32:1928–47.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  69. Bradshaw B, Thompson K, Frank U. Distinct mechanisms underlie oral vs aboral regeneration in the cnidarian Hydractinia echinata. Elife. 2015;4:e05506.

    PubMed  PubMed Central  Google Scholar 

  70. Hilman D, Gat U. The evolutionary history of YAP and the hippo/YAP pathway. Mol Biol Evol. 2011;28:2403–17.

    Article  CAS  PubMed  Google Scholar 

  71. Sopko R, McNeill H. The skinny on Fat: an enormous cadherin that regulates cell adhesion, tissue growth, and planar cell polarity. Curr Opin Cell Biol. 2009;21:717–23.

    Article  CAS  PubMed  Google Scholar 

  72. Hensel K, Lotan T, Sanders SM, Cartwright P, Frank U. Lineage-specific evolution of cnidarian Wnt ligands. Evol Dev. 2014;16:259–69.

    Article  CAS  PubMed  Google Scholar 

  73. Kusserow A, Pang K, Sturm C, Hrouda M, Lentfer J, Schmidt HA, Technau U, von Haeseler A, Hobmayer B, Martindale MQ, Holstein TW. Unexpected complexity of the Wnt gene family in a sea anemone. Nature. 2005;433:156–60.

    Article  CAS  PubMed  Google Scholar 

  74. Rigo-Watermeier T, Kraft B, Ritthaler M, Wallkamm V, Holstein T, Wedlich D. Functional conservation of Nematostella Wnts in canonical and noncanonical Wnt-signaling. Biol Open. 2012;1:43–51.

    Article  CAS  PubMed  Google Scholar 

  75. Philipp I, Aufschnaiter R, Ozbek S, Pontasch S, Jenewein M, Watanabe H, Rentzsch F, Holstein TW, Hobmayer B. Wnt/beta-catenin and noncanonical Wnt signaling interact in tissue evagination in the simple eumetazoan Hydra. Proc Natl Acad Sci U S A. 2009;106:4290–5.

    Article  PubMed  PubMed Central  Google Scholar 

  76. Lee PN, Pang K, Matus DQ, Martindale MQ. A WNT of things to come: evolution of Wnt signaling and polarity in cnidarians. Semin Cell Dev Biol. 2006;17:157–67.

    Article  PubMed  CAS  Google Scholar 

  77. Guder C, Pinho S, Nacak TG, Schmidt HA, Hobmayer B, Niehrs C, Holstein TW. An ancient Wnt-dickkopf antagonism in hydra. Development. 2006;133:901–11.

    Article  CAS  PubMed  Google Scholar 

  78. Takahashi M, Fujita M, Furukawa Y, Hamamoto R, Shimokawa T, Miwa N, Ogawa M, Nakamura Y. Isolation of a novel human gene, APCDD1, as a direct target of the beta-Catenin/T-cell factor 4 complex with probable involvement in colorectal carcinogenesis. Cancer Res. 2002;62:5651–6.

    CAS  PubMed  Google Scholar 

  79. Shimomura Y, Agalliu D, Vonica A, Luria V, Wajid M, Baumer A, Belli S, Petukhova L, Schinzel A, Brivanlou AH, et al. APCDD1 is a novel Wnt inhibitor mutated in hereditary hypotrichosis simplex. Nature. 2010;464:1043–7.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  80. Nakamura Y, Tsiairis CD, Ozbek S, Holstein TW. Autoregulatory and repressive inputs localize Hydra Wnt3 to the head organizer. Proc Natl Acad Sci U S A. 2011;108:9137–42.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  81. Moran Y, Praher D, Schlesinger A, Ayalon A, Tal Y, Technau U. Analysis of soluble protein contents from the nematocysts of a model sea anemone sheds light on venom evolution. Mar Biotechnol (NY). 2013;15:329–39.

    Article  CAS  Google Scholar 

  82. Vadon-Le Goff S, Hulmes DJ, Moali C. BMP-1/tolloid-like proteinases synchronize matrix assembly with growth factor activation to promote morphogenesis and tissue remodeling. Matrix Biol. 2015;44–46:14–23.

    Article  PubMed  CAS  Google Scholar 

  83. Lin J, Patel SR, Cheng X, Cho EA, Levitan I, Ullenbruch M, Phan SH, Park JM, Dressler GR. Kielin/chordin-like protein, a novel enhancer of BMP signaling, attenuates renal fibrotic disease. Nat Med. 2005;11:387–93.

    Article  CAS  PubMed  Google Scholar 

  84. Klug M, Tardent P, Smid I, Holstein T. Presence and localization of chitinase in hydra and podocoryne (cnidaria, hydrozoa). J Exp Zool. 1984;229:69–72.

    Article  CAS  Google Scholar 

  85. Mali B, Mohrlen F, Frohme M, Frank U. A putative double role of a chitinase in a cnidarian: pattern formation and immunity. Dev Comp Immunol. 2004;28:973–81.

    Article  CAS  PubMed  Google Scholar 

  86. Dunn DF, Liberman MH. Chitin in sea anemone shells. Science. 1983;221:157–9.

    Article  CAS  PubMed  Google Scholar 

  87. Garnham CP, Roll-Mecak A. The chemical complexity of cellular microtubules: tubulin post-translational modification enzymes and their roles in tuning microtubule functions. Cytoskeleton (Hoboken). 2012;69:442–63.

    Article  CAS  Google Scholar 

  88. Nachury MV, Loktev AV, Zhang Q, Westlake CJ, Peranen J, Merdes A, Slusarski DC, Scheller RH, Bazan JF, Sheffield VC, Jackson PK. A core complex of BBS proteins cooperates with the GTPase Rab8 to promote ciliary membrane biogenesis. Cell. 2007;129:1201–13.

    Article  CAS  PubMed  Google Scholar 

  89. Goetz SC, Anderson KV. The primary cilium: a signalling centre during vertebrate development. Nat Rev Genet. 2010;11:331–44.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  90. Lancaster MA, Schroth J, Gleeson JG. Subcellular spatial regulation of canonical Wnt signalling at the primary cilium. Nat Cell Biol. 2011;13:700–7.

    PubMed  PubMed Central  Google Scholar 

  91. May-Simera HL, Kelley MW. Cilia, Wnt signaling, and the cytoskeleton. Cilia. 2012;1:7.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  92. Reddien PW. Constitutive gene expression and the specification of tissue identity in adult planarian biology. Trends Genet. 2011;27:277–85.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  93. Rentzsch F, Fritzenwanker JH, Scholz CB, Technau U. FGF signalling controls formation of the apical sensory organ in the cnidarian Nematostella vectensis. Development. 2008;135:1761–9.

    Article  CAS  PubMed  Google Scholar 

  94. Fritzenwanker JH, Saina M, Technau U. Analysis of forkhead and snail expression reveals epithelial-mesenchymal transitions during embryonic and larval development of Nematostella vectensis. Dev Biol. 2004;275:389–402.

    Article  CAS  PubMed  Google Scholar 

  95. Magie CR, Pang K, Martindale MQ. Genomic inventory and expression of Sox and Fox genes in the cnidarian nematostella vectensis. Dev Genes Evol. 2005;215:618–30.

    Article  CAS  PubMed  Google Scholar 

  96. Koinuma S, Umesono Y, Watanabe K, Agata K. Planaria FoxA (HNF3) homologue is specifically expressed in the pharynx-forming cells. Gene. 2000;259:171–6.

    Article  CAS  PubMed  Google Scholar 

  97. Adler CE, Seidel CW, McKinney SA, Sanchez Alvarado A. Selective amputation of the pharynx identifies a FoxA-dependent regeneration program in planaria. Elife. 2014;3:e02238.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  98. Srivastava M, Mazza-Curll KL, van Wolfswinkel JC, Reddien PW. Whole-body acoel regeneration is controlled by Wnt and Bmp-Admp signaling. Curr Biol. 2014;24:1107–13.

    Article  CAS  PubMed  Google Scholar 

  99. Stefanik DJ, Wolenski FS, Friedman LE, Gilmore TD, Finnerty JR. Isolation of DNA, RNA and protein from the starlet sea anemone Nematostella vectensis. Nat Protoc. 2013;8:892–9.

    Article  PubMed  CAS  Google Scholar 

  100. Dobin A, Davis CA, Schlesinger F, Drenkow J, Zaleski C, Jha S, Batut P, Chaisson M, Gingeras TR. STAR: ultrafast universal RNA-seq aligner. Bioinformatics. 2013;29:15–21.

    Article  CAS  PubMed  Google Scholar 

  101. Anders S, Huber W. Differential expression of RNA-Seq data at the gene level – the DESeq package. 2016. p. 1–23.

    Google Scholar 

  102. Ernst J, Bar-Joseph Z. STEM: a tool for the analysis of short time series gene expression data. BMC Bioinf. 2006;7:191.

    Article  CAS  Google Scholar 

  103. Lawrence IKL. A concordance correlation coefficient to evaluate reproducibility. Biometrics. 1989;45:255–68.

    Article  Google Scholar 

  104. da Huang W, Sherman BT, Lempicki RA. Systematic and integrative analysis of large gene lists using DAVID bioinformatics resources. Nat Protoc. 2009;4:44–57.

    Article  CAS  Google Scholar 

  105. Huang DW, Sherman BT, Tan Q, Collins JR, Alvord WG, Roayaei J, Stephens R, Baseler MW, Lane HC, Lempicki RA. The DAVID gene functional classification tool: a novel biological module-centric algorithm to functionally analyze large gene lists. Genome Biol. 2007;8:R183.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  106. Huang DW, Sherman BT, Tan Q, Kir J, Liu D, Bryant D, Guo Y, Stephens R, Baseler MW, Lane HC, Lempicki RA. DAVID bioinformatics resources: expanded annotation database and novel algorithms to better extract biology from large gene lists. Nucleic Acids Res. 2007;35:W169–175.

    Article  PubMed  PubMed Central  Google Scholar 

Download references

Acknowledgements

We thank Inbar Plaschkes for her help in RNA-Seq results analysis, Adi Faigenboim for help with gene clustering and Yehu Moran and Arthur Schaffer for critical reading of the paper. This research was funded by the Ministry of Science and Technology grant # 3-8774 to UG.

Availability of data and materials

The RNA-Seq data of this study including the accession numbers (Gene i.d.) of the DE genes, the adjusted time effect p-values, the normalized and raw reads are given in Additional file 1: Table S1. The transcriptome sequencing data has been submitted to NCBI’s Sequence Read Archive (SRA) with accession number SRP079101.

Authors’ contributions

UG and AS conceived and designed the RNA-Seq screen. AS performed the experimental procedures including the construction of the Illumina libraries and qPCR with the participation of MB. VC and AS performed the bioinformatics analyses. KL performed some of the qPCR studies. UG, AS and MB analyzed the results. UG and AS wrote the manuscript and all authors read and approved the contents.

Competing interests

The authors declare that they have no competing interests.

Consent for publication

Not applicable.

Ethics approval and consent to participate

Not applicable.

Author information

Authors and Affiliations

Authors

Corresponding author

Correspondence to Uri Gat.

Additional files

Additional file 1: Table S1.

The differentially expressed (DE) genes of the regeneration screen. The list contains the gene id, followed by Uniprot or Blast2GO annotation and FDR adjusted p-value of the time effect. This is followed by the average reads for each time point, DESeq normalized read count per each duplicate time point and raw read count per each duplicate time point. The annotation (−−) appears when Blast2GO did not retrieve any blast hit in the RefSeq protein database. (XLSX 2933 kb)

Additional file 2: Figure S1.

Clusters of DE oral and physal genes. List of all 20 clusters of oral and aboral expression produced by STEM (see Methods). Clusters are ordered by number of genes in each cluster. Numbers in the top left corner represent the cluster number. Numbers in the bottom left corner represent the number of genes in each cluster. A list of DE genes and the cluster they belong to can be found in Additional file 3. (PDF 213 kb)

Additional file 3: Table S2.

DE genes and their CCC. List of all DE genes and their time effect FDR adjusted p-value and CCC. In addition, the side of higher expression in most of the time points and the side and hour of the highest level of change in expression per gene are given. The cluster the gene belongs to (from Fig. 2c and Additional file 2: Figure S1) is also indicated, as is whether the gene is consistently expressed throughout regeneration (Fig. 3a) (XLSX 524 kb)

Additional file 4: Table S3.

Enriched GO terms in low CCC genes. GO terms containing similar genes are colored similarly. Additional tabs contain a list of GO annotated low CCC genes belonging to the GO terms mentioned in Fig. 4 and a list of all the corresponding GO terms associated with each gene. (XLSX 575 kb)

Additional file 5: Table S4.

Metalloproteinase (MP) groups and inhibitors. A list of the major DE MPs identified in the screen with ADAMTS, Astacins, other MPs and their inhibitors. (XLSX 676 kb)

Additional file 6: Figure S3.

Time course graphs of metalloproteinase (MP) group genes and inhibitors. The FDR-adjusted p-values for the time effect (DESeq GLMs) for all genes are given. (XLSX 23 kb)

Additional file 7: Table S5.

Homeobox and other TF groups that are highly DE. (XLSX 790 kb)

Additional file 8: Table S6.

Wnt family genes clusters. A list of the cluster groups of the Wnt family genes and a list of the genes in each cluster shown in Fig. 5a. (XLSX 11 kb)

Additional file 9: Figure S4.

Time course graphs of the major DE Wnt factors and Wnt pathway members. (XLSX 81 kb)

Additional file 10: Table S7.

DE genes of different developmental pathways. A list of DE genes of major the developmental pathways that were manually detected (see Methods) and have a CCC < 0.6. The time-course graphs of prominent pathway genes are also shown. (XLSX 625 kb)

Additional file 11: Figure S5.

Expression of chitin upon regeneration. A. Expression patterns of representative chitin genes in our screen showing the difference in their oral and aboral expression levels at the hour 0 time point. B. A scheme illustrating the possible gradient patterns in intact Nematostella as suggested by qPCR on the different parts. These gradients can explain the early difference in expression level and the later return to the original polarized expression levels. (PDF 1722 kb)

Additional file 12: Table S8.

Time course graphs of the major DE tubulin and cilia associated genes. (XLSX 40 kb)

Additional file 13: Table S9.

Comparison of polarized gene expression in regeneration between Nematostella and Planaria. The anterior and posterior enriched genes of Planaria are listed with their oral or aboral Nematostella orthologs and homologs when present. The color coding is the same as in Fig. 7. (XLSX 17 kb)

Additional file 14: Table S10.

A list of the primers that were used for qRT-PCR experiments. (XLSX 11 kb)

Rights and permissions

Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated.

Reprints and permissions

About this article

Check for updates. Verify currency and authenticity via CrossMark

Cite this article

Schaffer, A.A., Bazarsky, M., Levy, K. et al. A transcriptional time-course analysis of oral vs. aboral whole-body regeneration in the Sea anemone Nematostella vectensis . BMC Genomics 17, 718 (2016). https://doi.org/10.1186/s12864-016-3027-1

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: https://doi.org/10.1186/s12864-016-3027-1

Keywords