Conserved processes of metazoan whole-body regeneration identified in sea star larvae

1


INTRODUCTION 1 2
Questions of how regenerative processes evolve have stimulated researchers for centuries. Examples of 3 species with a capacity for restorative regeneration are distributed throughout the metazoan tree of life 4 ( Figure 1A), however the extent to which any animal is capable of regenerating is quite varied. It is only 5 vaguely understood which features of regeneration are common to species possessing robust regenerative 6 abilities and in what ways these features are limited in species with more restricted regenerative capacities. 7 While many attempts have been made to synthesize regenerative phenomena in disparate taxa (Bely and 8 Nyberg, 2010; Sánchez Alvarado and Tsonis, 2006), or to provide an evolutionary context to genes utilized 9 during regeneration within a model (Kao et al., 2013;Petersen et al., 2015), few studies have made direct 10 comparisons between highly regenerative organisms broadly across the metazoa. This is in part because we 11 are as yet missing detailed descriptions of regeneration from key taxa. A goal of making such comparisons is 12 to address to what degree conserved mechanisms underlie regenerative abilities, which has significant 13 implications for whether and how regeneration can be induced in organisms with more limited potential. 14 15 Some of the best known models for understanding the mechanisms of regeneration are species of Cnidaria 16 (e.g. Hydra vulgaris) and Planaria (e.g. Schmidtea mediterranea), which are able to regrow all body parts 17 following amputation, in what has been termed whole-body regeneration (WBR) (Bely and Nyberg, 2010). 18 From extensive studies of regeneration in these organisms, it seems that WBR in these contexts involves key distinction between these models lies in the source of new cells to be differentiated. While planaria 23 utilizes a pool of somatic stem cells (neoblasts) to generate a proliferative blastema that is essential for 24 successful regeneration, hydroids generally do not require cell proliferation for successful regeneration but 25 rather rely on transdifferentiation of existing cells to replace those lost by injury. However, hydra do possess a 26 pool of somatic stem cells (interstitial cells or i-cells) that are both undifferentiated precursors of various hydra 27 cell types (David, 2012) and found to proliferate following injury (Chera et al., 2009). One particularly important 28 question is whether regeneration recapitulates developmental gene regulatory networks (GRNs), and, to date, 29 the embryonic developmental of these animals is poorly described and hence this questions remains 30 unresolved. Species of Hydra belong to the Phylum Cnidaria, the radial animals, and planaria are bilateria 31 protostomes. A similarly well studied organism from within the other major grouping of animals (i.e. 32 Deuterostomes) has yet to emerge. This greatly limits our ability to identify mechanisms common to all 33 metazoan regeneration. 34  (Oulhen et al., 2016). However, to date no comprehensive survey of gene expression changes during 25 larval echinoderm regeneration have been reported. The potential impact of studying this process in larval 26 echinoderms is great given that this is one of the few deuterostome taxa with species capable of WBR along 27 with the availability of detailed developmental gene regulatory network information for these species 28 (Davidson et al., 2002;Hinman and Davidson, 2004). 29 30 Here, we provide a characterization of the molecular and cellular series of events that occur during larval 31 echinoderm regeneration and then assess the expression patterns of orthologous genes in other WBR 32 models. We employ regenerating larval asteroid Patiria miniata. We first characterize the regeneration time 33 course including the landmark events of increased cellular proliferation and reestablishment of gene 34 expression patterns associated with axis polarity specification. Then, sampling RNA from regenerating larval 35 segments and surveying gene expression changes by RNA-seq, we find broad classes of genes similarly 36 expressed in both anterior and posterior regenerating halves. Finally, through identification of orthologous 1 genes between P. miniata and published datasets of regenerating hydra and planarian models ( Figure 1A), we 2 find sets of genes that have similar temporal expression profiles in these different regenerating organisms. 3 These results highlight similarities in the regeneration programs of a bilaterian deuterostome, a 4 lophotrochozoan, as well as a basal eumetazoan. Such similarities suggest regeneration features common to 5 the base of the metazoans. To make an informed comparison to other regenerative models, we first characterized the broad stages of 12 larval regeneration in P. miniata. 7 day old bipinnaria larvae were bisected midway along the transverse 13 Anterior-Posterior (AP) axis ( Figure 1B). Both resulting larval segments were completely regenerative, 14 restoring all lost tissues and organs over the course of two weeks, similar to previous reports of larval asteroid 15 regeneration (Oulhen et al., 2016;Vickery and McClintock, 1998). Anterior segments were delayed in their 16 regeneration rate compared to their posterior counterparts. 17

18
Two aspects of the regenerative process were assessed for the presence and timing of potential landmark 19 events during regeneration. First, the pattern of cell proliferation was observed using 6 hour pulses of EdU to 20 mark actively cycling cells. The normal pattern of EdU incorporation is widely distributed, but notably 21 concentrated in the ciliary bands, mouth, stomach and coeloms ( Figure 2A). We infer from this result that the 22 larvae are actively growing in these regions. Upon dissection, the number of EdU + cells decreased by 3 days 23 post bisection (dpb), but were otherwise similarly distributed. Beginning at 6 dpb there is a marked change 24 with almost no EdU + cells localized in either the ciliary band or the gut and instead a dramatic localization of 25 EdU + cells proximal to the wound site in both posterior and anterior regenerating segments (Figure 2A). This 26 population of EdU + cells is localized at the regenerating edge, corresponding to the area where tissues form 27 during later stages of regeneration. We therefore infer that the larva has stopped normal growth and has 28 concentrated proliferation to the site of regeneration. Regeneration blastemas in other models are typically 29 described as mixtures of undifferentiated epidermal and mesenchymal cells that proliferate proximal to the 30 wound site in order to produce cells that will become differentiated and replace those lost by injury. 31 Therefore, we consider this region the regeneration blastema, although we have yet to define the cell types or 32 differentiation state of these proliferative cells. 33 Regenerating anterior segments (top row) and posterior segments (bottom row) demonstrate similar initial distributions of proliferation, although the number of EdU + cells decreased by 3 dpb. Beginning at 6 dpb, EdU + cells were concentrated near the wound site in both anterior and posterior regenerating segment in a putative regeneration blastema (bl). (B) WMISH staining shows that PmFoxQ2 is expressed in the apical domain (*) of intact larvae and regenerating anterior segments (top row). Posterior segments (bottom row) were negative for PmFoxQ2 staining early, but were stained at the anterior pole of regenerating posterior segments (*') by 5 dpb. (C) PmNk1 WMISH staining reveals staining in the post-oral ciliary band (**) in both intact larvae as well as regenerating posterior segments (bottom row). Anterior segments (top row) were not stained for PmNk1 early in the regeneration timecourse, but by 5 dpb there is staining on the posterior aspect of the anterior segment.
We next wanted to understand how genes with known roles in establishing the embryonic axis were 1 expressed during the restoration of these axes during regeneration. A great deal of work has established the 2 GRN for the specification of cell types along the P. miniata embryonic axis. PmFoxQ2 is expressed  uncut larvae as well as posterior regenerating segments but not in anterior segments during early 10 regeneration ( Figure 2C). Nk1 was first expressed in regenerating anterior segments at 5 dpb. This shows that 11 as early as 5 dpb both regenerating anterior and posterior segments have, at least in part, re-established 12 gene expression patterns corresponding to normal AP axis polarity. as those that are specific to the regeneration of each segment. Using RNA from non-bisected, age-matched, 35 sibling larvae as a control diminished the import of expression changes due to continuing larval development 36 or genetic differences between cultures. Genes significantly differentially expressed between each 1 regenerating segment and the control larvae were identified for each sampled timepoint. In total we identified 2 9,211 differentially expressed genes (DEG) across all assessed regeneration vs control contrasts. 3 4 We implemented a hierarchical clustering approach to distinguish segment-specific expression patterns from 5 those expression changes that are shared in both regenerating segments ( Figure 3A). In total, five separate 6 expression clusters were identified: (I) genes upregulated early in both anterior and posterior segments, (II) 7 genes downregulated early in both segments, (III) genes up in the anterior, down in posterior, (IV) genes up in 8 the posterior, down in anterior, and (V) genes up-regulated later (i.e. by 6 dpb) in both segments ( Figure 3A). 9 Thus, we have identified three subsets of DEG that have similar expression profiles across our time course 10 independent of which segment is measured (i.e. Clusters I, II, and V) and two subsets that are strongly 11 segment-specific (i.e. Clusters III and IV). To provide further insight into the shared functions of genes 12 assigned to each cluster, we examined the enrichment of Gene Ontology (GO) terms ( Figures 3B and S2). 13

14
Genes in Clusters I and II, those genes that are either up-or down-regulated early, respectively, in both 15 regenerating segments, are enriched for GO terms that indicate a robust wound response. Up-regulated 16 genes are enriched for terms that include cell signaling pathways (e.g. "MAPK cascade" and "calcium channel 17 activity"), "response to wounding", and "immune system process" (Figures 3B and S2). Furthermore the 18 enrichments for terms such as "neuron projection development" and "motile cilium" among up-regulated such as "ribosome biogenesis", "gene expression", as well as primary metabolism more broadly (e.g. 23 "mitochondrion" and "metabolic process"). These term enrichments are consistent with an early response to 24 the bisection insult that involves down-regulation of highly energetic cellular processes and up-regulation of 25 functions that are specific to the injury response (e.g. wound healing and immune signaling). 26

27
Clusters III and IV are composed of genes whose profiles are highly segment specific; these are genes 28 specifically down-regulated in one segment relative to the control larvae and unaffected or up-regulated in the 29 other. As we describe in more detail below, many of these genes are asymmetrically expressed along the AP 30 axis, therefore bisection results in the loss of posterior-specific gene expression from anterior segments and 31 vice versa. This is evidenced by the GO term enrichments for these clusters. Cluster III is 32 enriched for genes annotated with functions specific to anterior larval segments, such as "head development" 1 and "G-protein coupled receptor signaling pathway" (Burke et al., 2006), while Cluster IV is enriched for genes 2 associated with posterior fates, such as "peptidase activity" and "Wnt signalling pathway" (McCauley et al., 3 2013). Restoration of these gene expression patterns is a demonstration of functional regeneration and, as we 4 discuss later, we find that several genes within these clusters recover normal expression levels over this Although Cluster V is comprised of a relatively small number of genes, it is perhaps the most functionally 8 coherent cluster as the GO term enrichments are among the most statistically significant and most 9 reproducible across the three sources of functional annotations tested ( Figure 3B and Figure S2). Genes 10 assigned to Cluster V are enriched for terms related to the cell cycle, DNA replication, and ECM remodeling, 11 among others. Thus Cluster V genes, up-regulated by 6 dpb in both regenerating segments, are likely 12 associated with the onset of localized cellular proliferation observed during this same period of regeneration 13 ( Figure 2A). Importantly, these genes are up-regulated in both regenerating segments over control larvae that 14 also clearly contain many cycling cells. Thus, genes assigned to Cluster V may represent a regeneration-15 specific increase in expression of proliferation associated genes that is distinct from the normal, growth-16 associated proliferation. Finally, the transcriptome-based observations of AP axis polarity gene expression 17 restoration (i.e. Figure 4) prior to the up-regulation of cell cycle associated genes reinforces our whole-mount 18 staining results (i.e. Figure 2). 19 20 Several other trends are suggested by the analysis of GO term enrichments among these clusters. First, we 21 see enrichment of processes that are common to other models of regeneration such as "response to 22 wounding", tissue homeostasis (e.g. "cell adhesion", "extracellular matrix", "metalloproteinases"), cell 23 proliferation, and even "regeneration". Furthermore, there is enrichment of several specific cell signaling 24 pathways that are implicated in other regenerative contexts, including intracellular Ca 2+ signaling in wounding, 25 MAP kinase, WNT, and TGFβ/BMP signaling. Notably, several of these pathways are enriched in distinct 26 clusters, suggesting a differential employment of these pathways during specific phases of larval RNA was sampled only from the distal tip of regenerating aboral tissues during the 48 hours it takes to 11 achieve complete head regeneration. As blastemal proliferation is not a feature of hydra regeneration, this 12 characteristic cannot be used to consider synchronization of the regenerative phases in this study to the other 13 datasets. Nonetheless, these published datasets provide a reasonable basis for comparison to our sea star 14 dataset. 15 16 We identified orthologous genes between the sea star and each of the other datasets. In order to identify 17 which orthologs share similar temporal dynamics, the reported expression values from each dataset was 18 clustered. For each dataset we sought very coarse cluster classifications such that assigned genes were 19 either up-regulated early and down later or vice versa in their respective time course. The result is three 20 clusters each for the S. mediterranea and H. magnipapillata datasets ( Figure S3 and S4). Finally, we identified 21 genes in each sea star cluster with an ortholog in each of the defined planaria and hydra clusters. Using this 22 approach, we find statistically significant overlaps between genes differentially expressed early in all three 23 datasets as well as genes in the posterior-specific sea star cluster with clusters indicating segment specificity 24 in each of the other organisms. In the following sections we describe how this allowed us to identify not only 25 broad groupings of conserved expression patterns but also specific orthologs similarly expressed across 26 regeneration in these metazoans. 27 28 29 30 Expression of genes in Clusters III and IV, the anterior and posterior-specific clusters respectively, suggest 31 that these clusters are composed of genes asymmetrically expressed along the AP axis. Thus we impute that 32 early down-regulation of anterior genes in posterior segments is a result of the removal of anterior structures, 33

Axis respecification involves orthologous genes with distinct axial relationships in each model
and vice versa. The identities of many of the down-regulated genes support this as, for example, genes 34 normally expressed in anterior domains of the larvae (e.g. PmFoxQ2, PmFrizz5/8, and others) are initially 35 down-regulated in the posterior segments relative to uncut larvae but are unaffected in anterior segments 36 (solid lines, Figure 4; Cluster III, Figure 3). Correspondingly, genes normally expressed in the posterior domain 1 (e.g. PmWnt8, PmFrizz9/10, PmNk1, and others) are down-regulated in anterior segments but unaffected in 2 posterior segments (dashed lines, Figure 4; Cluster IV, Figure 3). We noted that many of these genes are 3 components of the previously described developmental GRN and, indeed, we detect a significant enrichment 4 of genes with orthologs in the GRN in the posterior-specific sea star cluster (Cluster IV, hypergeometric 5 p=9.97x10 -8 ). Restoration of normal larval gene expression levels along the AP axis is clearly a component of 6 regeneration (e.g. Figure 2B-C). Several genes in each of these clusters do recover normal larval expression 7 levels within the 6 days sampled, although the recovery of many posterior-specific genes is apparently 8 delayed in regenerating anterior segments compared to the recovery of anterior-specific genes in posterior 9 segments (Figure 4). Similar differences have been observed in the rate of planarian regeneration depending 10 on the position of the wound along the AP axis (Reddien and Sánchez Alvarado, 2004). 11 12 We examined the genes assigned to these segment-specific clusters to identify genes with similar expression 13 trends in the other datasets. In the planarian data, one cluster is comprised of genes asymmetrically up-14 regulated in the anterior relative to the posterior segments (i.e. Cluster 2, Figure S3), thus this cluster consists 15 of genes typically expressed in anterior regions of the worm. Paradoxically, a significant number of these 16 genes overlap the posterior-specific sea star cluster (Cluster IV, hypergeometric p = 1.4 x 10 -2 ) ( Figure S3). 17 The genes from the sea star posterior-specific set (Cluster IV) found to be overlapping this planarian anterior-18 specific cluster include PmFrizz9/10, PmWnt7, and PmSix1/2. Sea star Wnt7 corresponds to the planarian 19 Wnt2-1, and both SmWnt2-1 (Gurley et al., 2010;Petersen and Reddien, 2008) and SmSix1/2 (Lapan and 20 Reddien, 2011) are expressed in anterior domains in planaria. Given that we expect this planarian cluster 21 consists of anterior-specific genes, we also considered the genes that overlap the anterior-specific sea star 22 cluster (Cluster III). One of these overlapping genes is PmZic, which is also anteriorly expressed in these 23 worms (Vogg et al., 2014). 24

25
Since the hydra dataset reports on RNA sampled from oral regions of the regenerating aboral body stalk, we 26 reason that the signature of late stage up-regulation would reflect oral, or head, gene expression recovery (i.e. 27 Cluster 1, Figure S4). Here as well, there is a significant overlap of posterior-specific sea star genes (Cluster  Bode, 1999). HmWnt5 has not been associated with hydra axial patterning, but rather is associated with 35 tissue evagination involved in bud and tentacle formation (Philipp et al., 2009), and as such is expected to be 1 expressed in oral regions near the hypostome. The observed overlap in asymmetrically expressed genes between these datasets suggests a commonality of 4 the regulatory toolkit used for axis respecification, especially WNT ligands and receptors, during regeneration 5 in all of these models. Strikingly however, the orientation of the axes is not conserved, which likely reflects 6 developmental usage. 7 8 Expression of proliferation-associated genes occurs in distinct temporal phases in these models 9

10
The pattern and use of cellular proliferation is one aspect in which these models differ considerably. Both sea 11 star and planaria exhibit concerted wound-proximal, blastemal proliferation that coincides with the final time 12 points sampled here, 6 dpb for sea star larvae and 3 dpb for planaria. There is also a burst of global neoblast 13 proliferation observed early during planarian regeneration (i.e. within 6 hours) (Elliott and Sánchez Alvarado, are strongly associated with cellular proliferation and this is correlated with the emergence of blastemal 22 proliferation. We examined these genes for orthologs with similar expression dynamics in the other datasets. 23 None of the planaria or hydra clusters have a significant number of orthologous genes that overlap the genes 24 in this proliferation associated sea star cluster. The planaria and hydra clusters corresponding to up-25 regulation at later stages in each time course, i.e. planaria Cluster 1 ( Figure S3) and hydra Cluster 3 ( Figure  26 S4), consist of very few orthologs to the comparable sea star cluster. There is a stronger, though not 27 statistically significant, overlap between the genes up-regulated late in sea star and those up-regulated early 28 in planaria (e.g. Cluster 3, Figure S3) and hydra (e.g. Cluster 1, Figure S4), and many of the overlapping genes 29 in these comparisons are clearly associated with cycling cells (e.g. DNA polymerase subunits, MCM genes, 30 SMC genes, PmOrc3, PmRrm1, PmPlk, and PmTtk). Thus, gene expression associated with blastemal 31 proliferation in the sea star are more similar to the genes associated with observed proliferation at early 32 stages in both regenerating planaria and hydra models. 33 34 We also examined the expression patterns of genes implicated in blastema regulation in other models and 35 identified genes with expression patterns suggesting a related function in sea star larvae. Msx genes are 36 Figure 5. Conserved proliferation-associated genes. These data show sea star log fold change values for genes differentially expressed at later stages in regenerating segments compared with non-bisected sibling control larvae (i.e. sea star cluster V). All genes assigned to cluster V are plotted in grey. Several genes, either referenced in the text or representative of functions considered, are indicated with colored lines. The genes plotted with dashed lines are not part of cluster V, but may be associated with blastemal proliferation, as discussed in the text. Next to the key for each gene is an indication (i.e. "+") of whether an ortholog for that gene was found in an analogous cluster in either the planaria (S.m.) or hydra (H.m.) datasets. Indicators in brackets (e.g. "[+]") are those where no overlapping ortholog was identified by our analyses, but genes with the same name were implicated by published datasets. profile. Ultimately we find very little evidence for conserved temporal regulation of cell proliferation and 23 blastema formation, underscoring that this is an aspect in which these models may be quite divergent. This 24 discrepancy could reflect inherent differences in stem cell regulation between these models; for example, 25 unlike both hydra and planaria, there may not be a pool of pluripotent progenitors primed for proliferation 26 early following injury in larval sea stars. 27

28
The evolutionarily conserved early regeneration response 29 30 The strongest overlap between these WBR models is among genes differentially expressed early in each 31 dataset. We identified a significant number of orthologs that are up-regulated at early stages in both the sea 32 star and planaria, as well as the sea star and hydra datasets (hypergeometric p = 4.5 x 10 -3 and p = 8.8 x 10 -9 , 33 respectively) ( Figure S3 and S4). We also found a significant number of orthologs commonly down-regulated 34 at early stages in the sea star and planaria datasets (hypergeometric p = 3.3 x 10 -4 ). Genes down-regulated 35 early in both sea star and planaria are enriched for GO terms such as "ncRNA processing" and "ribosome", 36 indicating a common early repression of the energetically expensive process of ribosome biogenesis. 1 Meanwhile, genes that are up-regulated early across these datasets are enriched for GO terms including 2 "cilium", "calcium transport", and "signaling". Recent proteomic data indicate that calcium signaling is involved in the anterior regeneration in planaria 7 (Geng et al., 2015). MAPK signaling, through both ERK and JNK pathways, is important in the neoblast 8 control and blastema differentiation in planaria (Tasaki et al., 2011a(Tasaki et al., , 2011b. Studies in hydra have 9 demonstrated that wound-responsive MAPK signaling is necessary for early specification of the head 10 organizer, and thus functional regeneration, and that this may be a shared feature of highly regenerative 11 organisms (DuBuc et al., 2014). The observation of similar signatures in the genes up-regulated early in larval 12 sea star regeneration point to these signaling pathways being highly conserved aspects of the response to 13 wounding in these models of regeneration. 14 15 The enrichment of cilium associated genes (e.g. Ccdc11, Rsph3, Iqcd, and Iqub) ( Figure 6A) indicates that in 16 all three models a key set of genes up-regulated early are components of motile cilia. While this feature has 17 not been reported in either planaria or hydra, the role of cilia in wound response and regeneration has been 18 noted in mammals (Joiner et al., 2015) and zebrafish (Hellman et al., 2010), and a recent transcriptomic study 19 of regeneration in a related cnidarian (Nematostella vectensis) suggested a similar involvement of ciliogenesis 20 in that model (Schaffer et al., 2016). It is plausible that the up-regulation of cilia-associated genes reflects an 21 early need for sea star larvae to re-establish the ciliated bands, which are critical for effective swimming and 22 feeding. However, the common up-regulation in all of these models suggests that the generation of cilia may 23 be involved more generally in wound healing, perhaps in re-epithelialization of the wound or in communicating 24 cell signals. also up-regulated early in the other models, and down-regulation of ribosome biogenesis genes may be 34 associated with this response. There is also an early signature of general cell cycle arrest in the hydra 35 transcriptome (Petersen et al., 2015). While planarian neoblasts continue to proliferate at sites distal 36 Figure 6. Evolutionarily conserved early regeneration response. These plots show sea star gene log fold change values for genes differentially expressed early in both anterior and posterior regenerating segments compared with non-bisected sibling control larvae. Genes up-regulated in both segments (top row) correspond to cluster I and genes down-regulated in both segments (bottom row) correspond to cluster II. All genes assigned to each cluster are plotted in grey. Several genes, either referenced in the text or representative of functions considered, are indicated with colored lines. Next to the key for each gene is an indication (+) of whether an ortholog for that gene was found in an analogous cluster in either the planaria (S.m.) or hydra (H.m.) datasets. Indicators in brackets (e.g. "[+]") are those was no overlapping ortholog identified by our analyses, but genes with the same name were implicated by published datasets.
to the injury even during blastemal proliferation, inactivation of planarian PTEN gene homologs resulted in 1 defective regeneration due to neoblast hyperproliferation (Oviedo et al., 2008). These results indicate that a 2 common early feature of WBR in these systems is the modulation of regulators of cell proliferation. Genes associated with cell death pathways are also commonly differentially expressed early in these models. 5 We note conserved down-regulation of several genes in the autophagy pathway (e.g. PmAtg16L1, PmAtg12, 6 PmAtg10, PmAtg14, and PmUvrag) ( Figure 6B) and an increase in the expression of genes involved in 7 regulation of apoptotic and autophagic death (e.g. PmFadd, PmBirc6, and PmUlk1). An early repression of insult. It is plausible that these genes are a key shared characteristic between highly regenerative species is 32 8a specific response to injury that permits the regeneration program. However it is equally potential that the 33 key to understanding the commonality among these highly regenerative species is in how they interpret this 34 common early response and connect this to specific downstream regenerative processes. While the capacity for larval sea stars to undergo WBR has been appreciated for over two decades, there has 3 not yet been a systematic characterization of the cellular and molecular processes involved. In the present 4 study we demonstrate that larval sea stars exhibit many stereotypical characteristics found in other models of 5 WBR. Through our transcriptome analyses, we detect an early wound-response phase involving significant 6 alterations in the expression of stress response genes, signaling genes (including MAPK, Ca 2+ ) and a broad 7 shut-down of energetically expensive anabolic processes (e.g. ribosome biogenesis). The first few days 8 following bisection are marked by a global decrease in the number and distribution cycling cells compared to 9 what is typically observed in growing, non-bisected larvae. This coincides with, or perhaps slightly precedes, 10 the re-establishment of developmental axes, specifically the AP axis ablated by bisection. Re-patterning of 11 the AP axis is observed both through in situ hybridization as well as transcriptome measurements. These 12 observations are facilitated by our extensive prior knowledge of sea star developmental patterning programs 13 and, indeed, genes described by the developmental gene regulatory network are enriched in these clusters. 14 Notably, through both our transcriptome and in situ experiments we observe that axis respecification occurs 15 prior to the onset of wound-proximal cell proliferation, which is the last phase assayed in the present study. 16 This is the first description of concerted, wound-proximal cell proliferation in regenerating sea star larvae. 17 Given that this wound-facing region in both regenerating segments is the primordium from which larval 18 tissues will regenerate, we define this proliferative zone as the regeneration blastema. In this study we have 19 only monitored the first half of the regeneration process up until the emergence of this blastema. Complete 20 regeneration in these larvae takes a total of 10-14 days (Oulhen et al., 2016;Vickery and McClintock, 1998). 21 22 In this work we sought to leverage the power of comparing regeneration in a variety of contexts to identify 23 common features. For example, we clustered gene expression levels to identify genes similarly differentially 24 expressed in both anterior and posterior regenerating sea star larval segments. These patterns were then 25 used as a basis for comparison to planaria and hydra regeneration datasets. In the present study we 26 compared regeneration in species that last shared a common ancestor approximately 580 million years ago, 27 at the base of the metazoa. This is the broadest direct comparison of regeneration yet described, 28 encompassing three of the major groupings of animals (Deuterostome, Protostome, and basal Eumetazoan). 29 We find evidence for conservation of both broad functional classes as well as specific orthologs involved with 30 the regenerative process among these models. These commonalities are summarized in Figure 7. The most 31 remarkable signature of shared genes and processes is among genes both up-and down-regulated early in 32 these models. We are potentially most empowered to detect such an overlap among early genes as temporal 33 synchrony between the models likely diverges later in the time course. Nonetheless early changes to Ca 2+ and 34 Figure 7. Summary of similarities between WBR models. The reported features of regeneration at early, middle, and late stages of regeneration, with respect to the datasets considered in this study, are indicated. Features detected in the sea star model in our study that are shared with the other two models are highlighted in red. Some aspects are considered in common based on shared gene expression (e.g. MAPK signaling) whereas others are based on cytological observations (e.g. blastema proliferation).
MAPK signaling pathways, upregulation of ciliogenesis genes, upregulation of tumor suppressor genes, 1 downregulation of autophagy genes, and activation of a suite of immediate early genes are common aspects 2 of regeneration in these models. There is also a set of conserved genes that we hypothesize are commonly 3 involved in axial respecification, most notably genes in the WNT signaling pathway. Although the temporal 4 order of axial respecification prior to blastemal proliferation is common among these models, the genes 5 underlying the proliferative response are poorly conserved. 6 7 These commonalities between highly diverged WBR models highlights a deep conservation in regeneration 8 mechanisms among the metazoa. This work also underscores the potential power of such comparative 9 inquiries in identifying the core components of the regenerative response and, potentially, how these 10 components are altered in non-regenerative species.

RNA-seq, Read Mapping, and Transcriptome Assembly 4
For transcriptome measurements, larvae were grown and bisected as described in results. RNA was collected 5 from pools of approximately 300 sibling individuals of regenerating anterior segments, regenerating posterior 6 segments, as well as uncut control larvae. Two biological replicate samples were prepared for each timepoint. Echinobase were used to guide transcript assembly. Reads uniquely mapping to a gene (locus) from this 16 Cufflinks transcriptome assembly were counted (HTSeq-count v0.6.1p1, (Anders et al., 2015)). Read counts 17 were normalized and genes detected with more than 3 reads per million, corresponding to 50-120 uniquely 18 mapping reads depending on the sample, in at least two samples were retained for further analyses, 19 corresponding to 31,798 expressed genes. Raw and processed sequencing reads have been deposited into 20 the NCBI Gene Expression Omnibus (in progress) and analysis scripts are available upon request. 21

Gene Ontology term annotation and Ortholog identification 23
The newly assembled sea star genes were annotated in three ways: by identifying the reciprocal best BLAST 24 hit (rBBH) between the sea star transcript and either sea urchin or mouse genes and using Blast2GO. 9,027 25 (28.4%) loci have an rBBH match to a sea urchin protein, 7,212 (22.7%) loci have an rBBH match to a mouse 26 gene, and 9,617 (30.2%) assembled loci were annotated using Blast2GO. GO terms for each sea urchin and 27 mouse gene were assigned to their respective rBBH match in the sea star set and these were used for 28 enrichment analyses. Overall the results based on all three annotation methods are highly similar ( Figure 3B  29 and Figure S2). Reciprocal best BLAST hits (rBBH) were also used to identify putative orthologs between the 30 sea star genes and the planaria and hydra transcripts. We found 5,220 S. mediterranea transcripts and 6,091 31 H. magnipapillata transcripts with an rBBH match to a sea star transcript. 32

Differential Expression Testing and Hierarchical Clustering 34
Expression levels in biological replicate samples are highly correlated (pearson correlation coefficient = 35 0.985). Regenerating segments were compared to age-matched sibling uncut control larvae and differential 36 expression was assessed using a generalized linear model quasi-likelihood F test (edgeR, (Lund et al., 2012;1 Robinson et al., 2010)), controlling for sample batch. Differentially expressed genes (DEG) were defined as 2 those changes detected below a p-value of 0.05 and with a fold change greater than 2-fold in either direction. 3 Using these criteria there are 9,211 total DEG in at least one regenerating segment compared to the control 4 larvae and at least one of the timepoints sampled, which represents 28.97% of all of the expressed genes 5 detected. 6 7 The fold-change values for all 9,211 DEG relative to control larvae were clustered by first computing the 8 euclidian distance matrix and then these values were then clustered using the "ward.D2" method provided as 9 part of the R hclust function. The optimum number of clusters was determined by cutting the resultant 10 dendrogram at various heights and empirically determining at which height the number of clusters began to 11 plateau (h=42). The result was 8 distinct clusters. However, we noted that several clusters shared similar 12 overall patterns ( Figure S1). As the similar clusters shared very similar GO enrichments and expression 13 patterns over the time course, we further grouped these into the final 5 clusters reported in the text. The 14 grouping of clusters did not alter the enrichment of GO terms or our other downstream analyses ( Figure S1). 15 16 For the planaria and hydra regeneration datasets, data was obtained from supplemental tables associated 17 with each publication. The planarian data were reported as normalized read counts for the 15,422 transcripts 18 detected. These counts were log 2 -transformed and then scaled to z-scores, or the number of standard 19 deviations from the mean value for each transcript, and only those transcripts considered differentially 20 expressed as reported by the authors were considered. This resulted in 7,975 transcripts that were then 21 clustered in the same way as described above for the sea star transcripts. The hydra data were reported as 22 binned z-scores for the 28,138 transcripts detected corresponding to lower, mid, and upper third of 23 expression range for each transcript. We only clustered transcript values for which a positive reciprocal match 24 was detected, leaving 5,779 transcripts for our analyses. The euclidian distance matrix was calculated, as 25 with the other datasets, but to accommodate the binned nature of these data the hierarchical clustering was 26 performed using the "average" method provided with the hclust R function. A fine-grained resolution of 27 common gene expression dynamics across these species is not warranted without more closely aligning 28 experimental designs, including sampling time points and normalization strategies. Therefore, for each of 29 these datasets we sought very broad cluster classifications such that assigned genes are either either up-30 regulated early and down later or vice versa in their respective time course.The result is three clusters each 31 for the S. mediterranea and H. magnipapillata datasets (Figures S3 and S4). 32

33
Zheng for helpful feedback during the preparation of the manuscript. This work was funded by Charles 1 Kaufman Award from the Pittsburgh Foundation.