A Comprehensive Transcriptomic and Proteomic Analysis of Hydra Head Regeneration.

The cnidarian freshwater polyp Hydra sp. exhibits an unparalleled regeneration capacity in the animal kingdom. Using an integrative transcriptomic and stable isotope labeling by amino acids in cell culture proteomic/phosphoproteomic approach, we studied stem cell-based regeneration in Hydra polyps. As major contributors to head regeneration, we identified diverse signaling pathways adopted for the regeneration response as well as enriched novel genes. Our global analysis reveals two distinct molecular cascades: an early injury response and a subsequent, signaling driven patterning of the regenerating tissue. A key factor of the initial injury response is a general stabilization of proteins and a net upregulation of transcripts, which is followed by a subsequent activation cascade of signaling molecules including Wnts and transforming growth factor (TGF) beta-related factors. We observed moderate overlap between the factors contributing to proteomic and transcriptomic responses suggesting a decoupled regulation between the transcriptional and translational levels. Our data also indicate that interstitial stem cells and their derivatives (e.g., neurons) have no major role in Hydra head regeneration. Remarkably, we found an enrichment of evolutionarily more recent genes in the early regeneration response, whereas conserved genes are more enriched in the late phase. In addition, genes specific to the early injury response were enriched in transposon insertions. Genetic dynamicity and taxon-specific factors might therefore play a hitherto underestimated role in Hydra regeneration.

Introduction

Different from plants, many animals have lost the ability to regenerate. At present, however, any comprehensive picture of the developmental mechanisms that are required to restore a missing or damaged tissue during animal regeneration is lacking (Gilbert 2014). It is even unclear whether different animal systems have evolved different strategies to restore missing structures (Brockes and Kumar 2008). Here, we studied head regeneration in the freshwater polyp Hydra sp. (Hydrozoa, Cnidaria), which is a paradigm for animal regeneration.

Trembley (1744) discovered the phenomenon of animal regeneration when he was cutting a Hydra polyp into two halves and observed each of them regenerating an intact individual within 2 days (Trembley 1744; Bode 2003; Holstein et al. 2003; Bosch et al. 2010). To initiate regeneration in Hydra, an injury signal is essential (Newman 1974; Kobatake and Sugiyama 1989; Guder et al. 2006). A first response on the cellular level is a significant cell cycle dynamics with a downregulation of mitosis followed by a later upregulation (Holstein et al. 1991) even though regeneration can also start without mitosis by a repatterning of existing tissue (“morphallaxis”; Hicklin and Wolpert 1973; Dübel and Schaller 1990). Also on the cellular level, a rapid reorganization of the epithelial cell layers is necessary for an effective wound closure. This process includes the rapid formation of new septate and gap junctions already 1 h after head removal (Bibb and Campbell 1973; Wood and Kuda 1980a, 1980b). The tissue reorganization later involves also the synthesis of a new mesoglea separating the ectodermal and endodermal tissue layers (Sarras et al. 1993). On the molecular level, Wnt signaling was described as major pathway acting throughout the entire regeneration process (Hobmayer et al. 2000; Technau et al. 2000; Guder et al. 2006; Chera et al. 2009; Lengfeld et al. 2009). The contribution of other pathways like transforming growth factor (TGF) beta signaling (Reinhardt et al. 2004; Rentzsch et al. 2007) is not well understood and so far a comprehensive molecular understanding of Hydra regeneration is lacking.

Recent whole genome and transcriptome studies of Hydra have revealed that almost the complete molecular toolkit is shared by bilaterian and cnidarian genomes (Chapman et al. 2010; Steele et al. 2011; Hemmrich et al. 2012; Krishna et al. 2013; Wenger and Galliot 2013; Juliano et al. 2014). In a recent review, a complex immune response was described on the transcriptome level for injured Hydra (Wenger et al. 2014). Similarly, in the starlet sea anemone Nematostella vectensis, an upregulation of immune response genes has been observed after injury of adult polyps (Dubuc et al. 2014).

The regeneration capacity varies largely among cnidarians. In N. vectensis, adult polyps can regenerate, but their regeneration capacity is limited. For example, it is not possible to dissociate intact N. vectensis polyps into single cells and restore intact new animals from such cell suspensions. And in many jellyfish, the ability to regenerate is completely absent (Tardent 1978). By comparison, when Hydra is dissociated into single cells, reaggregates of dissociated cells regenerate polyps size dependent within several days. Hence, Hydra also serves as a model for studying regeneration under conditions of de novo pattern formation without any preexisting polarity of the tissue (Gierer et al. 1972; Technau et al. 2000). However, the underlying molecular mechanisms of this variable regeneration capacity in cnidarians are completely unknown so far.

Here, we performed a comprehensive proteome and transcriptome analysis and established a large data set (available on http://hydra.cos.uni-heidelberg.de genome browser, last accessed April 3, 2015), which should be suited to monitor immediate, early and late responses in the regeneration process of Hydra. We performed a RNAseq study of head regenerating tissue, and developed a metabolic labeling approach for quantitative proteomics in Hydra by using stable isotope labeling by amino acids in cell culture (SILAC; Ong et al. 2002) which has the potential to identify small changes in protein abundance to be expected in processes like regeneration (Boser et al. 2013). By integrating different proteomic and transcriptomic data sets using a recently proposed multiple coinertia method (Meng et al. 2014), we identified as the two most significant response patterns an early injury response and a subsequent, signaling driven patterning of the regenerating tissue. The major contributors to injury are responses in global scale metabolic networks including Erk (MAPK), c-Jun N-terminal kinase (JNK), and Gsk3, whereas Wnt signaling seems to orchestrate the downstream signaling pathways in the patterning response, which represents an evolutionary conserved regulatory network in animal regeneration. Our data are underpinned on the functional level by an additional phosphoproteomic analysis of the regeneration process.

In our phylogenomic analysis, we also made the discovery that genes involved in the early phase of Hydra regeneration are evolutionarily more recent and their flanking regions are more dynamic than those involved in the patterning response. This analysis is based on the presence and age of repetitive sequences in the very dynamic genome of Hydra (Chapman et al. 2010) as well as the concept of an age index (Domazet-Lošo and Tautz 2010) applied to regeneration-specific genes. Transposable element insertions can act either on the present transcriptional units (resulting in deleterious or neutral mutations) or contribute to the origin of new regulatory elements (e.g., Chenais et al. 2012). We propose that estimating the transposable element burden in the genes active during the regeneration process may provide a new clue about the origin and evolutionary dynamics of regeneration.

Results

SILAC-Proteome and RNAseq-Transcriptome of Hydra Head Regenerates

Quantitative proteome analyses provide the most accurate and earliest information about a biological system because the measurements directly focus on the biological effector molecules (Gygi et al. 1999). We therefore combined this approach with RNAseq transcriptome profiling for a complete picture of the molecular dynamics of Hydra head regeneration (fig. 1).

F

Quantitative approach to study Hydra regeneration. (A) SILAC labeling. Yeast Lys deficient strain was labeled with [13C6]lysine and fed to Moina macrocopa copepods, which served as a freshwater food source for Hydra polyps. (B) Head regeneration was initiated by decapitation and regenerating tip tissue was isolated at the times indicated. (C) SILAC of an intact Hydra culture was used as “spike-in” standard to compare different regeneration samples. Parallel samples were used for the proteome analysis and for the RNAseq transcriptome analysis. See Materials and Methods and text for further details.

For the proteome approach, we adopted SILAC metabolic labeling to Hydra (fig. 1A). Production of SILAC Hydra required the establishment of a food chain, where labeled yeast was fed to freshwater cladoceran (Moina macropora) cultured in medium enriched with highly unsaturated fatty acids (HUFA). Moina macropora served as efficient food source for Hydra growth (fig. 1A). By this approach, we obtained a labeling efficiency of greater than 90%, which is in the range of other non-cell culture systems (Geiger et al. 2011).

To ensure a short regeneration time, we cut the animals at 80% body length (fig. 1B). Under these conditions gastric tissue can form a new head within 36–48 h (MacWilliams 1983). We isolated regenerating tips from 0, 0.5, 3, 6, 12, 24, and 48 h with three biological replicates each (fig. 1C). Twelve hours after head removal regenerating tips have the full morphogenetic potential of intact heads and at later stages the prepatterned regenerating tissue simply develops into a fully differentiated head (MacWilliams 1983). Therefore, the proteome study was restricted to the first 12 h after head removal (fig. 1C).

To quantify protein expression in regenerating animals, we used the SILAC proteome from intact and regenerating Hydra polyps that served as the heavy Lys6-labeled reference (fig. 1C, blue bar) as well as unlabeled regenerating tip tissue. After cell lysis, equal amounts of labeled and unlabeled proteins of regenerating tips at a given time point were mixed (fig. 1C, green bar). The peptide ratios from labeled and unlabeled regenerating tissue were then calculated resulting in a quantification of 3,410 peptide predictions from the proteome data set (fig. 1C; see also Materials and Methods). The SILAC labeling was also used to identify the dynamics of protein phosphorylation and dephosphorylation during regeneration.

In the transcriptome analysis, all raw reads were quality assessed and reads from time point 0–12 h were used to assemble a de novo transcriptome, which consists of 36,338 high quality transcripts (Materials and Methods). To select the best reference for the differential expression analysis, we mapped all transcriptomic reads to the Hydra genome (Chapman et al. 2010). Mapping rates were significantly lower when compared with the de novo transcriptome assembly (50–60% vs. 80%, respectively). The drop in coverage is most likely due to the fragmented genome assembly, but also polymorphism. We therefore have chosen our transcriptome assembly as reference for all subsequent analyses. We have used open reading frame (ORF) prediction from EMBOSS package (Chen et al. 2005) accompanied by the longest ORF selection for each transcript based on the similarity to the Hydra proteome (Chapman et al. 2010). This resulted in 27,913 peptide predictions, which were also used for peptide quantification (fig. 1C).

Our quantitative data on the transcriptome, proteome, and phosphoproteome are summarized in supplementary tables S1–S8, Supplementary Material online.

Decoupling between Transcriptomic and Proteomic Regeneration Responses

In order to detect the earliest regeneration responses on the proteome level and to disentangle the subsequent dynamics between the proteome and transcriptome levels we combined both data sets. So far only few approaches to study molecular changes at both proteomic and transcriptomic levels exist (Zhang et al. 2010; Looso et al. 2013). To reveal these dynamics in an unbiased approach, we applied the recently established multiple coinertia analysis (MCIA; Meng et al. 2014; fig. 2). Similar to a principal component analysis, MCIA is able to identify the most highly correlated patterns and their contribution to explain the total variance, but for several data sets simultaneously. We thus analyzed only those contigs that were shared between both data sets. Similar transcripts (mostly allelic variants) were merged based on overlapping peptide hits. This resulted in a total of 3,410 contigs.

F

Correlation analysis of proteome and transcriptome. (A) Multiple coinertia representation of the data, x axis and y axis are the first and second principal components, respectively. (B) Contribution of transcripts and (C) peptides to each axis, loci contributing significantly to each axis are color-coded. (D) Venn diagram shows little overlap among the different contributing loci, color-code corresponds to panel B. (E) Heat map showing the dynamics of both peptide and transcript levels for each locus (rows) and at different time points (columns). Only the most significant genes are represented. Heat map colors: green: higher amount; magenta: lower amount.

The four retained axes of the MCIA had the following variance distribution: first axis explains 65%, second axis 14% (the subsequent axes explain 11% and 10%, respectively). Figure 2A shows that axis 1 represents all time points except 0.5 h, whereas axis 2 has the highest contribution from 0.5 h. We presume that axis 2 corresponds to the injury response and initial wound healing, whereas the “general” regeneration response corresponds to axis 1. This result provides support for the notion of two different molecular programs acting during immediate and early regeneration responses.

We then determined the genes that significantly contributed to either axis 1 or 2 by computing the participation ratio, which provides an estimate for the number of variables (genes) that are sufficient to describe the main pattern within the data set. We selected those genes that contributed positively to axes 1 and 2 (fig. 2B and C and table 1; supplementary tables S1–S4, Supplementary Material online). At the transcriptome level (fig. 2B), the early response (156 upregulated genes) was more pronounced than the regeneration response (42 upregulated genes). This pattern appears reversed for the peptides (fig. 2C), as we observe only 29 peptides specific to the injury response and 385 during regeneration (supplementary table S12, Supplementary Material online).

Table 1.

MCIA Contribution Analysis.

MF	Prot Inj	Prot Reg	Tran Inj	Tran Reg	All
MF00001 receptor	0	3	2	5	560
MF00008 ligand-gated ion channel	0	1	0	0	51
MF00016 signaling molecule	0	1	1	0	399
MF00036 transcription factor	0	4	0	0	547
MF00040 cell adhesion molecule	1	1	1	0	210
MF00042 nucleic acid binding	5	45	20	4	2405
MF00048 DNA ligase	0	1	0	0	42
MF00054 RNA methyltransferase	0	1	0	0	18
MF00077 chaperone	0	2	9***	2	138
MF00082 transporter	0	7	2	1	499
MF00087 transfer/carrier protein	0	3	4*	0	87
MF00091 cytoskeletal protein	2	33***	11	5	754
MF00093 select regulatory molecule	4	32***	6	2	721
MF00107 kinase	0	14	2	0	401
MF00113 phosphatase	0	2	0	0	243
MF00118 synthase and synthetase	1	1	1	1	121
MF00123 oxidoreductase	0	2	4	1	465
MF00131 transferase	0	7	6	0	679
MF00141 hydrolase	0	5	0	0	567
MF00153 protease	0	3	6	1	400
MF00157 lyase	0	1	0	0	156
MF00166 Isomerase	0	3	4	0	106
MF00170 ligase	0	7	2	0	218
MF00173 defense/immunity protein	0	3	1	0	83
MF00178 extracellular matrix	0	0	1	2	144
MF00188 select calcium binding protein	1	2	2	1	110
MF00197 miscellaneous function	2	17	5	5	639
MF00267 membrane traffic protein	0	17**	2	0	342

Note.—MF, molecular function.

aThe enrichment of different annotation classes in axis 1 (regeneration; Reg) and axis 2 (injury; Inj) for the proteome (Prot) and transcriptome (Tran) data sets where analyzed by using the Fisher’s exact test and Bonferroni multiple test correction by comparing the counts of a given annotation within the group to the total counts of a given annotation in the rest of the loci. This resulted in enriched categories defined by PANTHER (Protein ANalysis THrough Evolutionary Relationships) Molecular Function (MF) categories. Significance levels indicate different P-values (*<0.1; **<0.01; ***<0.0001).

The Venn diagram shows that there was very little overlap between the genes contributing most to either axis in both data sets (fig. 2D). This indicates a complex molecular response pattern during regeneration, not necessarily involving the same genes in both proteome and transcriptome data sets. We then sought to investigate the dynamics of those transcripts and peptides. For this, the expression and peptide counts were plotted on a heat map (fig. 2E). The four significant groups were used to define the order of the rows (genes), therefore no heat map specific clustering was done. This method shows the common patterns within each group. Interestingly, the 42 regeneration transcripts (fig. 2E yellow, right side) show gradual upregulation over time. By comparison, the group of 385 regeneration peptides (fig. 2E blue, left side) shows a high contribution at both, 0.5 and 12 h, suggesting their involvement in both early response and prepatterning.

In conclusion, comparing the proteomic and transcriptomic data sets, we find distinct and independent levels of regulation during the different phases of hydra head regeneration.

Regeneration-Specific Arrest to Cell Cycle and Cell Signaling

Of special interest were early injury proteome and transcriptome responses (axis 2). We could identify 29 proteins that were strongly upregulated within the first hour after head removal (fig. 2E, red, left side; tables 1 and 2; supplementary table S1, Supplementary Material online). This group was almost exclusively comprised by proteins acting in cell cycle control and cell survival, and included several stress response factors. Among the cell cycle control proteins we found Rif1, a major cell cycle arrest factor that controls DNA replication by interfering with the phosphorylation of the minichromosome maintenance (MCM) complex (Hiraga et al. 2014; Mattarocci et al. 2014). We also found a transient downregulation of Cyclin D up to 6 h after head removal. Besides the cell cycle arrest factors, we also found Fas-associated factor Faf1, which is involved in inhibition of Wnt-signaling by β-catenin degradation (Zhang et al. 2011) and in apoptosis (Chu et al. 1995). Other proteins include GTPase activating protein, and vacuolar sorting-associated protein, indicating involvement of endocytosis (Lamb et al. 2013), Grasp involved in endosome recycling (Hoogenraad et al. 2010), and Dicer involved in miRNA processing (Botchkareva 2012). For the early injury transcriptome, we found mostly stress response genes to be enriched (table 1; supplementary tables S3 and S6, Supplementary Material online). This category includes, for example, Hsp90 but also 14-3-3 proteins (a known activator of protein kinase A; Kent et al. 2010) and a glucose-regulated protein (Grp78/94; Kozutsumi et al. 1988).

Table 2.

List of Injury Proteome Response Factors and their Putative Functions.

Injury Proteomea	Putative Function	References
Injury response
Chymotrypsin-like elastase	Enzyme	McDonnell et al. 1999
E3 Ubiquitin ligase 1, 2, 3	Protein ubiquitination	Berndsen and Wolberger 2014
Periostin	Injury response	Conway et al. 2014
Grip1 associated protein 1 (GRASP1)	endosome recycling	Hoogenraad et al. 2010
TBC-1	Silencing neuro. progenitors	Falace et al. 2014
Cingulin	Tight junction	Citi et al. 2014
Trichohyalin	Wound closure	Steinert et al. 2003
RasGEF	Membrane traffic	Mizuno-Yamasak et al. 2012
Cell cycle control, cell death, and survival
GTPase	IMAP member, membrane remodeling	Van Aelst and Symons 2002
Rif1	Cell cycle arrest	Mattarocci et al. 2014
Smc	Condensin complex, mitosis	Thadani et al. 2012
vWillebrand/Kazal	ECM	Ozbek et al. 2010
Dicer1	Cleavage of dsRNA	Foulkes et al. 2014
Nucleic acid binding
Complexed with Cdc5 protein Cwf19	Splice factor	Ohi et al. 2002
Tudor-domain	DNA repair	Lu and Wang 2013
RNA-binding protein 34	RNA-mediated gene reg.	Kwon et al. 2013
Nucleolin-like	Synthesis of ribosomes	Durut et al. 2015
Cytokinesis, endocytosis
α-Tubulin	Cytoskeleton	Cho and Cavalli 2012
Dynein heavy chain	Cytoskeleton	Rishal and Fainzilber 2014
Reps1 associated Eps domain	Endocytosis of GF receptors	Dergai et al. 2010
Rho GTPase activating protein 1 (srGAP1)	neuronal cell migration	McKerracher et al. 2012
Vps13/Chorein	Actin regulation	Ueno et al. 2001
ARF (ADP-ribosylation factor)	Membrane traffic	Donaldson and Jackson 2011
Cell signaling
Faf1 (Fas associated factor)	Suppression of β-Catenin	Zhang et al. 2011
Fibroblast Growth Factor (FGF) receptor 3-related	FGF-Gremlin inhibition	Verheyden and Sun 2008
Inositol-3-phosphate synthase	Inositol synthesis	Deranieh et al. 2013
Novel regeneration-specific Hydra genes/proteins
comp27033_c0_seq1	Hydra regeneration specific 1
comp27457_c0_seq1	Hydra regeneration specific 2
comp27949_c2_seq5	Hydra regeneration specific 3

aBLAST orthology.

In the regeneration transcriptome (axis 1) no specific enrichment of any category was detectable (table 1; supplementary tables S4 and S7, Supplementary Material online). This group comprises genes for several membrane bound metalloproteases, which can act both as receptors and enzymes (Nistico et al. 2012). We also found genes encoding discoidin domain receptor tyrosine kinases, which are activated upon binding to collagen (Vogel et al. 1997; Tran et al. 2004; table 1; supplementary table S2, Supplementary Material online). Further upregulated categories include signaling genes in morphogenesis and phagocytosis, for example, genes encoding endocytosis, JNK, Erk (MAPK), and mTOR pathways, Ran- and Rho-activating proteins, guanine nucleotide exchange factors, including DOCK proteins activating the Rho-family GTPases Rac and Cdc42 (Cote and Vuori 2007; Kulkarni et al. 2011) and genes encoding proteins involved in apoptosis, for example, the death receptor-associated protein Daxx (Chu et al. 1995) and Gsk3.

Early and Late Molecular Signatures

To further investigate the dynamics of genes including those that were missed by proteomics, we analyzed the transcriptome data set alone. This also allowed us to trace the dynamics up until 48 h after injury when the regeneration process has been finished. The response at 24 and 48 h constitutes a set of genes, which is partly new or a subset of the patterning regeneration response, indicating an end to the prepatterning stage by 12 h and a transition to the differentiation program.

Figure 3 represents PANTHER molecular function categories (at the highest hierarchy level) identified in the Hydra transcriptome. The same categories were found in the up- and downregulation data sets (fig. 3A). For example, several members of the PANTHER Signaling Molecules category (MF00016) including genes that encode TGF beta ligands (see below) and Wnt inhibitors, like the secreted frizzled-like proteins, are initially downregulated (fig. 3B), whereas members of the Wnt ligand secretion complex are significantly upregulated from 0.5 h onward (fig. 3D; see also below). Unexpected was the high regulation of nucleic acid binding genes (this category excludes transcription factors). While having several ribosomal genes that are downregulated at 0.5 h (e.g., L2, fig. 3C, correlating with the translational arrest), it also contains several transposable element derived loci, which become upregulated from 0.5 to 12 h. This expression returns to the preinjury levels only after 12 h (fig. 3C) and shows significant correlation with the activation of DNA methyltransferases and isomerases (fig. 3A).

F

Transcriptome dynamics. (A) Heat map of the number of genes significantly up (upper panel) and down (lower panel) regulated at a given time point (see text for details) in different PANTHER Molecular Function categories. Heat map colors: green: higher amount, magenta: lower amount. (B and C) Differential up- and downregulation within a single molecular function group is shown on panel B for Wnt pathway and on panel C for nucleic acid binding proteins (ribosomal protein L2 and transposons). (D) In situ hybridization (ISH) showing the expression pattern of differentially upregulated members of the Wnt pathway.

As earliest upregulated factors, we identified genes encoding the two AP-1 transcription factors Jun and Fos-related antigen (table 3). Both proteins form a heteromeric complex that is activated by the mitogen-activated protein kinase phosphorylation cascade (Eferl and Wagner 2003; Brenner 2014). Factors with a putative stem cell function are encoded in orthologous genes of the Ets family (Elk1/Ets2), which act in oncogenic transformation and early embryonic patterning in mice (Sharrocks 2001; Georgiades and Rossant 2006; Islas et al. 2012), but also in cell cycle progression (e.g., beside Elk1 there were further early-activated cell cycle factors, for example, Runx and the suppressor of hairless (Drosophila) homolog Rbpj (CBF1). Of interest was also a gene encoding the zinc finger protein FEZF (forebrain embryonic zinc finger), which functions as an intracellular factor attenuating Bone Morphogenetic Protein (BMP) signaling in sea urchins and vertebrates (Deheuninck and Luo 2009; Yaguchi et al. 2011).

Table 3.

Molecular Signatures in Hydra Regeneration.

Factor	Cell Type	Transcriptome Comp ID	Regeneration Dynamics	Amniota	Anamnia	Planaria
Regeneration factors
Jun	Ecto, endo	16707_c0_seq1	Early up, 0.5 h	+	+	+
Fos-like	Ecto, endo	21748_c0_seq1	Early up, 0.5 h	+	−	+
Brachyury	Ecto (?)	19584_c0_seq1	Early up, 3 h	+	+	+
Elk/Ets2	n.d.	25232_c0_seq1	Early up, 3 h	+	−	+
Fez	n.d.	22188_c0_seq1	Early up, 3 h	−	−	−
FoxA2	endo	23424_c0_seq2	Late up, 24 h	−	−	+
FoxN2/3	n.d.	22969_c0_seq1	Late up, 24 h	−	−	+
FoxQ2b3	i-cells	16080_c0_seq1	Late up, 24 h	−	−	+
Otp	n.d.	18117_c0_seq1	Late up, 24 h	−	−	+
Hydramacin	n.d.	28179_c0_seq1	Late up, 24 h	−	−	−
Ncoa5	i-cells	17486_c0_seq1	Late up, 24 h	−	−	+
Lhx1/5	endo	23270_c0_seq1	Late up, 24 h	−	−	+
FoxA/budhead	endo	41348_c0_seq1	Early down; late up,12 h	+	−	+
Gsc	endo	6976_c0_seq1	Early down, 3 h	+	+	−
OtxB	i-cells, endo	17873_c0_seq1	Early down, 3 h	−	−	+
CnVas1	i-cells, endo	18466_c0_seq2	Early down, 3 h	−	−	−
FoxQ2 /FoxF	n.d.	14694_c0_seq1	Early down, 3 h	−	−	−
Piwi	i-cells, ecto	09094_c0_seq	Early down, 3 h	−	+	+
Lmx1	n.d.	26323_c0_seq2	Late down, 12 h	−	+	−
Periculin	i-cells	21062_c0_seq1	Late down, 12 h	−	−	−
Dmbx (manacle)	Ecto	22739_c0_seq1	Late down, 12 h	−	−	−
FoxQ2b4	n.d.	20704_c0_seq1	Late down, 24 h	−	−	−
Unc4.1	n.d.	06824_c0_seq1	Late down, 24 h	−	−	−
Stem cell factorsa
Zfp46	Ecto	10542_c0_seq1	Early up, 0.5 h	−	−	−
Ets1	Ecto	31467_c0_seq1	Early up, 3 h	+	−	−
Sp4	endo	18844_c0_seq1	Early up, 3 h	−	−	+
Tcf	All	16698_c0_seq1	Early up, 6 h	+	+	+
FoxK1	i-cells	19812_c0_seq1	Late up, 24 h	+	+	−
Sox14	i-cells, ecto	23939_c0_seq1	Late up, 24 h	−	−	−
Srf	i-cells, ecto	19948_c0_seq1	Late up, 24 h	+	−	−
c-myc (myc1)	i-cells ecto	19544_c0_seq3	Late up, 24 h	+	+	−
Rfx3	i-cells, ecto	19747_c0_seq1	Late up, 24 h	−	−	−
FoxO	All	26218_c0_seq1	Late up, 24 h	+	+	−
Otx1/2	i-cells, endo	31182_c0_seq1	Early down, late up,12 h	−	−	−
Ski	All	25113_c0_seq1	Early down, 3 h	+	−	−
Zic2, Zic3	endo	22460_c0_seq1	Early down, 3 h	+	−	−
Klf3	i-cells	20529_c0_seq1	Early down, 6 h	+	−	−
Klf8, Klf11	Ecto	20529_c0_seq1	Early down, 6 h	+	−	−
Hmgxb3	i-cells	28252_c0_seq1	Late down, 12 h	−	−	−
Cux1	i-cells	25132_c0_seq1	Late down, 24 h	+	−	−

Dynamics of lineage-specific putative stem cell factors identified in (Hemmrich et al. 2012). To avoid redundancy, early upregulated genes encoding for Stem cell factors were not listed extra in Regeneration factors.

Among the early downregulated factors are genes encoding the vasa-related zinc finger transcription factor Cnvas1, which is involved in germ cell formation in Hydra and higher metazoans (Mochizuki et al. 2001), the Hydra Goosecoid ortholog Cngsc, which was described to get downregulated during Hydra head regeneration (Broun et al. 1999), FoxQ2, which (besides other functions) is linking the Wnt and Nodal signaling in sea urchins (Yaguchi et al. 2008), and an orthologous gene of Otx defining the anterior end in bilaterians (Smith et al. 1999).

Cell Type Specific Signatures in Hydra Regeneration

We next analyzed Hydra’s three stem cell populations for regeneration-specific signatures. Previous work in the H. vulgaris strain AEP has identified stem cell specific molecular signatures in each of the three lineages (Hemmrich et al. 2012). Table 3 shows that among 17 lineage-specific putative stem cell factors, 4 were strongly upregulated during early regeneration. The earliest gene encodes for the zinc finger factor Znf436, which is a negative transcriptional regulator of MAPK signaling (Li et al. 2006), followed by Ets1 (member of the E26 transformation-specific family of transcription factors that are associated with the malignant transformation of hematopoietic stem cells in humans [Sharrocks 2001; Dittmer 2003]), Sp4 (member of the Sp/Krüppel-like factor gene family [Suske et al. 2005]), and Tcf, encoding the crucial transcription factor in canonical Wnt signaling (Cadigan and Waterman 2012). The earliest downregulated gene was Ski, a proto-oncogene and potent negative regulator of TGF beta signaling (Luo 2004; Deheuninck and Luo 2009).

Cell type specific transcription factors that were downregulated between 6 and 48 h included several zinc finger transcription factors. We tested the cell type specificity of three related Zn Finger transcription factors by fluorescence-activated cell sorting (FACS) using cell type specific marker genes (i.e., transgenic Hydra strains with a nanos::GFP marker for the interstitial cell lineage, and strains in which either the ectoderm or endoderm was labeled with GFP or RFP, respectively). HyZic1 was expressed in the interstitial cell lineage, HyZic2 was specifically expressed in the endodermal cell lineage, and Krüppel-like Factor 3 predominantly in the ectodermal cell lineage (supplementary fig. S1, Supplementary Material online). HyZic1 was described to be an early switch in nematocyte differentiation (Lindgens et al. 2004). However, our data demonstrate a global downregulation of Krüppel-like Zn finger transcription factors in all three germ layers. This might indicate a function in proliferation and the survival of stem cells, similar to vertebrates (Nandan et al. 2014).

Of particular interest was the stem cell factor FoxO, which has been described as a main Hydra stemness factor (Boehm et al. 2012) and stress-response gene (Bridge et al. 2010). Surprisingly, FoxO became upregulated only after 24 h. Also Cnox-2, which was described to be involved in neurogenesis and to attenuate regeneration when blocked (Miljkovic-Licina et al. 2007; Brockes and Kumar 2008; Galliot et al. 2009), was not significantly regulated during regeneration.

Signaling Pathways Activated during Hydra Head Regeneration

We found a number of different pathways to be regulated during regeneration and most of them have not been characterized in Hydra regeneration so far (supplementary figs. S2–S4 and tables S5–S7, S9, and S10, Supplementary Material online). Among the earliest upregulated genes during regeneration were those involved in cell cycle and cell death regulation as well as pathways involved in MAP kinase signaling and proteolysis (supplementary fig. S2, Supplementary Material online).

Of particular interest were the dynamics of regulatory networks that have been reported to be active during regeneration in bilaterians (Poss 2010) and Hydra (Lengfeld et al. 2009). We found members of the Wnt pathway including Wnt ligands to be continuously upregulated over time (fig. 3B; supplementary fig. S3 and tables S14 and S16, Supplementary Material online). One of the first components that increased already at 0.5 h was HyWnt3 followed by HyWnt7, -9/10 a, -9/10 c, and -16 at later stages. This upregulation corresponds to a significant elevation in the expression levels of the downstream components (e.g., Frizzled5/8, Dvl3, β-catenin, and Tcf/Lef) and it continued past 12 h. We have validated these expression dynamics by quantitative polymerase chain reaction (qPCR; supplementary fig. S4 and table S11, Supplementary Material online). We also found essential members of the Wnt pathway that have not been identified so far in Hydra, for example, the LRP5/6 co-receptor, which is funneling a number of different Wnt ligands, is also strongly upregulated (supplementary fig. S3, Supplementary Material online). Of particular interest were several orthologous antagonists and agonists of Wnt signaling, that were early upregulated, for example, the extracellular Wnt target gene wingful (wf) (named Notum by fly base) encoding for a pectin acetyl esterase that acts as a Wnt feedback inhibitor, the cytoplasmic inhibitor Naked cuticle (Nkd), and the Wnt secretion factor Wntless (Wls/Evi). These factors are all involved in the long-range control of Wnt gradients and in Wnt-dependent patterning. In all cases we found a sharp and early upregulation (fig. 3D). For the main enzyme of canonical Wnt signaling, Gsk3, we observe a sharp downregulation in the transcriptome and proteome data sets after 3 h. β-catenin activation by Gsk3 inhibition is an essential step in the activation of Wnt signaling (Li et al. 2012). Gsk3 is also a negative regulator of proteins involved in metabolism, transcription, translation, cell cycle, anti-apoptosis, and signal transduction.

The dynamics of TGF beta/Smad signaling pathways were quite different (supplementary table S15 and fig. S6, Supplementary Material online). The expression of BMP2/4 and BMP5/8c was downregulated (supplementary fig. S5, Supplementary Material online), despite a transient increase in the amounts of BMP2/4 at 12 h. Also, the expression of genes encoding Activin ligands (supplementary figs. S5 and S6, Supplementary Material online) was downregulated at 24 and 48 h (supplementary figs. S2 and S3, Supplementary Material online). The expression of two main inhibitors of BMP ligands, that is, Noggin and Chordin (Rentzsch et al. 2007) was strongly upregulated. These data suggest that BMP and Activin signaling gets downregulated during Hydra head regeneration (see Discussion). However, on the level of the cytoplasmic transducers, the receptor Smads (Smad1/5/8 and Smad2/3) and Co-Smad (Smad4) were significantly upregulated during regeneration (supplementary fig. S5, Supplementary Material online). The reason for this upregulation is currently not understood, but we presume that it is related to the interaction of Smads, which serve as a node for signal integration, with several other signaling pathways, including Ras-MAP-kinase signaling, Erk (MAPK), CDKs, proteolysis, and apoptosis (Massague et al. 2005).

Phosphoproteome Reveals Activated Proteins

In our SILAC phosphoproteome, we analyzed four biological replicates each measured in two technical replicates, either on an LTQ-Orbitrap XL or QExactive mass spectrometer. In total, we could identify 6,629 phosphorylation sites with a false discovery rate (FDR) of 0.01. Four thousand six hundred fifteen sites could be quantified, whereas 780 sites were significantly regulated during 12 h regeneration (see Materials and Methods for further details; supplementary table S8, Supplementary Material online). As shown in figure 4, four distinct clusters, that is, phases of phosphorylation events could be identified, an immediate phosphorylation/dephosphorylation phase at 0.5 h and later phosphorylation/dephosphorylation events from 3 to 12 h of head regeneration. To identify potential kinases, which are responsible for these changes, we assigned kinase substrate motifs (supplementary fig. S7, Supplementary Material online) to the identified phosphorylation sites and tested for specific enrichment of certain motifs in each cluster.

F

Heat maps of significantly differentially regulated phosphorylation sites. Heat map colors: green: higher amount, magenta: lower amount. Significance was tested using multiple-sample ANOVA, and phosphorylation sites were accepted to be significantly regulated with a FDR of 0.01. The corresponding sites were clustered according to their dynamic regulation. Kinase substrate motifs were assigned to each phosphorylation sites and tested for significant enrichment in each cluster using Fisher exact test with an FDR of 0.04. Significantly enriched substrate motifs are shown next to the cluster, whereby motifs, which appear uniquely in one cluster, are highlighted in bold letters.

Interestingly, Casein kinase 1 and 2 motifs seem to be uniquely overrepresented in the immediate regeneration responses (fig. 4). Although the activity of Ck1 is reduced, Ck2 becomes activated within 0.5 h of Hydra head regeneration. Both kinases are well known to play a central role in the activation and control of the Wnt signaling pathway (see Discussion). Among the dephosphorylation motifs the family of cyclin-dependent kinases (CDKs) was specifically enriched at 0.5 h (fig. 4). Protein phosphorylation by CDKs controls the cell cycle by modifying key regulators for DNA replication and mitosis including Cyclins. The strong dephosphorylation of various CDKs during the onset of Hydra head regeneration could explain the cell cycle arrest observed at the site of regeneration (Holstein et al. 1991). In line with this is the dephosphorylation of Bard1/Brct domain binding motifs. Brct domain proteins are involved in cell cycle control and DNA damage including the name giving family-member Brca1/2 (Bork et al. 1997).

Among the phosphorylated proteins, we identified an immediate and a highly specific response (peak at 0.5 h) representing a large number of proteins that are known to be involved in the stress and injury response, in signaling and wound healing, and the regulation of autophagy and apoptosis (supplementary table S13, Supplementary Material online). Members of this first cluster were Uromodulin, which has a function in wound healing of sea anemones (Dubuc et al. 2014), histone deacetylase 4 (Hdac4), and various G protein modulators. Of interest was also the combined phosphorylation of Sgk3 (serum and glucocorticoid-regulated kinase 3) and PI3K (phosphatidylinositol-3-kinase). Phosphorylated Gsk3 is part of the stress induced PI3K/Act signaling pathway and can provide a link to β-catenin/Tcf signaling by inhibiting Gsk3 (McCormick et al. 2004).

The immediate dephosphorylation response (cluster 2) included various autophagy apoptosis inhibitors (e.g., Foxk), GTPases, and exchange factors including Rho and Ras GEFs, but also RNA processing proteins. In terms of patterning genes, we found that Trim33/TIF1γ, a factor promoting ectodermal and dorsal fate in vertebrates, was dephosphorylated. Trim33/TIF1γ is a monoubiquitin ligase for Smad4, the signal transducer of BMP signaling, and thereby acts as an antagonist for BMP gradient formation (Dupont et al. 2005, 2009; Wisotzkey et al. 2014). The dephosphorylation of Trim33/TIF1γ might be in line with the inhibition of BMP signaling observed in our SILAC analysis, although data on the function of phosphorylation for this factor are missing.

Cluster 3, which contains factors getting phosphorylated in the late stage of regeneration, is mainly characterized by two groups of proteins: factors involved in membrane trafficking and cytoskeletal reorganization, and factors regulating stemness. The regulation of vesicular trafficking by endocytosis-associated factors as Epsin and Auxillin and of general stress-responsive proteins involved in cytoskeletal organization like Adducin and p21-activated kinases could be associated with the ongoing apoptosis/autophagy response initiated at the early regeneration phase while phosphorylation of CDK1-regulated Filamin and of Talin, a factor in integrin signaling, indicate the onset of mitosis and epithelial reorganization. Of interest are several stem cell factors such as Buttonhead, which maintains neural stem cell identity in Drosophila by promoting the generation of neural progenitors (Jiang and Reichert 2014), Hdac1 and Rif1. Similar to the Wnt-induced activity of Sprouty, an Fibroblast Growth Factor (FGF) signaling inhibitor, which was also found in this cluster, their activity might be dedicated to the growth regulation of progenitor cells.

Cluster 4, which comprises late dephosphorylation events, is probably characterized by reciprocal events to cluster 3, because also here factors regulating membrane trafficking and cytoskeletal reorganization such as Nexin, Enthoprotin and proteins involved in rho-GTPAse regulation, are found. Notably, we also identified Piwi in this cluster, a recently characterized factor of Hydra stem/progenitor cells (Juliano et al. 2014).

Evolutionary Novel Genes and Genome Architecture

In the light of highly variable regeneration mechanisms among metazoans, we sought to characterize the genomic dynamics of genes involved in Hydra regeneration. We approached this by studying novel genes and by analyzing the genomic architecture of all genes involved in regeneration.

Novel genes have already been implicated in regeneration (Brockes and Kumar 2008). We therefore searched for novel genes in Hydra. Ideally, these are genes that have arisen when Hydra evolved within the hydrozoan lineage at least about 60 Ma (Martinez et al. 2010). However, the current genomic sampling of hydrozoans is insufficient to obtain a good set of novel genes for that phylogenetic node. We have therefore used the cnidarian node, and identified 382 transcripts that share no detectable hit with species outside cnidarians. Out of them 15 genes are either active during proteome or transcriptome regeneration response. This constitutes a specific enrichment of novel transcripts during regeneration (Fisher’s exact test P-value 0.003, or 2.7% vs. 1.1%). Interestingly, all were exclusively active during proteome regeneration response (MCIA axis 1, proteome; fig. 2C). Only one locus was also active during the transcriptome injury response. Interestingly, we found that the majority of novel genes (9/15) have neighboring, non-overlapping transcript mappings, mostly (7/9) at the N-terminal end (supplementary fig. S8, Supplementary Material online). The annotations of those additional peptides may provide clue about the function of the novelty. For example, we found an addition of an extra domain to hepatocyte growth factor-regulated tyrosine kinase substrate, Tropomodulin, and a cell adhesion protein afadin. The remaining six novelties constitute single novel genes, with multiple exons and complete ORFs. Additional genomic sampling using related Hydra genomes will help to reveal their evolution.

F

Genome architecture and transposable element impact on differential expression. Charts show genes upregulated at different time points (x axis). (A) Genes that are upregulated at 0.5 h have much smaller intron sizes (left), while still having similar exon sizes (right). (B) Transposable element counts (all classes) in the different gene structures, normalized by the total length (in bp). Higher repeat content is observed specifically in 5′-regions of genes active at 0.5 and 3 h and lower counts are observed for genes at 6 and 12 h. No information is available for intronic transposons due to the smaller intron sizes in early-regulated genes. (C) Transposable element age (Jukes–Cantor distance) distribution shows a shift toward younger repeat copies in the 0.5 h upregulated genes.

Using a previously described phylo-stratigraphy method (Domazet-Lošo and Tautz 2010), we further sought to characterize the distribution of novelties among differentially expressed genes (supplementary fig. S9, Supplementary Material online). We find that the earliest (wound response) events are mostly enriched in genes shared only among hydrozoans while in later stages (especially 24 and 48 h) genes with highest conservation are expressed.

We additionally investigated the genomic structure of genes involved in the regeneration response. By selecting genes with at least one intron and no exonic repeat content, we found that in early response genes the intron size was significantly decreased compared with late response genes while the flanking region contained younger repeat copies (fig. 5A and B). Previous reports (Castillo-Davis et al. 2002) indicate that short intron size genes are advantageous to transcribe due to the reduced transcriptional burden (see Discussion).

Discussion

We report a combined analysis of the molecular response reflected by the global transcriptome and proteome during Hydra head regeneration (see also summary fig. 6). Our data show a modest correlation between the two separate approaches, SILAC and RNAseq. These findings complement previous studies performed in yeast (Foss et al. 2007), plants (Fu et al. 2009), and mice (Ghazalpour et al. 2011). On the biological level, the observed differences can be explained by diverging transcriptional and translational efficiencies, a different turnover of transcripts and proteins, and a wide array of posttranslational modifications in proteins that are not reflected by the transcriptome. As demonstrated recently for the Wnt/STOP pathway, around 10% of all cellular proteins can be posttranscriptionally stabilized upon Wnt pathway activation in a β-catenin independent manner (Acebron et al. 2014). Also, Ghazalpour et al. (2011) report in their comparative study in mice a divergent extent of genetic regulation for transcripts and peptides. On the technical level, the divergence of the data sets is mainly due to a lower sensitivity of the proteomic analysis. This is reflected by the fact that several transcriptionally regulated genes are not identified in the proteome due to low abundance. In addition, genetic variation within the analyzed sample can principally contribute to a reduced depth of analysis in the proteomic approach.

F

Scheme showing the evolutionary and molecular dynamics of regeneration-specific genes and proteins (magenta) during Hydra head regeneration.

In our systematic investigation of Hydra head regeneration, we identified the most highly activated regeneration-specific genes by combining a global expression analysis (RNAseq) with an analysis on the proteome level uncovering quantitative changes in protein synthesis and degradation along with changes in the phosphorylation status. Our unbiased and statistically supported multiple coinertia analysis identified a fast response (at 0.5 h), correlating with an initial injury reaction, and a subsequent slow response that continues up to 12 h correlating with the prepatterning of the regenerating tissue. By integrating total and phosphorylated protein time series, we also identified multiple functional sites in early regenerates. Later stages (24 to 48 h) are correlated with the differentiation of the prepatterned tissue. The general picture emerging from these combined early proteome data indicates that a first step in the regeneration process encompasses anti-inflammatory responses and a broad arrest in cell cycle and cell signaling, which is followed by a global activation of metabolic and signaling process.

Our in-depth analysis identified an unexpected richness of factors ranging from responses in global scale metabolic networks (e.g., Gsk3, JNK, Erk/MAPK) to signaling (Wnt, TGF beta, Hedgehog, Notch), and transcriptional regulators (Elk, Ets1, Jun, Sp4), but also a number of novel factors. The increased production of membrane trafficking proteins (e.g., Syntaxin, Synaptotagmin, and Epsin) and cytoskeletal proteins (e.g., tubulin, afadin, and Villin) might be correlated to an extensive tissue remodeling at the side of injury. Our data also show that for factors shared among proteome and transcriptome (such as Gsk3 and mTOR) the main initial changes in the signaling pathways occur only at the proteome level. A starting point in the regeneration process appears to be the broad arrest in cell cycle and cell signaling, which is followed by a global activation of metabolic and signaling processes.

Comparing our results to previously published studies in cnidarians (Dubuc et al. 2014; Wenger et al. 2014), we find on the transcriptional level shared pathways that are acting primarily on the level of an injury and wound response and include factors required for apoptosis, clotting, cytoprotection, proteasome activity, and stress response. Because these factors were also found in wounded (punctured) N. vectensis tissue (Dubuc et al. 2014) they are likely injury-specific responses, similar to those that have been previously discovered in several bilaterian models (Li et al. 2007; Jiang et al. 2009; Lee et al. 2009). Our proteome study, however, additionally revealed factors involved in cell survival and stress response (table 2) that were not transcriptionally regulated. None of these early proteome response factors identified in our combined proteome-transcriptome data set has been discovered in previous regeneration studies on H. vulgaris (Wenger et al. 2014) or N. vectensis (Dubuc et al. 2014), and only a few are represented in the other data sets, that is, 2 in proteome regeneration, 3 in transcriptome injury, and 3 in transcriptome regeneration. This list includes factors involved in tissue remodeling, that is, matrix metalloproteinase (proteome regeneration, transcriptome regeneration), Leucine Rich Repeats (LRR) proteins (proteome regeneration), Calreticulin (transcriptome regeneration), Discoidin (transcriptome regeneration), and factors involved in stress response like DnaJ (transcriptome injury), Hsp90 (transcriptome injury), and Caspase (transcriptome injury) (Wenger et al. 2014). Moreover, among our early protein responses were factors that are well-known to induce a broad and general arrest of cell cycle and signaling, for example, Rif1, a major cell cycle arrest factor (table 3) or Faf1, a novel inhibitor of Wnt-signaling (Zhang et al. 2012). These factors not only explain the dramatic cell cycle dynamics in regenerating Hydra tissue (Holstein et al. 1991), they are also candidates for linking the injury response with the re-patterning of the regenerating tissue.

Currently, only a handful proteomic studies from other metazoans exist (Looso et al. 2012, 2013) mostly focusing on the role of stem cells in regeneration at a single time point. By comparing the SILAC data of proteins changed in the epithelial cell fraction (X1/Xin) from Schmidtea sp. (Boser et al. 2013) to their orthologs in Hydra, we found the highest correlation (Pearson correlation coefficient 0.25, F-test P-value 0.001) to the proteome at 0.5 h. Other time points showed smaller correlation (<0.2). This suggests a similar proteomic response in injury-induced networks in Hydra and stem cells in Schmidtea sp., including key regulators such as Gsk3, cell cycle proteins, and Dock. More coordinated cross-species studies are necessary to reveal common patterns in proteomic response.

Our analysis also reveals an understanding of the contribution and dynamics of the various pathways in Hydra regeneration, which was unclear so far. The transcriptome data demonstrate that an upregulation of Wnt signaling and an inhibition of BMP/Activin signaling can globally characterize both signaling cascades. Of particular significance are our phosphoproteome data demonstrating the early dephosphorylation of a main intracellular inhibitor of BMP signaling Trim33/TIF1γ, a monoubiquitin ligase for Smad4, which is transducing a BMP signal to the nucleus (Dupont et al. 2005, 2009; Wisotzkey et al. 2014). Despite the fact that it is not known so far how the activity of Trim33/TIF1γ can be regulated (phosphorylation vs. dephosphorylation), our data show that this is a sensitive checkpoint for the patterning of the newly regenerating head. We have recently shown that a Chordin-like inhibitor of BMP signaling is strongly upregulated in the early phase of regeneration suggesting an early inhibition of BMP signaling (Rentzsch et al. 2007). In line with these findings we could not identify any significantly increased phosphorylation or dephosphorylation of the receptor Smads. Only the co-Smad Smad4 was initially phosphorylated, which might be due to its interaction with other signaling pathways (Massague et al. 2005). Thus, the pivotal function of Wnt signaling for Hydra head regeneration (fig. 3; supplementary figs. S2 and S3, Supplementary Material online; Hobmayer et al. 2000; Lengfeld et al. 2009) is highly likely, although future work will show to what extent the regulatory switch between metabolism and patterning is governed by β-catenin/Gsk3 protein dynamics.

At this moment it is not clear how the injury response is translated into the Wnt controlled patterning process. In our initial work (Hobmayer et al. 2000), we found a 5- to 6-fold upregulation of β-catenin transcripts 0.5-1 h after head removal suggesting a direct linkage of injury and the activation of Wnt signaling. Our new RNAseq and qPCR experiments show, however, a gradual increase of β-catenin transcripts (supplementary fig. S4, Supplementary Material online). We presume that β-catenin levels are highly dynamic during the initial phase of regeneration suggesting further pathways operating. In liver regeneration, elevated levels of β-catenin proteins were found only 5 min after hepatectomy, followed by a decrease at 15 min onward (Monga et al. 2001). Accordingly, our proteome data set reveals a strong upregulation of a large set of proteins after 0.5 and 6–12 h, which is transcriptionally independent. And our phosphoproteome shows a rapid activation/inhibition of Casein kinases 1 and 2, both regulators of the Wnt pathway. Ck2 has a clearly defined positive role in the Wnt pathway (Song et al. 2003; Bryja et al. 2008; Bernatik et al. 2011; Cruciat 2014; de Groot et al. 2014), whereas different isoforms of Ck1 can have a complementary activity (Taelman et al. 2010; Acebron et al. 2014). We therefore assume that injury mediated activation/inhibition of these kinases initiates β-catenin transcription, which in turn activates β-catenin target genes including Wnt3 (Nakamura et al. 2011) that are all crucial for the later patterning process. In addition, previous studies using pharmacological inhibitors of MAP kinase have shown that several protein kinases are required for the development of the head organizer in Hydra (Fabila et al. 2002; Arvizu et al. 2006).

The fact that Wnt signaling can stabilize proteins in a global manner (Acebron et al. 2014) can explain the transcription independent protein upregulation after 0.5 and 6–12 h of Hydra head regeneration. The injury response leads to a ligand independent activation of the Wnt pathway, potentially mediated by Ck1/2 and/or Gsk3 which then activates β-catenin-dependent transcription of Wnt ligands and subsequently a ligand dependent activation of the pathway after 6–12 h. Thus, an altered state of kinase activities could lead to an initial and ligand independent activation of Wnt signaling. A putative candidate is the Wnt/STOP pathway, which has been recently shown to induce a global stabilization of proteins in response to Wnt activity. A β-catenin and ligand-dependent activation of Wnt signaling could be the essence of a second phase of Wnt activity at 6–12 h (see also Lengfeld et al. 2009). A mechanistic and biochemical proof of this hypothesis is beyond the scope of this study.

Our data can also shed light on an on-going debate on the cell-type specificity of the regeneration response in Hydra regeneration. Similar to the function of neurogenesis in newt (Kumar et al. 2007; Brockes and Kumar 2008), it was proposed for Hydra that de novo neurogenesis is required for normal head regeneration (Miljkovic-Licina et al. 2007; Brockes and Kumar 2008; Galliot et al. 2009). However, Cnox-2, which was described to be involved in neurogenesis and to attenuate regeneration when blocked (Miljkovic-Licina et al. 2007; Brockes and Kumar 2008; Galliot et al. 2009) was not significantly regulated during regeneration. Of particular interest was the stem cell factor FoxO, which has been described as a main Hydra stemness factor (Boehm et al. 2012) and stress-response gene (Bridge et al. 2010). Surprisingly, FoxO became upregulated only after 24 h. Our new data therefore strongly support previous findings demonstrating no major function of the nervous system or the interstitial cell lineage in Hydra regeneration (Marcum and Campbell 1978a,b; Sugiyama and Fujisawa 1978). This is underpinned by our finding that only protein clusters that denote late phosphorylation/dephosphorylation events are enriched by stem cell factors as Piwi.

Considering the initial onset of regeneration, we also made the unexpected discovery that genes specific to the injury response exhibit a small intron size and higher transposable element (TE) insertion counts, whereas genes that are highly active during the 3–12 h range are depleted in transposon insertions (fig. 5B). There are two scenarios to explain this result: either there is a strong selective pressure during the later regeneration program against TE insertions as they are likely to disrupt function of the core regeneration network, or the TE insertions have been used as a source of material for the evolution of new enhancers to drive early injury response. We presume that short intron sizes potentially reduce transcriptional burden of the early response genes and that the higher transposable activity in their flanking elements is a reflection of their evolutionary dynamicity. Hydra magnipapillata has one of the most repetitive metazoan genomes (Chapman et al. 2010) and its repeat content is estimated to be at around 50%, varying in closely related species. We expect that transposon insertions and their interference with evolutionary conserved stress response and patterning genes represent a major factor that can either foster or inhibit regeneration. Without further experimental evidence and functional testing we are unable to distinguish between the two scenarios. Interestingly, an upregulation of (retro-) transposons has been identified in salamander limb (Zhu et al. 2012) and radial organs in echinoderm (Mashanov et al. 2012a) regeneration. Their exact role, however, remains elusive (Mashanov et al. 2012b). Investigating transposon dynamics may thus help to explain why some closely related species significantly differ in their ability to regenerate (Scadding 1977; Brockes et al. 2001; Brockes and Kumar 2008). The enrichment of novel genes in the early injury phase underscores this novel aspect in regeneration indicating a taxon-specific adaptation process that might provide an explanation for the unusual regenerative potential of Hydra.

Materials and Methods

Hydra Regeneration Experiments

Experiments were carried out with H. magnipapillata (105 strain; Sugiyama and Fujisawa 1977) using 24 h starved, budless animals. Batches of ten animals were decapitated at 80% body length, transferred into fresh medium, and allowed to regenerate. Regenerative tissue (upper 10% of the body column) was isolated at 0, 0.5, 3, 6, 12, 24, and 48 h, immediately snap frozen (proteome) or dissolved in TRIzol (RNAseq). Three biological replicates were prepared for each time point, with 200–100 regenerates each. From these samples either proteins or RNA was extracted, and cDNA libraries were prepared (see supplementary materials and methods, Supplementary Material online).

Transcriptomics

The biological triplicates were sequenced for each time point, having 52989101/61382893/67032329 reads for time point 0 h, 31119141/27073313/21361307 reads for time point 0.5 h, 56640600/ 58074389/63795310 reads for time point 3 h, 51502370/36194810/86804837 reads for time point 6 h, and 40596586/57830400/23501305 reads for time point 12 h (single end reads). Time point 24 h had 28519851/28432965/25318135 reads and the 48 h time point had 31448683/31520325/30189391 reads. Reads were mapped against our RNAseq originated de novo assembly. For de novo assembly, reference mapping, differential expression analysis, qPCR, and in situ hybridization see supplementary materials and methods, Supplementary Material online.

SILAC Hydra and proteomics

For the proteomics approach, we used the auxotroph [13C6]lysine strain YAL6b of Saccharomyces cerevisiae, which was labeled with 1 mM [13C6]lysine (Gruhler et al. 2005). Because Trypsin is the prominent protease in shotgun proteomics, we first tried to label Hydra with the heavy version of Arginine (Arg10) and Lysine (Lys8), but observed an unnatural peptide isotope pattern. These deleterious effects were not observed in Lys6 labeled Hydra. Labeled yeast was than fed to freshwater cladocerans (Moina macrocopa), which in turn served as food source for the SILAC Hydra culture (fig. 1A). Twenty Hydra polyps were fed every second day with SILAC labeled and HUFA enriched M. macrocopa for 60 days up to a culture size of 2,000 animals (supplementary material and methods, Supplementary Material online). For the “Spike-in” standard preparation a mixture of regenerating and intact Hydra was used.

For each biological sample 200 regenerative tips were lysed, quantified and mix with the internal standard in a ratio of 1:1.25. Subsequently, the samples were separated by one dimensional (1D) Sodium dodecyl sulfate (SDS) Polyacrylamide gel electrophoresis (PAGE), cut in 27 fractions and digested by Lys-C. Before mass spectrometry analysis the samples were purified by C18 stage tips. See supplementary methods, Supplementary Material online, for further details.

For the phosphoproteomic approach, 300 µg unlabeled regeneration samples were mixed 1:1 with internal standard 2, digested with Lys-C and desalted with homemade solid phase extraction (SPE) cartridges using oligo R3 material. For phosphopeptide enrichment, we adopted the immobilized metal ion affinity chromatography (IMAC) based IMAC-IMAC method (Ye et al. 2010) with some modifications. Using this method we could separate multi- from single phosphorylated peptides, however only the single phosphorylated fraction was analyzed. Detailed information of SPE cleanup and phosphopeptide enrichment is described in supplementary materials and methods, Supplementary Material online.

Mass Spectrometry

Proteome samples were analyzed on an LTQ-Orbitrap XL mass spectrometer coupled to a nanoAcquity ultra performance LC system and phosphoproteome samples were analyzed on a LTO-Orbitrap XL as described above and on a QExactive mass spectrometer coupled to Agilent 1200 nanoflow LC system (Agilent Technologies; supplementary materials and methods, Supplementary Material online).

Data Analyses

LC-MS/MS raw data were analyzed with MaxQuant version 1.4.0.3 (Cox and Mann 2008) with default settings. Proteins were identified using an ORF-prediction of our de novo transcriptome assembly. Cysteine carbamidomethylation was used as a fixed modification, and methionine oxidation; protein N-terminal acetylation; and serine, threonine, and tyrosine phosphorylation (for the phosphoproteome data set) as variable modifications. For the identification, the FDR was set to 0.01 for peptides, proteins, and sites (the minimal peptide length allowed was six amino acids).

Protein ratios were normalized by its own Spike-in ratio as described in Looso et al. (2012). Using Perseus, protein ratios were additionally normalized for systematic errors by dividing with the median ratio. Protein IDs were filtered for proteins with at least three quantitative values in at least one time point.

The phosphoproteome data were analyzed using Perseus. Phosphorylation sites were filtered for a localization probability of 0.5 and additionally filtered for at least two quantitative values for each time point.

To define significantly differentially regulated phosphorylation sites multiple-sample analysis of variance (ANOVA) testing was performed using a cutoff FDR (fault discovery rate) of 0.01. Significant hits were isolated, normalized by Z-scoring and clustered using default settings. Linear phosphorylation motifs were assigned to each phosphorylation site, and Fisher’s exact testing was performed to identify significantly enriched phosphorylation motifs in each cluster. The sequence logo of each enriched linear motif was illustrated using WebLogo3 (Crooks et al. 2004).

To utilize the de novo transcriptome for high throughput SILAC-based proteomics approaches, we performed a multiple coinertia analysis (MCIA). We translated all transcripts in six potential ORFs and kept all ORFs longer than 100 amino acids as described earlier (Looso et al. 2013). Transcriptome data were additionally normalized and centered around zero with R’s (http://www.r-project.org/index.html, last accessed April 3, 2015) scale function to account for the different read counts in the different data sets. SILAC normalized protein counts were taken as is. MCIA was performed with the MCIA function from omicade4 package (http://www.bioconductor.org/packages/2.14/bioc/html/omicade4.html, last accessed April 3, 2015). Heat map was plotted with heatmap.2 function from gplots package (http://cran.r-project.org/web/packages/gplots/index.html, last accessed April 3, 2015).

In the Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis ≥5 members of a pathway must be present in the transcriptome to provide sufficient statistical support for the activation of a given pathway, at least three on the MCIA significance list.

Genome Structure and Transposons

Assembled transcripts were mapped to the Ringer-Phrap (RP) Hydra genome assembly (Chapman et al. 2010) using BLAT (Kent 2002) at 95% identity. Only one locus was allowed per transcript. The BLAT output was formatted to the General Feature Format (GFF) using blat2hints.pl tool from AUGUSTUS package (Stanke et al. 2006), and a custom perl script was used to compute the overlap of transcripts with repeats (based on a masked genome with RepeatMasker; Smit et al. 2007). Only genes with no exonic overlap with repeats and at least one intron with a length of 100 bp or longer were selected (to avoid counting recently inserted “pseudogenes”).

Supplementary Material

Supplementary tables S1–S16, figures S1–S8, and materials and methods are available at Molecular Biology and Evolution online (http://www.mbe.oxfordjournals.org/).