The prog-1⁺ (early progeny) and AGAT-1⁺ (late progeny) postmitotic cell populations constituting the first neoblast lineage described in planarians (Eisenhoffer et al., 2008) have been widely used as an assay for neoblast differentiation (Fraguas et al., 2011; Pearson and Sánchez Alvarado, 2010; Scimone et al., 2010; Wagner et al., 2012). The zeta-class neoblasts are a major subclass of planarian neoblasts identified molecularly by a specific gene signature, mainly two TFs, zfp-1 and soxP-3 (van Wolfswinkel et al., 2014). Zeta neoblasts are generated from the collectively pluripotent sigma-class neoblasts, and have recently been shown to generate prog-1⁺ and AGAT-1⁺ cells as well as other populations spanning the epidermis (van Wolfswinkel et al., 2014), suggesting that they are part of an epidermal lineage (Figure 1A). This lineage is a useful model to dissect the choreography of neoblast progeny differentiation because it is an abundant cell population undergoing rapid turnover, and both prog-1⁺ and AGAT-1⁺ cells are molecularly and spatially distinct. Prog-1⁺ (early progeny) cells likely represent a very transient cell population because they are lost about two days after a lethal dose of irradiation, whereas AGAT-1⁺ (late progeny) cells are not completely lost until about seven days post-irradiation (Eisenhoffer et al., 2008). Therefore, the precise molecular relationships between these cell populations and the mechanisms that control the progression and maturation of the zeta-class epidermal lineage still must be resolved.

Figure 1 with 1 supplement see all

Download asset Open asset

To further characterize the molecular transitions implicated in epidermal cell maturation, we compared whole worm expression profiles from two different RNAi knockdown conditions known to drastically reduce the population of AGAT-1⁺ cells. Previous studies have shown that the chromatin-remodeling factor chd4 and the tumor suppressor gene p53, both expressed in the neoblast compartment, are required for tissue homeostasis and regeneration (Pearson and Sánchez Alvarado, 2010; Scimone et al., 2010). Both chd4 and p53 are also expressed in AGAT-1⁺ cells. Chd4 is expressed broadly throughout the animal parenchyma as well as in the brain and ventral nerve cords, and p53 is expressed in discrete mesenchymal cells that also include smedwi-1⁺ and prog-1⁺ cells (Figure 1B). RNAi knockdown of both chd4 and p53 results in a dramatic loss of AGAT-1⁺ cells (Figure 1C). However, it is not well understood whether these mechanisms of action are direct or indirect. Based on their close proximity to the basement membrane and on post-irradiation kinetics, AGAT-1⁺ cells likely represent a relatively stable transition at which multiple fate decisions could be made. Therefore, we reasoned that the union of chd4 and p53 RNAi whole-worm RNA-seq down-regulated datasets would contain genes enriched in AGAT-1⁺ cells, as well as genes with enriched expression in epidermal cell populations that arise from AGAT-1⁺ progeny cells, potentially revealing critical factors required for post-mitotic differentiation.

We compared whole worm gene expression profiles of chd4 and p53 RNAi animals to those of control animals through multiple time points of RNAi treatment (Figure 1D). To verify the specificity and sensitivity of the data sets, we examined a number of control genes (smedwi-1, prog-1, AGAT-1 and chd4, p53) with known kinetics of disappearance after RNAi knockdown and found strong correlation (Figure 1—figure supplement 1A,B), indicating that known marker genes display predicted behaviors in our RNA-seq data.

To identify common down-regulated genes in chd4 and p53 RNAi-treated animals, we adopted criteria where such candidate genes would be required to have significant adjusted p-values (p.adj <1e-13) weighted towards multiple later time points after RNAi treatment (for chd4: d12, d15; for p53: d12, d15, d18). Using these parameters, a total of 1,250 genes were designated to be down-regulated after chd4 RNAi treatment, 1,174 genes were down-regulated after p53 RNAi treatment, and a total of 587 common genes were found to be significantly down-regulated in both data sets (Figure 1D—figure supplement 1C,D). Of the 587 common down-regulated genes, 70% (411/587) are predicted to encode homologs of proteins found in other organisms (Supplementary file 1,2 for up-regulated genes). Earlier studies have reported genes involved in creatine metabolism (AGAT-1, AGAT-2, AGAT-3), polyamine biosynthesis (odc-1), and monooxygenase activities (cyp1a1), as well as others with unknown function, to be expressed in similar cell populations (Eisenhoffer et al., 2008; van Wolfswinkel et al., 2014; Zhu et al., 2015). We analyzed the common gene set by assigning gene ontology (GO) terms and found that small molecule metabolic processes and transporter activity were the most overrepresented biological processes and molecular functions (Supplementary file 1).

We selected approximately 150 genes from the common down-regulated data set to screen by whole-mount in situ hybridization (WISH) (Supplementary file 1). We also included some genes that were uniquely down-regulated in the p53 RNAi data set (d18) because chd4(RNAi) animals that were to be collected for expression analysis died prematurely (before d18). A wide net was cast for candidate gene cloning, including genes with homologs predicted to be involved in a diverse array of biological processes and those with no known homologs, together spanning a wide range of expression levels. We performed an in situ expression screen to further narrow down our common gene list to those that are enriched in an AGAT-1-like mesenchymal or a similar epidermal pattern.

The majority of genes tested displayed distinctive WISH expression patterns (Supplementary file 3). We subsequently narrowed our candidate gene list to 29 unique genes based on a combination of their representative expression pattern, relatively strong signal intensity, and predicted gene function, for further analysis (Figure 1E). Based on colorimetric WISH, we classified the gene expression patterns of the 29 candidate genes into four main categories: discrete, AGAT-1-like sub-epidermal mesenchymal expression throughout the animal (AGAT-1 mesenchymal); mesenchymal cells that appeared more dense than AGAT-1 and/or discrete expression in the epidermis (AGAT-1 mesenchymal and epidermal); gut-enriched; and expression spanning multiple tissues exhibiting discrete mesenchymal/epidermal cells and in the gut/pharynx (Figure 2A). A wide assortment of biological processes were represented by the 29 candidate genes, including creatine metabolism (AGAT-1, gatm), monooxygenase activity (cyp3140A1, cyp3142A1, cyp3143A1, cyp3G1, cyp3141A1), fatty-acid metabolism (pla2, acsl-1, acsl-2), other metabolic processes (odc-1, mpv17, adss), protein-binding (traf-4), transporter activity (ttpal, slc25a-19), cytoskeletal (vim-3), DNA-binding/zinc-finger (egr-5, litaf, nkx-2.2, zfp-2, zfp-3), and multiple novel genes with unknown function. These novel genes with no known homology, predicted to encode small proteins (∼100 amino acids), contain a conserved signal peptide sequence. Because their expression patterns are down-regulated in zfp-1 RNAi animals (see below), we have named them zeta-class protein of unknown function (zpuf), followed by a unique designation number.

Figure 2 with 3 supplements see all

Download asset Open asset

We were surprised that many candidate genes from our screen were expressed prominently in the gut, as neither chd4 nor p53 have been previously reported to affect intestinal gene expression. For example, nkx-2.2 is a homeodomain TF that is enriched in the gut and has previously been shown to be required for proper intestinal regeneration (Forsthoefel et al., 2012). To verify the RNA-seq data by WISH, we selected a representative gene from each of the major expression categories and monitored their expression in chd4, p53 and zfp-1 RNAi backgrounds (Figure 2—figure supplement 1). The gene expression patterns for all five genes tested (zpuf-3, zpuf-6, vim-3, nkx-2.2, acsl-2) were significantly reduced in chd4 and p53 RNAi animals, whereas discrete mesenchyme and epidermal expression, but not gut expression, were reduced in zfp-1 RNAi animals (Figure 2—figure supplement 1A). We conclude that transcript levels in the RNA-seq data set are highly predictive of gene expression in vivo.

We next performed whole-mount fluorescent in situ hybridization (FISH) to determine whether the candidate genes exhibiting discrete mesenchymal cell expression patterns overlapped with AGAT-1 expression (Figure 2B). The majority of genes displayed substantial overlap with AGAT-1, but with varying degrees of signal intensity (Supplementary file 3). Several genes (zpuf-6, cyp3140A1, ttpal and ascl-2), despite showing substantial overlap with AGAT-1 ( >98%), were also expressed in additional mesenchymal cells (Figure 2B, yellow arrowheads). vim-3 was expressed in fewer sub-epidermal cells and showed very little overlap with AGAT-1.

After confirming that candidate genes with discrete mesenchymal cell patterns exhibit overlapping expression with AGAT-1 by whole-mount FISH, we next focused our attention on genes that are expressed in additional AGAT-1-negative (AGAT-1^neg) mesenchymal cells to determine any potential molecular relationships. We first characterized zpuf-6 (also known as NB.36.10A) because it has a robust expression pattern, is expressed in AGAT-1⁺ mesenchymal cells, and is also expressed in discrete epidermal cells (Figure 2C). We determined that the zpuf-6⁺ AGAT-1^neg mesenchymal cells do not overlap with collagen or prog-1 (Figure 2—figure supplement 2), suggesting that they are neither muscle cells nor early progeny cells. Rather, these zpuf-6⁺ AGAT-1^neg mesenchymal cells likely mark a unique population linked to epidermal cells.

We conducted triple whole-mount FISH for all other candidate genes with zpuf-6 and AGAT-1 to resolve whether they also overlap with zpuf-6⁺ AGAT-1^neg mesenchymal cells or represent some other heterogeneous cell type. Interestingly, all genes tested displayed significant or partial overlap with zpuf-6 in the mesenchyme and epidermis (Supplementary file 3). For example, pla2 and odc-1 both showed significant overlap with zpuf-6⁺ AGAT-1^neg mesenchymal cells and are very weakly expressed in zpuf-6⁺ epidermal cells (Figure 2—figure supplement 3). Furthermore, mpv17 expression overlaps with zpuf-6⁺ AGAT-1^neg mesenchymal cells but is not expressed in the epidermis (Supplementary file 3, data not shown). Acsl-2, adss and vim-3 were the only genes found to be expressed in additional epidermal cell types that did not overlap with zpuf-6. Adss, which encodes the enzyme adenylosuccinate synthase and is involved in purine biosynthesis, was the one exception that appeared to be expressed in additional AGAT-1^neg zpuf-6^neg mesenchymal cells. Taken together, our in situ expression analysis revealed that the candidate genes identified by RNA-seq are all expressed in similar, overlapping mesenchymal and epidermal cells marked by AGAT-1 and zpuf-6, and that these cell populations may all be constituents of the same lineage.

The spatiotemporal disappearance of various morphological markers in vivo can be used to identify cell populations that belong to the same lineage. Irradiation is a powerful tool that has been demonstrated to rapidly and specifically eliminate neoblasts, and subsequently, their immediate progeny cells (Eisenhoffer et al., 2008). The prog-1⁺ and AGAT-1⁺ cell populations were identified based on their spatial expression domains and down-regulation kinetics after irradiation. That is, the temporal order of prog-1 (early) and AGAT-1 (late) down-regulation correlates with the spatial distribution of these markers: the more peripheral the location of the cell, the longer the marker persists after irradiation. We were intrigued by a potential relationship between zpuf-6 and vim-3, also identified in our screen (Figure 2), because it is expressed in fewer sub-epidermal mesenchymal cells that do not overlap with AGAT-1, but more epidermal cells than zpuf-6. Therefore, both zpuf-6 and vim-3 from our candidate gene set expand the spatial distribution of known progeny markers into the epidermis (Figure 3A). The distinct expression domains of prog-1, AGAT-1, zpuf-6, and vim-3 suggest that these progeny cells undergo outward migration during the course of epidermal differentiation.

Figure 3 with 2 supplements see all

Download asset Open asset

We therefore extended this paradigm, that more differentiated cells have slower turnover kinetics, to determine the spatiotemporal irradiation kinetics of these newly characterized markers. After irradiation, both zpuf-6 and vim-3 displayed slower down-regulation kinetics compared to AGAT-1 and exhibited expression loss in a similar ventral-to-dorsal and anterior-to-posterior manner (Figure 3B, Figure 3—figure supplement 1). We also monitored the expression of rootletin, a gene expressed in the ciliated epidermis and tubule cells of the protonephridia (Glazer et al., 2010; Scimone et al., 2011), as a marker for differentiated cells and confirmed that its expression is not significantly affected after irradiation (Figure 3B, bottom panel). Therefore, the irradiation kinetics data support our interpretation that zpuf-6 and vim-3 mark post-mitotic cells further downstream of the same prog-1 and AGAT-1 epidermal lineage. In other words, AGAT-1⁺ cells likely give rise to zpuf-6⁺ and vim-3⁺ expressing cells.

We also looked for a conservation in AGAT-1 and zpuf-6 spatiotemporal expression kinetics in additional contexts. In both chd4 and p53 RNAi backgrounds, AGAT-1⁺ and zpuf-6⁺ cells are progressively lost in a ventral-to-dorsal, anterior-to-posterior fashion, with zpuf-6⁺ cells also undergoing a more delayed loss compared to AGAT-1 (Figure 3—figure supplement 2A). Interestingly, in p53(RNAi) animals, there is a characteristic banding pattern of AGAT-1⁺ and zpuf-6⁺ cells on the ventral surface of the animals (Figure 3—figure supplement 2A bottom panels, 2B) as cells are progressively lost over the course of RNAi treatment. This banding pattern of cells is also present for vim-3 on the ventral epidermis in p53(RNAi) animals at later time points (not shown), further suggesting that AGAT-1, zpuf-6, and vim-3 are markers exhibiting strong spatial correlative patterns.

Lineage tracing experiments using the thymidine analog bromodeoxyuridine (BrdU), in combination with whole-mount FISH, have revealed spatial and temporal regulation of neoblasts and the distribution of their differentiating progeny cells (Eisenhoffer et al., 2008; Newmark and Sánchez Alvarado, 2000). A single-pulse of BrdU has also revealed the turnover dynamics of prog-1⁺ and AGAT-1⁺ and epidermal cells, revealing that prog-1 and AGAT-1 become incorporated markedly earlier than epidermal cells, consistent with the notion that they mark an intermediate stage of epidermal differentiation (van Wolfswinkel et al., 2014).

We performed BrdU pulse-chase analysis to examine the temporal kinetics of BrdU-labeling from neoblasts to prog-1⁺, AGAT-1⁺ and zpuf-6⁺ cells during normal tissue turnover. The percentage of prog-1⁺ BrdU⁺ cells reached maximum levels around 10 days post-BrdU, followed by AGAT-1⁺ BrdU⁺ cells peaking at 14 days post-BrdU, and zpuf-6⁺ BrdU⁺ cells reached maximum levels at 22 days post-BrdU (Figure 4A). Confidence tests of the difference between individual data points show that through days 1-–14, the difference between AGAT-1⁺ BrdU⁺ and zpuf-6⁺ BrdU⁺ are only marginally significantly different (most significant data point p = 0.04 for d6). However, at 22 days post-BrdU the difference between the points of these same curves is highly significant (p <0.01). These results lend further support to our hypothesis that prog-1, AGAT-1, and zpuf-6 are distinct markers representing progressive stages of epidermal progeny differentiation.

Because zpuf-6 labels the first transition state into the planarian epidermis, we wondered whether zpuf-6⁺ epidermal cells embody a specific differentiated epidermal cell type or whether they still possess the potential to differentiate further. To address this question, we built upon the spatial and molecular relationships between AGAT-1, zpuf-6 and vim-3 defined by our previous irradiation kinetic studies by performing combinatorial whole-mount FISH. Triple whole-mount FISH of AGAT-1, zpuf-6 and vim-3 revealed that vim-3 overlaps with AGAT-1^negzpuf-6⁺ mesenchymal cells, zpuf-6⁺ epidermal cells, and there are additional zpuf-6^negvim-3⁺ epidermal cells present on the dorsal side of animals (Figure 4B). vim-3, predicted to encode an intermediate-filament like protein, also co-localizes with vim-1 (Figure 4—figure supplement 1A), which was shown to be down-regulated in zfp-1 RNAi animals (van Wolfswinkel et al., 2014). Intermediate filaments comprise a diverse class of molecules and are expressed in a variety of cell types, including epithelial cells. These filaments generally provide a scaffold to integrate components of the cytoskeleton and organize the internal cell structure (Snider and Omary, 2014). Thus, we wondered whether the expression overlap of zpuf-6 and vim-3 could reflect cells undergoing a morphological transition within the epidermis.

Figure 4 with 3 supplements see all

Download asset Open asset

The formation of cilia is a signature of terminal differentiation (May-Simera and Kelley, 2012). The ventral epidermis is lined with multi-ciliated cells responsible for gliding motility, while the more densely packed dorsal surface contains many non-ciliated and mucous-secreting, rhabdite-containing cells (Pedersen, 1976). Given that zpuf-6 expression is evenly distributed throughout the dorsal and ventral epidermis, we queried whether zpuf-6⁺ epidermal cells co-expressed markers for cilia genes. We used the planarian rootletin gene as a marker for ciliated epidermal cells and found that zpuf-6 and rootletin do not significantly overlap in expression (Figure 4C, top row), suggesting that zpuf-6⁺ epidermal cells are non-ciliated. However, the combination of vim-3 and rootletin in whole-mount FISH reveals that while cells expressing high levels of vim-3 display virtually undetectable levels of rootletin, there are cells with detectable vim-3 expression that do co-express rootletin (Figure 4C, bottom row), suggesting that vim-3 cells may begin to undergo ciliogenesis.

Under normal homeostasis conditions, zpuf-6⁺ epidermal cells appear randomly but broadly distributed on both the dorsal and ventral epidermal surface of intact animals (Figure 2). We reasoned that if zpuf-6⁺ epidermal cells still harbor the potential to differentiate further, then the density and distribution of zpuf-6⁺ cells would increase if epidermal integrity were severely perturbed. Therefore, we amputated wild-type animals and examined zpuf-6 expression in the regenerating blastema over 7 days (Figure 4—figure supplement 2). Consistent with our hypothesis, the density of zpuf-6⁺ epidermal cells as measured by fractional area is higher in the new undifferentiated blastema tissue compared to old tissue at 4 and 7 days post-amputation.

We also used whole-mount FISH to examine zpuf-6 expression with NB.22.1E and laminB, markers that are expressed in specific domains in the epidermis (van Wolfswinkel et al., 2014). NB.22.1E labels marginal adhesive gland cells at the body margin of the animal (Reddien et al., 2007; Tazaki et al., 2002) as well as cells lining the ventral mouth opening. Curiously, there was very little overlap between zpuf-6 and NB.22.1E at the edge, but NB.22.1E⁺ cells around the ventral opening co-expressed zpuf-6 (Figure 4—figure supplement 3A). LaminB is also expressed mainly in marginal adhesive cells and a few subepidermal cells near the edge, and displayed little overlap with zpuf-6 (Figure 4—figure supplement 3A). We wondered whether there was minimal co-expression between zpuf-6 and NB.22.1E and laminB because these marginal adhesive cells do not turnover as quickly as epidermal cells on the dorsal and ventral surface. To test this, we monitored the kinetics of NB.22.1E and laminB expression after a lethal dose of irradiation. After 14 days post-irradiation, NB.22.1E and laminB mesenchyme expression is completely lost, but only cells in the anterior region have disappeared, indicating that the marginal edge cells do not turnover as quickly (Figure 4—figure supplement 3B). Moreover, vim-3 overlaps considerably with both NB.22.1E⁺ and laminB⁺ cells at the animal body margin (Figure 4—figure supplement 1B), suggesting that both zpuf-6 and vim-3 are expressed in multiple epidermal cell types.

Taken altogether, transcriptional profiling of chd4 and p53 RNAi animals has identified additional epidermal progeny gene markers that are expressed broadly throughout the animal in discrete but overlapping cell populations in the mesenchyme and epidermis. zpuf-6⁺ epidermal cells can still undergo progressive differentiation into an unknown number of epidermal cell types (potentially including NB.22.1E⁺ cells lining the mouth and laminB marginal adhesive cells), by expressing markers of cytoskeletal morphogenesis (vim-1/vim-3) followed by cilia gene markers (rootletin). The addition of these new molecular markers in the planarian epidermal lineage will greatly facilitate the study of epidermal progenitor dynamics and differentiation during tissue homeostasis and regeneration.

Transcriptional changes are likely drivers of maturation within the epidermal lineage, which occurs in the absence of cell division. We therefore focused on characterizing putative transcription factors from our common down-regulated gene set as likely candidates important for post-mitotic epidermal differentiation. egr-5 is a planarian homolog of the early growth response (EGR1) family of C2H2-type zinc-finger TFs (Figure 5—figure supplement 1). EGR1 genes are known to be induced by extracellular signals including growth factors, hormones and neurotransmitters, and couple these signals to long-term responses by altering the expression of target genes (Thiel and Cibelli, 2002). Egr-5 is expressed in discrete mesenchymal cells located throughout the animal in a spatial pattern similar to AGAT-1 cells (Figure 2B). We examined egr-5 expression levels in greater detail by whole-mount FISH in combination with other post-mitotic markers of the epidermal lineage. Robust egr-5 expression appears to correlate with strong AGAT-1 expression, but egr-5 is also lowly expressed in AGAT-1^neg zpuf-6⁺ mesenchymal cells (Figure 5—figure supplement 2A, top panel). By whole-mount FISH metrics, egr-5 expression is also barely detectable in a subset of prog-1⁺ cells and in the epidermis (Figure 5—figure supplement 2A, bottom panels), demonstrating that egr-5 expression, although varied in signal intensity, spans the domain of all post-mitotic lineage markers. In addition, egr-5 transcripts are predominantly found in the FACS-dissociated Xins population comprised of post-mitotic cells and not in the X1 neoblast dividing fraction (Figure 5—figure supplement 2B), further suggesting that egr-5 mainly functions in post-mitotic cells.

To investigate the role of egr-5 during normal tissue turnover, animals were fed egr-5 dsRNA every 3 days and subsequently screened for gross morphological defects. egr-5(RNAi) animals displayed a range of deformities: loss of epidermal integrity, anterior blebbing, and complete loss of the anterior region (Figure 5A). Moreover, the severity of phenotypic progression leading to animal death caused by lysis is inversely correlated with animal size (Figure 5B). Because these gross morphological defects manifested in egr-5(RNAi) animals suggested problems with epidermal morphology or density, we quantified epidermal nuclear density and found that egr-5(RNAi) animals had a significantly reduced number of epidermal cells compared to control animals (Figure 5C).

Figure 5 with 4 supplements see all

Download asset Open asset

Next, we investigated epidermal morphology by looking at expression of zpuf-6 and rootletin using whole-mount FISH and ultrastructural features by scanning electron microscopy (SEM). We observed a significant reduction in zpuf-6 expression (see more below) in the epidermis between day 12 and day 24 of RNAi treatment, whereas rootletin expression did not diminish, although it appeared disorganized compared to controls (Figure 5—figure supplement 3A–B). Analysis of the ventral cilia by SEM qualitatively showed similar levels in abundance. However, the ventral epidermal cells in egr-5(RNAi) animals appeared more smooth and stretched out and lacked epidermal pores/pits as seen in control animals (Figure 5—figure supplement 3C). The smoothness of the ventral epidermis also contributed to difficulty in delineating individual epidermal cells. Our data suggest that the reduced epidermal density in egr-5(RNAi) animals could be caused by the failure of new epidermal progeny cells to properly differentiate and incorporate into the epidermis, resulting in stretching and apparent morphological defects. Consistent with this hypothesis, we also found that there was a significant reduction in the number of laminB⁺ cells in the anterior region of egr-5(RNAi) animals compared to control animals after 21 days of RNAi treatment (Figure 5—figure supplement 4).

Given that egr-5 knockdown causes a marked reduction in zpuf-6 expression in the epidermis, we looked at our panel of epidermal lineage markers to determine which particular transition may be affected in the course of differentiation. In RNAi conditions such as chd4, p53 and zfp-1 knockdown where AGAT-1⁺ cell populations are lost, there is a distinct progressive disappearance in the spatial distribution of cells, but the signal intensity of the gene marker expressed in the remaining cells is not significantly diminished (see Figure 2—figure supplement 1A; Figure 3—figure supplement 2B). This suggests that these factors are required for the generation or maintenance of AGAT-1⁺ cells, but not necessarily required for the expression of AGAT-1 or other markers in these cells. However, a TF can be responsible for the expression of AGAT-1 or some other marker, and not be required for the maintenance of that cell type, or both. To distinguish between the loss of cell type versus the loss of specific gene expression, we made RNA probe mixes consisting of multiple zpuf genes from our screen expressed in either AGAT-1 or zpuf-6 cells to increase overall signal intensity and detection.

We first investigated the spatial distribution of prog-1 and AGAT-1 progeny markers and found that the stereotypical sub-epidermal expression domains of both markers had expanded deeper into the mesenchyme in egr-5(RNAi) animals (Figure 5D), suggesting potential temporal defects in the earlier steps of epidermal lineage progression. We then simultaneously analyzed prog-1, AGAT-1 and zpuf-6 markers by whole-mount FISH and established that in egr-5(RNAi) animals, epidermal progeny cells in the mesenchyme are not specifically lost, but that AGAT-1 and zpuf-6 expression levels were significantly reduced (Figure 5E–F). Although zpuf-6 mesenchymal expression is very weak after egr-5 knockdown, it is virtually undetectable in the epidermis, potentially highlighting a key transitional defect. Given that egr-5 knockdown reduces both epidermal cell density as well as the expression of AGAT-1 and zpuf-6 markers, we suggest that egr-5 is responsible for both the maturation of epidermal progeny cells and the expression of progeny markers. Interestingly, prog-1 and AGAT-1 normally do not exhibit significant overlapping expression domains, but here they display significant overlapping expression patterns after loss of egr-5 (Figure 5E, left panels). Together, our results demonstrate an essential role for egr-5 in the proper spatial and temporal progression of progeny cells in the epidermal lineage and that misexpression of these marker domains likely causes improper differentiation and failed epidermal maturation.

Although egr-5 is expressed in post-mitotic epidermal lineage cells to coordinate their proper temporal transition states, we wondered if loss of egr-5 may have any non-autonomous effects on neoblast and progenitor dynamics. To examine neoblast division kinetics, we quantified histone H3Ser10 phosphorylation (H3P) during egr-5 RNAi knockdown and observed a significant increase in neoblast proliferation, followed by an eventual decline (Figure 6A–B).

Figure 6 with 1 supplement see all

Download asset Open asset

Because neoblasts constitute a mixed population of pluripotent stem cells and lineage-committed progenitors, we looked more closely at the composition of dividing neoblasts to determine if egr-5 knockdown expanded all dividing progenitors or only affected a specific subset. We picked a time point where H3P levels were high in egr-5(RNAi) animals (before the decline) and then analyzed the proportion of dividing progenitors for major tissues, including the epidermis (zeta-class) (van Wolfswinkel et al., 2014), protonephridia (pou2/3) (Scimone et al., 2011), gut (hnf4) (Wagner et al., 2011) and brain (pax6A) (Wenemoser et al., 2012) (Figure 6C–G, Figure 6—figure supplement 1A). Despite demonstrating a dramatic increase in overall dividing progenitors, egr-5(RNAi) animals displayed similar proportions of dividing tissue progenitors compared to their control counterparts. These results suggest that egr-5 RNAi knockdown causes an overall expansion of multiple lineage progenitors, including the epidermal progenitors (zeta-class). Therefore, failure to generate mature epidermal cells in egr-5(RNAi) animals is not caused by a failure to generate zeta-class epidermal progenitors. Rather, it is likely caused by a defect at the post-mitotic differentiation level. This decrease over time in the integration of functional, differentiated cells leads to the eventual breakdown of epidermal integrity and results in animal lysis.

To assess whether the expanded progenitor population in egr-5(RNAi) animals led to any pronounced defects in differentiation, we examined by colorimetric WISH the protonephrida (slc6a-13) (Vu et al., 2015), epidermal progenitors/prog-1/AGAT-1 (p53), and overall gut morphology (hnf4). Qualitatively, all markers exhibited prominent increases in number and expression in egr-5(RNAi) animals compared to controls (Figure 6—figure supplement 1B–D). We quantified the number of protonephridial proximal units (PU) and found that egr-5(RNAi) animals indeed display a significant increase in protonephridial density (Figure 6H–I). Based on general morphology and the observation that egr-5 RNAi animals did not exhibit defects in osmoregulation by bloating, we conclude that the protonephridia are functional. The strong increase in p53 expression builds upon our previous finding that loss of egr-5 increases the number of zeta-class epidermal progenitors (Figure 6C, Figure 6—figure supplement 1A), as well as prog-1 and AGAT-1 cells (Figure 5D–E). The enriched gut expression marked by hnf4 in egr-5(RNAi) animals also suggests the presence of an increased number of gut cells, without any discernible physiological consequences (e.g., animals can still eat). Together, our data suggest that the extensive increase in neoblast proliferation caused by egr-5 depletion is remarkably counterbalanced by a corresponding increase in the normal progression of multiple lineage progenitors that differentiate and integrate into functional tissues.

The potential non-cell autonomous effect of egr-5 knockdown on neoblast and progenitor dynamics raised the possibility that epidermal defects resulting from abnormal differentiation and defective homeostasis may induce some kind of global stress response. Therefore, we measured apoptosis using a whole-mount TUNEL assay (Pellettieri et al., 2010) and observed a remarkable increase in cell death throughout the course of egr-5 RNAi treatment compared to control animals (Figure 7A). Notably, the increase in cell death appeared to occur uniformly throughout the animal, with the majority of TUNEL-positive nuclei found in the mesenchymal tissue and not in the epidermis (not shown).

A global increase in both cell proliferation and cell death are general features of the initial stages of planarian regeneration (Pellettieri et al., 2010; Wenemoser and Reddien, 2010). Simple wounding and amputations causing loss of tissue are both capable of activating the expression of many wound-induced genes within the first 24 hr of insult (Wenemoser et al., 2012). We hypothesized that the epidermal defects caused by egr-5 RNAi knockdown may trigger the activation of wound-induced genes expressed in differentiated cells. Delta-1 is a putative Notch signaling pathway ligand whose expression is induced in the epidermis between 6-24 hr after injury, not necessarily localized to the site of wounding (Wenemoser et al., 2012). We monitored the expression of delta-1 in intact, uninjured egr-5(RNAi) animals and observed a marked increase in delta-1 in the anterior region of animals compared to their control counterparts through 21 days of RNAi treatment (Figure 7B, top row). We also measured the in vivo expression of other ‘immediate early genes’ that are activated after wounding including fos-1, jun-1 and egr-3 (Wenemoser et al., 2012), and observed notable increases in the anterior region of egr-5(RNAi) animals (Figure 7B, middle rows). Interestingly, the TGF-β inhibitor follistatin, whose expression is specifically induced for multiple days after a loss of animal tissue (Gavino et al., 2013), does not appear to be stimulated over the course of egr-5 RNAi treatment (Figure 7B, bottom row). In addition, the immediate early genes fos-1, jun-1 and egr-3 are enriched in our chd4 and p53 RNAi whole-worm RNA-seq up-regulated datasets (Supplementary file 2; Figure 7—figure supplement 1). Taken together, these data suggest that both defects in epidermal differentiation (egr-5 RNAi) and loss of epidermal progeny cells (chd4 and p53 RNAi) likely contribute to a breach in epidermal integrity, which consequently may activate a systemic wound response program.

Figure 7 with 1 supplement see all

Download asset Open asset

Zeta-class neoblasts, marked by zfp-1, were recently reported to represent a class of epidermal progenitor cells, giving rise to prog-1⁺, AGAT-1⁺ mesenchymal cells, and other markers of epidermal cells (van Wolfswinkel et al., 2014). However, their precise molecular relationships and the mechanisms underlying the differentiation of these progeny cells are poorly understood.

Our results extend upon these studies by implementing a whole worm RNA-seq approach to identify genes expressed in AGAT-1⁺ cells and their progeny and to characterize the post-mitotic transition states downstream of epidermal progenitors. We utilized a combination of assays, including whole-mount FISH, spatiotemporal kinetics after irradiation, and BrdU lineage tracing, to establish the temporal order of epidermal cell differentiation. Combinatorial FISH revealed that zpuf-6 is co-expressed at varying levels in virtually all AGAT-1⁺ cells, additional AGAT-1^neg mesenchymal cells (AGAT-1^negzpuf-6⁺), and cells evenly distributed throughout the dorsal and ventral epidermis (Figure 2). We examined the expression patterns of all other genes from our screen that displayed discrete mesenchymal cell patterns in combination with AGAT-1 and zpuf-6, and found that they all overlapped extensively. Many of these genes, such as pla2 and odc-1, exhibited moderate expression in AGAT-1^neg zpuf-6⁺ cells, but expression was noticeably reduced in zpuf-6⁺ epidermal cells (Figure 2—figure supplement 3), suggesting that these cells are turning off expression of those genes. In addition, AGAT-1 and zpuf-6 showed conserved spatiotemporal down-regulated kinetics after irradiation and in RNAi conditions leading to loss of AGAT-1⁺ cells (Figure 3, Figure 3—figure supplement 1,2), and a pulse of BrdU incorporated into dividing neoblasts showed a temporal progression through prog-1⁺, AGAT-1⁺ and zpuf-6⁺ cells (Figure 4A). Altogether, we interpret these observations to mean that prog-1, AGAT-1 and zpuf-6 are markers of three major spatially and temporally related cell populations representing distinct transitional stages of epidermal maturation (Figure 8).

Figure 8

Download asset Open asset

Differentiation and migration toward the mature planarian epidermis are temporally correlated, but what are the fates of cells once they pierce through the basement membrane and integrate into the epidermal layer? We propose that zpuf-6⁺ epidermal cells may represent a stable transition state with the potential to give rise to an unknown number of differentiated cell types (Figure 8). Zpuf-6⁺ epidermal cells are interspersed throughout the epidermis but do not express high levels of a cilia gene, rootletin (Figure 4C), which is a marker for differentiated ciliated epidermal cells. However, zpuf-6⁺ epidermal cells overlap extensively with a predicted intermediate filament (IF) gene, vim-3, which is expressed more broadly in the dorsal epidermis. Zpuf-6^neg epidermal cells express lower levels of vim-3 than zpuf-6⁺ cells, but these cells begin in turn to express rootletin (Figure 4C). IF proteins are a diverse class of cytoskeletal elements that provide cells with mechanical support, contribute to the structural integrity of cells, and act as markers of morphogenesis and differentiation (Kim and Coulombe, 2007). In addition, zpuf-6⁺ cells overlap with two specific subset of cells: one expressing the marker NB.22.1E surrounding the ventral mouth opening and the other expressing the laminB IF gene that marks adhesive cells along the animal body margin (Figure 4—figure supplement 3). Lastly, zpuf-6⁺ cells are present at higher cell density in the regenerating blastema than under normal homeostasis conditions (Figure 4—figure supplement 2), suggesting that they have not completed differentiation.

The decisions influencing a particular cell fate path for zpuf-6⁺ epidermal cells could be determined by the immediate, surrounding environment. Given that dorsal and ventral epidermal cells are morphologically distinct, positional cues and signals are likely to contribute to the cell fate specification process. For example, transplantation experiments involving a plug of tissue from a donor transplanted in the reverse dorsal-ventral (DV) orientation with regard to the host can trigger tissue outgrowths on both sides where the grafted tissue maintains its original DV identity (Kato et al., 1999). The sub-epidermal layer of body wall muscle cells express most of the polarity cues that organize the planarian body axes, including receptors, ligands and inhibitors of the BMP and Wnt signaling pathways (Witchley et al., 2013), and could potentially impact epidermal fates as well. The role of Notch signaling has not been well characterized in planarians but has been implicated in the regulation of cell fate determination in epithelial cells of the mouse intestine (VanDussen et al., 2012) and the lung (Rock et al., 2011; Rock and Hogan, 2011). We cannot exclude the possibility that zpuf-6⁺ epidermal cells simply reflect an abundant, differentiated cell population that undergoes rapid turnover. Additionally, there could be other progenitors (zeta-class dependent or independent) that give rise to specific epidermal cell types that have not been defined molecularly. To build upon our understanding of epidermal cell fate specification, it will be necessary to define all the types and numbers of cells that ultimately compose both the ventral and dorsal planarian epidermis, as well as the factors that contribute to their identity and maintenance.

In planarians, TFs play prominent roles in specifying neoblasts to a committed fate (Reddien, 2013). Key TFs that are important for specifying neoblasts during homeostasis and regeneration often remain expressed in mature cell types, suggesting that these TFs are also required for the maintenance of cellular identity. Studies of planarian eye regeneration currently provide the most detailed characterization of planarian lineage specification. During regeneration, eye progenitors express conserved TFs ovo, six-1/2 and eya at a distance from the eye primordium, forming a trail of cells and undergoing changes in gene expression as they migrate toward the primordium, where differentiation markers begin to be expressed (Lapan and Reddien, 2011; Lapan and Reddien, 2012). However, ovo, six-1/2 and eya are also expressed in mature eyes and are likely involved in the maintenance of the optic cup cells and photo-sensing neurons of this sensory organ.

Our findings also suggest that distinct molecular changes occur both spatially and temporally, and are necessary for the proper differentiation of epidermal progenitors into mature cell types. However, egr-5, which is required for the stable maturation of epidermal progenitors and their progeny, does not appear to be required for the maintenance of the differentiated state. In fact, the lack of egr-5 expression in the epidermis suggests that shutting off its expression may be an important step in promoting differentiation. Whether this regulation is intrinsic or extrinsic awaits further investigation. Moreover, because loss of egr-5 leads to the misexpression of AGAT-1 and zpuf-6 markers in prog-1⁺ cells, it could be directly or indirectly involved in turning off prog-1 expression and activating the expression of genes required for the AGAT-1⁺ transition stage. A key remaining question is how individual epidermal progeny cells undergo distinct differentiation transition stages, including responses to external signals or potential changes in chromatin structure.

Genes involved in small molecule metabolic processes (specifically organonitrogen processes) were enriched in our common chd4 and p53 RNAi down-regulated gene list. Many of the metabolic genes screened by in situ were expressed in AGAT-1⁺ and zpuf-6⁺ cells, including gatm, odc-1, cytochrome p450s (cyp), and pla2. Furthermore, genes including acsl-2, adss, slc25a-19 and zpuf-7 were expressed in both the gut and cells of epidermal lineage, which suggests that these tissues may share common functional roles.

There are four arginine:glycine amidinotransferases (AGAT-1, AGAT-2, AGAT-3, gatm) identified in planarians. These enzymes catalyze the first step in creatine biosynthesis by transferring the amidino group of arginine to glycine to yield ornithine and guanidinoacetic acid (Wyss, 2000). Creatine plays an important role in muscle energy homeostasis. In vertebrates, de novo synthesis of creatine mainly takes place in the kidney, pancreas and liver (Nabuurs et al., 2013). Thus, it is assumed that AGAT-1⁺ cells synthesize creatine, which is then released and taken up by neighboring muscle cells and neurons (Eisenhoffer et al., 2008), though this has not been formally demonstrated.

Polyamines are ubiquitous polycations that are essential for eukaryotic cell growth, and polyamine metabolism is frequently dysregulated in cancer (Pegg, 2009). Some of the elucidated roles for polyamines in cell growth include maintenance of chromatin conformation, gene regulation, ion channel regulation, and free radical scavenging (Casero and Marton, 2007). Ornithine decarboxylase (ODC) catalyzes the first rate-limiting step in polyamine biosynthesis by converting ornithine to putrescine, which then is converted to spermidine and spermine. In planarians, ODC activity has been reported to be induced early in regeneration near the wound site, and inhibition of ODC activity appears to cause a reduction in cell proliferation and differentiation (Saló and Baguna, 1989). However, mechanisms of action by odc-1 and the regulation of neoblasts require further characterization.

The cytochrome P450s (CYPs) are enzymes that use molecular oxygen and complex reaction chemistry to modify their substrates. They are involved in a large number of physiological processes, including steroid hormone synthesis and detoxification of xenobiotics (Anzenbacher and Anzenbacherova, 2001; Zhang and Yang, 2009). CYPs make up one of the most diverse eukaryotic gene families, and are classified into clans, families and sub-families based on phylogenetics and sequence identity (Nelson et al., 2013). We identified a number of cyp genes that were expressed in AGAT-1⁺ progeny cells as well as in the gut. Steroid metabolism and detoxification are two processes generally performed by the liver, raising the questions of what physiological roles these various planarian cyp genes play and whether AGAT-1⁺ cells may have an additional endocrine-like function.

Although p53 is best known for its central role as a tumor suppressor, it also regulates several aspects of cellular metabolism, including autophagy, central carbon metabolism, and lipid metabolism (Berkers et al., 2013). Fatty acids (FA) have essential roles in the cell as sources of energy, membrane components, and signaling molecules (Lopes-Marques et al., 2013). Long-chain fatty-acid coA-ligases (ACSL) are key enzymes involved in the initial steps of FA metabolism. We identified three enzymes involved in fatty acid metabolism, acsl-1, acsl-2 and pla2, that are expressed in the gut as well as in epidermal progeny cells and the mature epidermis. Taken together, the presence of multiple metabolic genes co-expressed in AGAT-1⁺ cells and their progeny suggests that despite these cells being in a temporal stage transitioning to mature epidermal cells, they may also play a dual functional/physiological role in planarian homeostasis.

Epithelia are a hallmark of multicellular organisms, forming the interface between the organism and the external environment and lining the cavities and surfaces of organs (Donati and Watt, 2015). The planarian epidermis provides the first response to amputation-induced injury by covering up the wound site within 30 min (Reddien and Sánchez Alvarado, 2004). It has been postulated that rhabdites are released at wound sites, where their contents produce a protective mucosal covering, possibly providing immunological functions (Reisinger and Kelbetz, 1964). Therefore, a breach in epidermal integrity could be the stimulus for initiating regeneration, though how this signal could be transduced to neoblasts is unknown.

The egr-5(RNAi) phenotype resulting in the expansion of neoblasts and multiple progenitor populations reveals that complex dynamics are at play between differentiating progeny and how neoblasts respond to loss of tissue integrity. It will be interesting to investigate the neoblast response in other RNAi or pharmacological contexts that directly affect epidermal integrity to determine how these perturbations affect neoblast function. However, it still remains possible that neoblast expansion could be a response to the global increase in apoptosis, although we still do not understand how the processes of cell proliferation and cell death are integrated. Alternatively, egr-5 may be responsible for the expression of a signaling molecule to which neoblasts respond in a feedback control process. Feedback signaling from differentiating progeny to neoblasts may be an important and robust strategy to control stem cell activity during homeostasis and regeneration. Therefore, unveiling the downstream targets of egr-5 will help discern these possibilities, as well as provide more mechanistic insight on the role of egr-5 and the progression of epidermal fate maturation.

Several key questions emerge from our findings. Why does epidermal lineage progression require multiple transition states? Do neoblast progeny contribute to the dynamics of the microenvironment, thereby playing a niche-like role? How are cell autonomous and non-cell autonomous mechanisms coupled for the maintenance of cellular equilibrium at the organismal level? How is the overall physiological status of the animal sensed, integrated and converted to the appropriate neoblast output? With a more nuanced understanding of their adaptive behavior, planarian neoblasts and the epidermal lineage have emerged as an important model system in which to study the dynamics between stem cell self-renewal and the orchestration of progeny differentiation.

Schmidtea mediterranea CIW4 asexual strain was maintained in 1X Montjuic salts supplemented with 50 μg/ml Gentamicin and fed homogenized calf liver paste as previously described (Gurley et al., 2008; Reddien et al., 2005a). Animals were on average starved between 7–10 days prior to starting experiments and ranged in size from 2 mm to 8 mm. Animals were exposed to 6000 Rads of gamma irradiation using a GammaCell 40 Exactor irradiator.

Genes in this study were cloned from a CIW4 cDNA library into pPR-T4P vector (J. Rink) as described elsewhere (Adler et al., 2014). Primer sequences are provided in Supplementary file 4. Cloned gene vectors were transformed into bacterial strain HT115 for dsRNA production. RNAi food was prepared by mixing 50 ml of pelleted culture with 250 μl of calf liver paste (2X) or 50 ml of culture with 125 μl of calf liver paste (4X). Animals were fed every 3 days, with the first day designated as Day 0 of RNAi treatment. The number of RNAi feedings performed for egr-5 knockdown varied depending on starting-size of animals and RNAi food concentration. RNAi feedings and time points analyzed are noted in figure legends as XFdY (3Fd18 = 3 RNAi feedings, day18). For all RNAi feeding experiments, unc22 dsRNA was used as the control for the same number of feedings as experimental RNAi animals. Images of live animals were captured using a Leica M205 stereoscope.

Whole-mount colorimetric and fluorescent in situ hybridizations were performed using a detailed protocol as previously described (King and Newmark, 2013; Pearson et al., 2009). Fluorescence-labeled animals were mounted in ScaleA2 solution (Hama et al., 2011). Immunostaining with anti-H3P (1:1000, Millipore, Billerica, MA) was performed following fluorescent in situ development and was detected using Alexa-conjugated secondary antibody (1:1000, Abcam, Cambridge, MA). Animals were fixed and stained for TUNEL using a method previously described (Pellettieri et al., 2010) with modifications: animals were bleached in 0.075% ammonia and 3% hydrogen peroxide and treated with ProteinaseK (2 μg/ml) in PBSTx (0.3% Triton) for 10 min followed by 4% formaldehyde incubation for 10 min prior to TdT reaction.

Cryosectioned animals were processed after fluorescent in situ development as previously described (Tu et al., 2012). For nuclear fast red staining, animals were processed for colorimetric ISH and were subsequently fixed overnight in 4% paraformaldehyde (PBS) at 4°C followed by dehydration in 30, 50, and 70% ethanol. Fixed specimens were embedded in paraffin and serial sectioned at 10 μm thickness and counter stained with nuclear fast red.

For scanning electron microscopy, animals were immersed in a relaxant fixative (1% HNO₃, 0.85% formaldehyde, 50 mM MgSO₄) as described elsewhere (Rompolas et al., 2013) for 5 min and were replaced with fresh fixative and rocked overnight at room temperature. Animals were then transferred to a solution containing 2.5% glutaraldehyde, 2% paraformaldehyde, 1% sucrose, 1 mM CaCl2 in 0.05 M NaCacodylate buffer pH 7.36 and left at 4°C until ready for processing. Animals were rinsed in ultrapure water, and secondary fixation was performed at 4°C overnight in 2% aqueous osmium tetroxide. Samples were dehydrated in a graded series of ethanol and dried in a Tousimis Samdri-795 critical point dryer. Samples were mounted on stubs and sputter coated with gold palladium. Imaging was done with a Hitachi TM-1000 tabletop SEM.

Colorimetric WISH images were acquired using a Zeiss Lumar V12 stereomicroscope equipped with an AxioCam HRc. Colorimetric WISH tissue sections were imaged using a Zeiss Axiovert206. Fluorescent images were acquired with either a Zeiss LSM-510 VIS or a Perkin Elmer Ultraview VOX spinning disk. Stitching and batch processing of images, and H3P quantifications per surface area were performed as previously reported (Adler et al., 2014). For H3P/progenitor colocalization thresholding, spots were selected based on signal intensity a multiple above the noise background and filtered based on size. For fluorescence intensity quantification, Z-stack images were acquired and z-projected by average in a fixed neighborhood around the bottom of the animal as determined by DAPI nuclear staining after background subtraction. For epidermal cell quantifications, single focal plane images were acquired and nuclei within a fixed sized rectangle in regions of uniform cell density were quantified using ‘Find Maxima’ with Fiji software. All macros and plugins are available at https://github.com/jouyun.

BrdU was administered by soaking animals for 24 hr in 20 mg/ml BrdU (Sigma, St. Louis, MO) dissolved in 3% DMSO and 1X Montjuic. Animals were then incubated in 5 g/L Instant Ocean supplemented with 50 μg/ml Gentamicin over the course of BrdU chase period. Animals were fixed and processed using the in situ hybridization protocol except bleached in 6% H₂O₂ in PBSTx (0.3% Triton) for 3–4 hr under direct light. After development, specimens were treated with 2N HCl for 45 min at room temperature. BrdU was detected using a rat anti-BrdU antibody (1:1000; Abcam, Cat. No. ab6326). Primary antibody was detected with HRP-conjugated anti-rat antibody (1:1000; Jackson ImmunoResearch, West Grove, PA). Images were acquired with a Perkin Elmer (Waltham, MA) Ultraview VOX spinning disk with a 20x 0.8 NA Plan Apochromat objective (Zeiss, Oberkochen, Germany) onto an Orca R2 (Hamamatsu Photonics, Hamamatsu, Japan) camera. Integration times were adjusted per worm to achieve well-saturated images. Z-stacks were collected at the anterior end of the worm over 250 µm in the Z dimension, with 2 µm steps, which contained nearly the entire depth of the worm. Images were cropped to a 400 pixel² (132 µm²) region posterior to the brain. The stack was also cropped in Z to the ventral half of the stack. Target-specific cells were first marked manually in Fiji followed by target and BrdU double-positive cells. At least 6 worms (up to 10) per time point were analyzed and roughly 6,500 cells total were quantified. The average percent positive cells were reported ((double-positive/total)*100). Representative images (Figure 4A) were first background subtracted in ImageJ with a 25 pixel radius rolling ball and then gaussian-blurred with a 1 pixel radius.

For RNA-seq analysis of control, chd4(RNAi) and p53(RNAi) animals, three biological replicates of 10 worms each were collected for RNA isolation. mRNAseq libraries were generated from 500ng of high quality total RNA, as assessed using the LabChip GX (Perkin Elmer). Libraries were made according to the manufacturer’s directions for the TruSeq Stranded mRNA LT– set A and B (Illumina, San Diego, CA; Cat. No. RS-122-2101 and RS-122-2102). Resulting short fragment libraries were checked for quality and quantity using a LabChip GX and Qubit Fluorometer (Life Technologies, Carlsbad, CA). Libraries were pooled, requantified and sequenced as 50 bp single reads on the Illumina HiSeq 2500 instrument using HiSeq Control Software 2.0.5. Following sequencing, Illumina Primary Analysis version RTA 1.17.20 and Secondary Analysis version CASAVA-1.8.2 were run to demultiplex reads for all libraries and generate FASTQ files. RNA-seq analysis was carried out by mapping sequence reads to a set of 36,035 S. mediterranea transcripts assembled from various sources including a previous transcriptome used in microarray studies (Adler et al., 2014), trinity assemblies from lab-generated data involving whole animals, embryos, and sorted X1 cells, a transcriptome from the Bartscherer lab (Boser et al., 2013), and the Dresden transcriptome assembly from PlanMine (http://planmine.mpi-cbg.de), reduced as a collection to unique representations of loci via CD-HIT (Fu et al., 2012). Sequences can be downloaded from http://smedgd.stowers.org. Reads were aligned using bowtie with the following parameters: --best --strata -v 2 -m 5, and read counts to genes were tallied from the SAM files with a custom script. Differential gene expression was evaluated using R and the edgeR library (Robinson et al., 2010). P-values were adjusted as previously described (Benjamini and Hochberg, 1995). The RNA-Seq data and the transcriptome against which it was quantified have been archived at GEO and is available under accession number: GSE72389. Gene Ontology (GO) (Gene Ontology Consortium, 2015) terms were assigned to each S. mediterranea gene based on homologous PFAM domains and significant Swissprot hits. GO term enrichment was performed using the R package topGO (Alexa and Rahnenfuhrer, 2010).

The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication.

We would like to thank all members of the Sánchez Alvarado lab for their support, especially Carrie Adler and Erin Davies for helpful comments on the manuscript. We thank Beth Duncan for kindly sharing X1 data and Sarah Elliott and Aurimas Gumbrys for RNA probes. We thank Eric Ross and Kirsten Gotting for bioinformatics assistance, Jay Unruh for statistical assistance, David Nelson (University of Tennessee Health Science Center) for planarian CYP nomenclature, and Mark Miller for illustration assistance. We also acknowledge Melania McClain for SEM assistance, Allison Peak for RNA-seq library assistance and all other members of the Histology, Microscopy, Molecular Biology and Planaria core facilities at the Stowers Institute for their technical support. This work was supported by NIH R37GM057260 to ASA. KCT was an Ellison Medical Foundation/AFAR Fellow of the Life Sciences Research Foundation. ASA is an investigator of the Howard Hughes Medical Institute.

© 2015, Tu et al.

This article is distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use and redistribution provided that the original author and source are credited.