Comparative transcriptomic analyses and single-cell RNA sequencing of the freshwater planarian Schmidtea mediterranea identify major cell types and pathway conservation

Background In the Lophotrochozoa/Spiralia superphylum, few organisms have as high a capacity for rapid testing of gene function and single-cell transcriptomics as the freshwater planaria. The species Schmidtea mediterranea in particular has become a powerful model to use in studying adult stem cell biology and mechanisms of regeneration. Despite this, systematic attempts to define gene complements and their annotations are lacking, restricting comparative analyses that detail the conservation of biochemical pathways and identify lineage-specific innovations. Results In this study we compare several transcriptomes and define a robust set of 35,232 transcripts. From this, we perform systematic functional annotations and undertake a genome-scale metabolic reconstruction for S. mediterranea. Cross-species comparisons of gene content identify conserved, lineage-specific, and expanded gene families, which may contribute to the regenerative properties of planarians. In particular, we find that the TRAF gene family has been greatly expanded in planarians. We further provide a single-cell RNA sequencing analysis of 2000 cells, revealing both known and novel cell types defined by unique signatures of gene expression. Among these are a novel mesenchymal cell population as well as a cell type involved in eye regeneration. Integration of our metabolic reconstruction further reveals the extent to which given cell types have adapted energy and nucleotide biosynthetic pathways to support their specialized roles. Conclusions In general, S. mediterranea displays a high level of gene and pathway conservation compared with other model systems, rendering it a viable model to study the roles of these pathways in stem cell biology and regeneration. Electronic supplementary material The online version of this article (10.1186/s13059-018-1498-x) contains supplementary material, which is available to authorized users.


Background
Investigations using model organisms such as Caenorhabditis elegans, Drosophila melanogaster, zebrafish, and mice continue to drive fundamental insights into the molecular mechanisms driving a variety of conserved biochemical processes [1]. However, much attention has recently turned to the use of non-traditional organisms as models to explore more specialized pathways. For example, while freshwater planarians (flatworms) have been used in a laboratory setting for more than 100 years due to their ability to regenerate following virtually any injury, the planarian Schmidtea mediterranea has emerged as a powerful model for dissecting the molecular basis of tissue regeneration [2,3]. Despite significant resources put forth to develop S. mediterranea as a model in the lab, systematic genome-scale investigations of gene function and conservation are lacking.
Much of the interest in planarians is driven by the fact that approximately 20% of their adult cells are stem cells (called neoblasts), at least some of which are pluripotent [4][5][6][7]. In addition, planarians are one of the only models that can be used to rapidly test gene function in adult animals through RNA interference (RNAi) screening. Placing gene function in an evolutionary context is critical not only to inform on the conservation of pathways related to stem cell biology and regeneration, but also because planarians represent a key member of the otherwise neglected superphylum Lophotrochozoa/Spiralia (subsequently referred to as Lophotrochozoa), and they can further be used to model closely related parasitic flatworm species (e.g., flukes and tapeworms), which infect an estimated hundreds of millions worldwide [8].
In attempts to complement ongoing genome sequencing efforts [9,10], several transcriptome datasets have been generated for S. mediterranea under various physiological conditions using a variety of experimental techniques [11][12][13][14][15][16][17][18]. In isolation, each set provides a snapshot of planarian gene expression under a specific condition; however, recent efforts have focused on integrating several transcriptomes to generate a more comprehensive overview of gene expression [9,19]. The SmedGD repository was generated by integrating transcriptomes from whole-animal sexual and asexual worms, whereas the PlanMine database serves as a repository for the published genome as well as existing transcriptomes from the community to be deposited and queried. However, they lack systematic and comparative evolutionary and functional genomics analyses, which are required for understanding the mechanistic basis of biological processes. Together these datasets comprise more than 82,000 "transcripts" with little assessment of "completeness" from an evolutionary perspective.
Typically, transcriptome datasets are generated from entire organisms or tissues [20][21][22]; however, such analyses can mask the contribution of specific cell subpopulations, which can be particularly problematic when attempting to elucidate, for example, pathways expressed during key cellular events. While cell sorting offers the capability to enrich for specific cell subpopulations, the emergence of single-cell RNA sequencing (scRNAseq) offers a powerful route for interrogating gene expression profiles from individual cells [23,24]. Applied to S. mediterranea, this technology is expected to yield molecular-level insights into the roles of distinct cell types, such as neoblasts, during homeostatic tissue maintenance and regeneration [7,[25][26][27]. Indeed, scRNAseq experiments have already been used to resolve neoblast heterogeneity and identify regulators of lineage progression [26][27][28][29][30].
In this study, we generate a high-confidence transcriptome pruned from an integrated transcriptome generated earlier in the lab [18], which, through combining transcriptomes from diverse physiological conditions and experimental techniques, leads to a large number of transcripts (n = 83,469) for S. mediterranea. Next, we apply systematic bioinformatic approaches to annotate and compare the complement with model organisms and other Platyhelminthes. This pipeline predicts putative functional annotations of the transcriptome, identifying a set of transcriptionally active transposons as well as extended families of cadherins and tumor necrosis factor (TNF) receptor associated factor (TRAF) proteins. Metabolic reconstruction further reveals an increased biochemical repertoire relative to related parasitic platyhelminths. In order to gain insights into the role of these pathways in planarian biology, high-throughput scRNAseq was performed, capturing the transcriptional signatures from~2000 cells. From the 11 distinct clusters of transcriptional profiles, we identified clusters corresponding to neoblasts, epithelial progenitors, muscle, neurons, and gut, among which neoblasts exhibit the most metabolically active profiles. We also identify a novel cluster: a cathepsin + cluster representing multiple unknown mesenchymal cells. Beyond giving us new insights into the evolution and dynamics of genes involved in regenerative pathways, the data and analyses presented here provide a complementary resource to ongoing genome annotation efforts for S. mediterranea. They are available for download from http://www.compsysbio.org/ datasets/schmidtea/.
Next, the protein-coding potential of the remaining transcripts was predicted using the error-tolerant ESTScan [39]. Finally, transcripts without matches to the above were parsed through a six-frame translation algorithm to identify the largest potential open reading frame (LongestORFs). ESTScan and LongestORFs predictions were further filtered such that only those predicted to have > 100 amino acid residues and also to co-localize on the genome with known S. mediterranea transcripts derived from complementary resources (EST database of the NCBI, SmedGD v2.0 [9] and the Oxford dataset [14]) were included in our final filtered dataset (Fig. 1a, b).
Together, this filtered set comprises 36,026 sequences, of which 28,583 map to 22,215 loci of the S. mediterranea genome assembly deposited in SmedGD v2.0 [9]; the remaining 7443 sequences could not be mapped. Of these unmapped transcripts, 1008 share significant sequence similarity, i.e., ≥ 80% sequence identity as assigned by the Basic Local Alignment Search Tool (BLAST) [40], with a known S. mediterranea protein, 106 to a protein from the closely related planarian D. japonica, and 65 to proteins from other Platyhelminthes. Such matches indicate that these sequences are likely bona fide transcripts that are missing from the current S. mediterranea genome assembly. Interestingly, among the 7443 unmapped transcripts, we also identified 794 with significant sequence identity (≥80% sequence identity as assigned by BLAST) to a non-metazoan protein in the UniProt database. Among these were 728 sequences matching sequences from Tetrahymena thermophila and a further 22 [10] and applying the F1 cutoff defined by the Spaln alignment tool (corresponding to~73% sequence identity and~73% coverage) [41] resulted in mapping 33,487 transcripts (~95% of the transcriptome) to 20,483 genomic positions (Fig. 1c, Additional file 2: Figure S1A). In contrast, using similar parameters resulted in the mapping of 38,186 PlanMine transcripts (~91.5% of the transcriptome) to 26,510 positions. Of these, 31,286 (~89%) Toronto transcripts overlap with 33,191 PlanMine transcripts (79.5%), corresponding to 14,145 positions. Although both transcriptomes map a substantial proportion of their transcriptomes to the reference genome, PlanMine maps a higher number of transcripts. However, it is noteworthy that the Toronto transcriptome contributes 2231 transcripts (~6%) that exclusively map to the reference genome. Interestingly, while PlanMine and Toronto transcripts that map to the same loci are of similar length, PlanMine transcripts that are either unmapped or map to unique regions are significantly longer than the equivalent Toronto transcripts (Additional file 2: Figure S1B). Analyzing the distribution of sequence similarity bit scores further reveals that the unmapped transcripts from both the Toronto and PlanMine transcriptomes consist of many high-scoring matches, suggesting their likely validity (Additional file 2: Figure S1C).
Comparisons with three previously generated transcriptomes: SmedGD v2.0 (n = 22,855, [9]), PlanMine (n = 41,475, [19]), and Oxford (n = 23,545, [14]), revealed a core set of 24,477 transcripts common to all four sets, together with 1820 transcripts unique to the Toronto set (defined as those with bit score < 40 for BLASTn [40] searches using a relaxed word size of 7 in order to maximize sensitivity); Fig. 1d). Of the unique transcripts, 371 (20.3%) share significant sequence similarity (BLAST, E-value <1e-08, % sequence identity ranging from 1.5% to 100%) to known proteins in UniProt and 1427 (78%) represent ESTScan predictions. Supporting the validity of these unique transcripts, we note that 1399 (~74%) map to the latest PlanMine genome dd_Smes_g4 [10]. To further assess transcriptome completeness, we performed a systematic comparison with the core eukaryotic and metazoan gene sets defined by BUSCO v1 [38], demonstrating that our high-quality transcriptome exhibits similar coverage (81% eukaryotic, 78% metazoan) as PlanMine (81% eukaryotic, 78% metazoan) and higher coverage than the Oxford (78% eukaryotic, 73% metazoan) and SmedGD (62% eukaryotic, 50% metazoan) datasets ( Fig. 1e). Additionally, the Toronto transcriptome features a lower fraction of partially recovered transcript sets. However, it is noteworthy that of the 348 BUSCO genes, representing single-copy genes from 310 different eukaryotes that were completely recovered by the Toronto dataset, 86 appear to possess paralogs in the Toronto dataset as compared to 112 in PlanMine. Such duplicates might represent either errors during transcript assembly or alternative spliceoforms.
Functional annotation of S. mediterranea proteome: expanded set of transposons and TRAFs Having compiled and validated a high-confidence set of transcripts, we next analyzed functional potential through a systematic annotation of protein domains inferred by the InterPro resource [36]. Gene Ontology (GO) assignments [42,43] based on domain annotations of predicted proteins revealed that transport, signal transduction, biosynthetic process, cellular nitrogen compound metabolic process, and cellular protein modification process are the five most abundant biological processes, consistent with other eukaryotes (Additional file 2: Figure S2).
To identify taxon-specific gene family expansions in S. mediterranea, we compared the 20 most abundant Pfam [44] annotations of predicted protein sequences in our dataset to the proteomes of Homo sapiens, Drosophila melanogaster, and Caenorhabditis elegans, as well as several parasitic flatworms for which genome sequence data are available (cestodes: E. granulosus, E. multilocularis, T. solium, H. microstoma; trematodes: Schistosoma mansoni, S. haematobium, C. sinensis, O. viverrini; monogeneans: Gyrodactylus salaris) (Fig. 2a). Consistent with the other metazoans, the most abundant domains are Pkinase (PF00069), 7tm (PF00001), and Ank (PF12796). Among the remaining 17 abundant domains, three represent lineage-specific expansions: transposase-related domains, DDE_1 (PF03184) and DDE_Tnp_1_7 (PF13843) (ranked 4th and 9th most abundant, respectively)which are significantly expanded only in S. mediterranea and not in other Platyhelminthesand the meprin and TRAF homology (MATH) domain (PF00917, ranked 8th most abundant)expanded in S. mediterranea in comparison to other Platyhelminthes. Another domain of interest is the cadherin domain (PF00028, ranked 16th most abundant), which is expanded throughout Platyhelminthes and also in humans, suggesting a more fundamental role for this domain.
Although S. mediterranea exhibits a larger (n = 290) repertoire of the transposase-related domains, DDE_1 and DDE_Tnp_1_7, relative to other helminths (Fig. 2a), the transcripts associated with these domains are expressed at relatively low levels: mean reads per kilobase per million mapped reads (RPKM) 1.22 +/− 0.04 and 1.10 +/− 0.42 for DDE_1 and DDE_Tnp_1_7, respectively; bottom 40% of expressed transcripts (Additional file 1). Transposable elements (TEs, sequences which can change position within a genome) are classed into two types: class I (retrotransposons), which operate via a copy-and-paste mechanism and include long and short interspersed nuclear elements (LINEs and SINEs, respectively), and class II (DNA transposons), which operate via a cut-and-paste mechanism [45]. DNA transposons are the most abundant elements for transcripts with both DDE_1 and DDE_Tnp_1_7 domains. To determine whether these elements may be functionally active in the S. mediterranea genome, we estimated the sequence divergence of each copy relative to the consensus (Fig. 2b, [46]). Of 1641 elements, we found that 180 (13%) of DDE_1 domains and 97 (25%) of DDE_Tnp_1_7 domains exhibit relatively low sequence divergence (< 5%), indicating that they may still be functionally active. Among DDE_1 domain transcripts, almost half represent the TcMar-Tigger element, thought to be a distant relative of Mariner [47], while for DDE_Tnp_1_7 domain transcripts, the majority represent the PiggyBac element.
Beyond transposons, we found that the MATH (121 domains) domain represents S. mediterranea-specific expansions. MATH domains are present in mammalian tissue-specific metalloendopeptidases (meprins) and TNF receptor associated factor (TRAF) proteins. BLAST searches of MATH-domain-associated proteins in S. mediterranea suggest they are likely TRAF proteins (Additional file 1), important regulators of signal transduction, cell death, and cellular responses to stress [48], immune response [49], and cellular degradation [50]. Many of these domains contain transcripts that are expressed at relatively high levels (mean RPKM 18.05 +/− 5.39; top 20% of expressed transcripts; Additional file 1), suggesting an important regulatory role. Another gene family with abundant representation in Platyhelminthes is the cadherins. Cadherin-domain-containing transcripts were moderately expressed (mean RPKM 4.22 +/− 1.26; top 40% of expressed transcripts; Additional file 1). Cadherins are transmembrane proteins involved in regulating cell-cell adhesion, morphogenesis, and cell recognition [51,52]. More than 100 cadherins have been characterized in vertebrates,  Figure S3). This latter cluster corresponds to calsyntenins (CLSTN), calcium-binding type I transmembrane proteins belonging to the cadherin superfamily, predominantly expressed in neurons. This cluster contains sequences from human (CLSTN1, CLSTN2), C. elegans (CASY-1), C. sinensis, O. viverini, and S. mediterranea (Smed-calsyntenin -SmedASXL_013539). Consistent with its expression in neurons in other organisms, Smed-calsyntenin is predominantly expressed in the brain and ventral nerve cords (with weaker expression detected in the gut), and it exhibits a high degree of co-localization with the cholinergic neuron marker chat (Fig. 2d, e). In the future it will be interesting to determine whether the expansion of TRAF proteins in comparison to other parasitic flatworms and the abundance of cadherins in S. mediterranea represent increased functional complexity in signal transduction and regeneration in planarians.

S. mediterranea expresses a diverse repertoire of transcription factors
We next investigated the repertoire of transcription factors in S. mediterranea in the context of other eukaryotes. Transcription factors were predicted for S. mediterranea, together with an additional 165 eukaryotes [53]. Our predictions suggest that 843 S. mediterranea transcripts encode transcription factors associated with 55 classes ( GATA, have been studied in S. mediterranea [54,55], several others with high abundances in S. mediterranea and vertebrates remain poorly characterized. These include CSD (cold-shock domain; involved in transcriptional repression and activation and in mRNA packaging, transport, localization, masking, stability, and translation) and bZIP_maf (acting as key regulators of terminal differentiation in many tissues, such as bone, brain, kidney, lens, pancreas, and retina, as well as in blood). These transcription factors have not been studied in S. mediterranea and are likely to be important candidates in the function of specific cell types. Two types of transcription factors found in 75% of eukaryotic species listed in the comparative genomics resource PhyloPro v2 [53] were not predicted in S. mediterranea: AF-4 (a transcriptional activator that has previously been implicated in childhood lymphoblastic leukemia, mental retardation, and ataxia [56]) and Myc_N (a leucine zipper-type transcription factor implicated in cell cycle progression, cell death, and transformation). The loss of this latter transcription factor in particular suggests that planarians may have adopted an alternate mechanism of regulating Myc's canonical roles in cell proliferation and cell death.

Metabolic reconstruction reveals biochemical pathways distinct from those of parasitic helminths
Genome-scale metabolic reconstructions provide a powerful route to interrogate the metabolic capabilities of an organism [57][58][59][60]. Here we applied an integrated pipeline, developed in house (see Methods), to compare the metabolic potential of S. mediterranea with those derived from other helminths, human, and mouse ( Fig. 3b). Among notable pathways present in S. mediterranea but not in other platyhelminths are several involved in fatty acid metabolism, branched chain amino acid metabolism, mucin-type O-glycan biosynthesis, and one carbon pool by folate. The loss of pathways involved in fatty acid metabolism in the parasitic flatworms may reflect their largely parasitic lifestyles. For example, schistosomes and cyclophyllidean tapeworms spend much of their life cycle in glucose-rich environments (blood and small intestine, respectively) and may therefore have adapted their metabolism to optimize glucose and glycogen as main sources of energy rather than lipids [61,62]. Focusing on amino acid pathways, S. mediterranea displays similar auxotrophies as other helminths; however, a notable exception is branched chain amino acid degradation, which is largely absent from other platyhelminths yet appears to function in S. mediterranea. Conservation of this pathway was surprising given its role in longevity in C. elegans, because S. mediterranea exhibits no evidence of aging and is believed to be immortal [63]. Beyond core metabolic processes, S. mediterranea appears unique among platyhelminths in possessing enzymes required for the production of core 1 mucin-type O-glycans. Such production is likely related to the formation of the mucous secretions that coat the planarian, enabling locomotion, predation, innate immunity, and substrate adhesion [64]. Finally, our comparisons report the presence of several enzymes required for folate interconversion which are otherwise absent in parasitic flatworms. These interconversions provide additional routes for the production of various folate intermediates that are used as co-factors in a variety of metabolic processes, such as tetrahydrofolates involved in nucleotide and amino acid biosynthesis [65].

Spatial annotation of S. mediterranea transcripts by whole-animal scRNAseq
In order to place the annotated transcriptome data in the context of different tissues, the functional information of these transcripts was integrated with spatial information derived from single-cell RNA sequencing (scRNAseq) data of dissociated planarians obtained using Drop-seq technology [66]. The scRNAseq data consist of 51,563 transcripts expressed in 2000 cells. Pruning this dataset to only consider transcripts from our definitive set resulted in a set of 25,168 transcripts expressed in 2000 cells. The R package Seurat [67], which uses an unsupervised clustering approach by combining dimensional reduction with graph-based clustering, was used to cluster the data and discover cell types and states. Based on the set of most variable transcripts in the dataset (n = 4586), Seurat clusters 1195 of the 2000 cells into 11 clusters (Fig. 4a). It is noteworthy that clustering based on the larger set of 51,563 transcripts identified as expressed in the cells recapitulated a similar clustering pattern. Clusters were found to correspond to specific tissues based on the expression of previously described tissue-specific genes (Fig. 4b). In this way, clusters representing epithelial, neural, gut, muscle, parapharyngeal, and stem cells (neoblasts) were identified. Four clusters could not be identified based on previously published planarian gene expression data; however, two of these clusters displayed high expression of the cathepsin homolog Smed-CTSL2 and were thus named cathepsin+ a and cathepsin+ b (Fig. 4c). Cluster 11 displayed enriched expression of Smed-egr-5 and is therefore likely an epithelial subtype (discussed further below; see Fig. 5). Cluster 1 was not specifically enriched for any markers and displayed scattered expression of both neoblast and differentiated tissue markers (Fig. 4b). Its central location on the t-distributed stochastic neighbor embedding (t-SNE) plot, linking the neoblast cluster to the various tissue clusters, led us to conclude that Cluster 1 likely represents transient cell states as neoblasts differentiate along different lineages, and this idea is consistent with recently published scRNAseq studies [29,30].
Differential expression analysis and in situ hybridization demonstrate that the cathepsin + a/b clusters represent mesenchymal populations including pigment cells For the 11 clusters identified by Seurat, cluster markers are identified on the basis of average differential expression. This identified a larger set of cluster markers, ranging from 23 for parapharyngeal cells to 627 for neoblasts (available on figshare https://doi.org/10.6084/ m9.figshare.6852896) [68]. In order to identify the most distinguishing markers, the set of highly differentially expressed genes in a cluster with respect to all other clusters was identified using pairwise assessments of differential expression using a Bayesian approach to single-cell differential expression analysis (SCDE) [69]. This approach builds probabilistic error models for individual cells, capturing both over-dispersion (greater variability than expected) as well as high magnitude outliers and dropout events, thereby providing a more robust approach for detecting differential expression signatures. The clean-up step in this approach is far more stringent than in Seurat, retaining only~60% of the cells compared to the Seurat pipeline (n = 712). For the 11 clusters identified by Seurat and 11,538 transcripts expressed in the cells, transcripts significantly differentially expressed (q value < 0.05) in 10 out of 11 clusters are considered putative markers for the cluster (available on figshare https://doi.org/10.6084/m9.figshare.6852896) [68]. Although there is a larger set of markers detected using Seurat, SCDE also identified unique markers (available on figshare https://doi.org/10.6084/m9.figshare.6852896) [68].
Differential expression analysis identified a significant enrichment for a cathepsin L homolog, Smed-CTSL2 (Sme-dASXL_018694), in the cathepsin + clusters. Cathepsin L is a lysosomal cysteine proteinase with roles in antigen processing and presentation in humans (http://www.uniprot .org/uniprot/P07711). Smed-CTSL2 is expressed across the entire length of the animal in a pattern of branched cells surrounding the gut (Fig. 4c). Interestingly, re-clustering only the cells in the cathepsin + clusters resulted in four distinct subclusters, each with a set of putative markers identified by Seurat (Fig. 4d, Additional file 2: Figure S4A). In situ hybridization of these putative markers demonstrated their unique expression patterns: Subcluster 1 was expressed throughout the mesenchyme (although these cells did not express piwi-1 by scRNAseq) and tightly surrounded the gut (Fig. 4e, Additional file 2: Figure S4B); Subcluster 2 had a punctate expression pattern throughout the animal with randomly localized cell aggregates (Fig. 4e, Additional file 2: Figure S4B); Subcluster 3 was expressed largely within the gut (Fig. 4e); and the final subcluster, interestingly, represented previously described planarian pigment cells based on the enriched expression of published pigment lineage  c Bright field images of Smed-egr-5(RNAi) animals exhibiting reduced eye pigmentation. d In situ hybridization of tyr-1 and immunohistochemistry for ARR in Smed-egr-5(RNAi) regenerating and intact animals. Smed-egr-5 knockdown animals regenerated significantly fewer tyr-1 + cells, while intact animals were comparable to controls. Images are 20-μm confocal z-stacks through the left eye. Scale bars = 10 μm. e The majority of Smed-egr-5(RNAi) animals displayed normal photoreceptor neuron regeneration and reinnervation by ARR staining following head amputation (70%), although some abnormalities were also observed. Images are 40-50 μm confocal z-stacks. Scale bars = 50 μm. f Smed-egr-5(RNAi) and control(RNAi) animals displayed comparable eye regeneration following right eye resection. Images are 20-μm confocal z-stacks. Scale bars = 50 μm. g At a lower dose of RNAi (3 feeds), Smed-egr-5(RNAi) animals exhibited some optic cup regeneration at 7 dpa (red arrowheads), which were subsequently lost at later time points. n = 10-20 for all experiments markers, such as pbgd-1 (Fig. 4e) [54]. Importantly, markers for each of these subclusters were found to be co-expressed to varying degrees in Smed-CTSL2 + cells by double fluorescent in situ hybridization (FISH), consistent with the scRNAseq data (Additional file 2: Figure S4C-F). Interestingly, Subcluster 3 cells also expressed the neoblast marker piwi-1 by scRNAseq (Additional file 2: Figure S4G). As an actively cycling population, the neoblast population is lost following a lethal dose of 6000 rads of irradiation. Likewise, the mesenchymal component of ctcfl (the Subcluster 3 marker) expression was found to be irradiation-sensitive, consistent with its partial expression in neoblasts (Additional file 2: Figure S4H).

Transcription factor analysis reveals cell type-specific expression
Mapping the 843 transcription factors to each cluster identified 30 exhibiting differential expression in specific clusters (significantly upregulated in 8/10 pairwise comparisons) (Additional file 4). Clusters that correspond to muscle, epithelial, and parapharyngeal cell types were associated with the most (7, 7, and 3, respectively) cluster-specific transcription factors, reflecting their generally higher number of differentially expressed transcripts (Additional file 4). Although neoblasts expressed a high number of transcription factors (n = 8), only 1 was cluster-specific. As expected, the most enriched transcription factor domains (zf-C2H2 and LIM) were also the most enriched in the cluster-specific transcripts. However, it is interesting to note that the Ets domain was associated with cluster-specific transcription factors in both epithelial progenitors and Cluster 11, with similar patterns of expression observed in epithelial progenitors and Cluster 11.
Aside from cluster-specific transcription factors, we identified five transcription factors that were abundant and ubiquitously expressed in all clusters (Additional file 4), comprising a Linker_histone domain involved in nucleosome assembly (SmedASXL_006919), and four CSDs, which are present in DNA-and RNA-binding proteins, and implicated in transcriptional regulation.
Analysis of differentially expressed transcription factors identifies the Cluster 11-specific Smed-egr-5 as a regulator of optic cup regeneration Expression of Smed-egr-5 was specific to the unidentified Cluster 11 (Fig. 5a, b). Previous work on Smed-egr-5 demonstrated a striking homeostatic phenotype in which worms exhibited tissue regression and ultimately lysed [70]. Consistent with previous reports, we observed Smed-egr-5 expression subepidermally across the animal with enriched expression on the dorsal side (Additional file 2: Figure S5A) and knockdown of Smed-egr-5 with a high dose of double-stranded RNA (dsRNA) RNAi food (2× dose) resulted in the previously described phenotype (Additional file 2: Figure S5B). dFISH revealed a very low degree of co-localization between Smed-egr-5 and the early epithelial progenitor marker prog-2, but nearly 95% of Smed-egr-5 + cells co-expressed the late epithelial progenitor marker AGAT-1 (Additional file 2: Figure  S5C). Because of the cluster specificity of Smed-egr-5, we sought to further characterize its function by using a lower dose of dsRNA (1× dose) to attempt to uncover further phenotypes. With our 1× RNAi food, we did not observe major defects in epithelial regeneration in Smed-egr-5 knockdown animals (Additional file 2: Figure  S5D); rather, we uncovered a new role for Smed-egr-5 in eye regeneration. After eight feeds of 1× RNAi food, the new head tissue in Smed-egr-5(RNAi) regenerating animals appeared to lack eyes (Fig. 5c). To determine the extent of the missing eye tissue, Smed-egr-5(RNAi) animals were amputated 3 days after the eighth RNAi feed (8fd3) and were allowed to regenerate for 10 days. Regenerating animals were then stained for the optic cup marker Smed-tyrosinase-1 (tyr-1) as well as anti-ARRESTIN (ARR), which marks the optic cup, photoreceptor neurons, and optic nerves. Smed-egr-5(R-NAi) animals regenerated significantly fewer tyr-1 + optic cup cells (p < 0.05), and the cells that did regenerate had noticeably weaker tyr-1 expression (Fig. 5d). There were no apparent eye defects in homeostatic animals (Fig. 5d). ARR staining, on the other hand, revealed largely normal regeneration and reinnervation of photoreceptor neurons, although tissue organization was disrupted in a minority of animals (Fig. 5d, e). Because tyr-1 and ARR staining in intact animals appeared largely normal, we hypothesized that Smed-egr-5 is required specifically during optic cup regeneration.
To test this hypothesis further, an eye scratch assay was performed in which the right eye was resected without significant injury to the surrounding tissue. Previous work has demonstrated that this injury is not sufficient to illicit a regenerative response from the neoblasts; alternatively, the missing eye is restored by maintaining homeostatic levels of new cell incorporation and decreasing the rate of cell death [71]. At 14 days following eye resection, Smed-egr-5(RNAi) animals and control(RNAi) animals had comparable levels of eye restoration, supporting the hypothesis that eye homeostasis is independent of Smed-egr-5 (Fig. 5f).
Interestingly, when Smed-egr-5(RNAi) animals were amputated at an earlier time point of 3fd3, optic cup regeneration was observed at 7 days post-amputation (dpa); however, these cells were subsequently lost at later time points post-amputation (Fig. 5g). The time-sensitive nature of this phenotype suggested that Smed-egr-5 may be involved during the earliest stages of optic cup differentiation: optic cup progenitors that are still remaining after three RNAi feeds are capable of differentiating, but at later time points this progenitor population becomes exhausted and optic cup regeneration ultimately fails. From these data we hypothesize that Smed-egr-5 plays a role in the production of optic cup progenitors. Thus, the lack of an observable homeostatic phenotype may simply be a consequence of the slow turnover of optic cup cells, and it remains possible that optic cup homeostasis may fail at later time points post-RNAi. Further studies at the neoblast level will help to elucidate the precise mechanisms by which Smed-egr-5 promotes proper optic cup regeneration.

Systematic analysis of enriched Gene Ontology terms recapitulates cluster cell types
To provide deeper insights into functional properties associated with each cluster, we performed a GO enrichment analysis. GO mappings for 5900 transcripts expressed in the clusters were obtained through sequence similarity searches of putative homologs with GO annotations from model organisms H. sapiens, Mus musculus, C. elegans, Danio rerio, and D. melanogaster. Although these 5900 transcripts capture only~10% of all transcripts identified in the scRNAseq data, statistically enriched terms were found to complement the previous marker gene analysis, with five of ten clusters consistent with previous cluster definitions: muscle, neural1, neural2, neoblast, and epithelial progenitors (Fig. 6a, Additional file 5). For example, the top ten enriched terms for muscle include terms such as structural constituent of muscle, muscle contraction, and muscle thin filament tropomyosin; neoblast is associated with many terms related to chromosomes and DNA replication, reflecting the high turnover associated with these cells; epithelial is enriched in terms related to endoplasmic reticulum, likely reflecting protein secretion associated with mucoid tissue [72]; and neural1 and 2, although displaying fewer enriched terms than the other tissues, are largely associated with neural functions. Our ability to identify similar consistent patterns of annotations in other clusters is probably limited due to the unavailability of specific GO terms for certain cell types (e.g., parapharyngeal) or due to lower numbers of cells (e.g., < 20 for gut cells) and significantly differentially expressed transcripts in these clusters.

Analyzing correlated gene expression across cell populations reveals transcriptional similarities between distinct cell clusters
In order to identify the set of known/novel subpopulations of cells sharing co-expressed sets of transcripts, we applied the Pathway and Geneset Overdispersion Analysis (PAGODA) component of the SCDE package [73]. This method identifies both the set of GO terms (assigned based on 1:1 orthologs of human) as well as de novo transcript sets consisting of well-correlated gene expression profiles. In this method, since multiple GO terms and de novo gene sets may comprise a common set of genes, clusters sharing the same set of genes are combined to arrive at a final set sharing coordinated variability in expression among the measured cells.
Our analysis reveals a set of four non-redundant clusters, two of which are shown in Fig. 6b. Note, while cell labels were not used during PAGODA, hierarchical clustering of the significantly correlated modules largely recapitulated the patterns of cell clustering generated by the Seurat analysis, especially for muscle, epithelial progenitor, and neural cells. Indeed, epithelial progenitor cells display the most distinct pattern of coordination, which PAGODA associates with Cluster 11 cells. The hierarchical clustering also places the gut and cathepsin + cells together, suggesting that they share transcriptionally co-regulated transcripts. One of the clusters corresponds to a set of cytoskeletal-related proteins in epithelial progenitor cells, as it is enriched in actins, dyneins, and FERM-domain-containing protein (found in several cytoskeletal-associated proteins [74]). The cluster also consists of several unannotated proteins, suggesting their likely involvement in cytoskeleton-related aspects. Although cytoskeletal-related proteins are found in all eukaryotic cells, they are likely to be enriched in epithelial cell types given the role of the cytoskeleton in epithelial cell polarity and intracellular trafficking [75,76]. Although the second "cathepsin+ specific" cluster consists of proteins annotated to be involved in the lipid metabolic process in the lysosome [77,78], phosphorylation/dephosphorylation [79], and cytoskeletal processes, it is unclear as to why these transcripts are co-expressed, opening up novel avenues for experimental interrogation. Reassuringly, Smed-CTSL2 and SmedASXL_009754 (encoding the cathepsin domain) are also identified in this cluster, emphasizing its abundant and unique expression in these cells.

scRNAseq data reveal tissue-specific patterns of metabolic pathway expression
The availability of cell-specific expression profiles generated through scRNAseq raises the intriguing possibility of identifying tissue-specific expression patterns for metabolic enzymes. Applying the hypergeometric test to mean enzyme expression (calculated using SCDE) for each cluster allowed the identification of significantly upregulated or downregulated metabolic pathways, as defined by the Kyoto Encyclopedia of Genes and Genomes (KEGG) [80] (Additional file 6). Consistent with expectations, neoblasts were identified as the most metabolically active cell type followed by muscle and epithelial progenitors (Fig. 7a). The most significantly upregulated pathways are glycolysis/gluconeogenesis in  Fig. 6 Co-expressed sets. a Heatmap depicting the top 20 GO terms significantly enriched in each cluster along with the average expression of transcripts per GO term. The total numbers of statistically significant comparisons and upregulated transcripts for each GO term are also indicated alongside. b Unlabeled hierarchical clustering of cells based on GO gene sets and de novo gene sets consisting of significantly co-expressed offsets of transcripts with very similar gene expression profiles, generated using PAGODA. Two of the most significantly co-expressed modules are indicated, along with the changes in their expression A B Fig. 7 Differential expression of metabolic pathways in clusters. a Distribution of significantly upregulated and downregulated enzymes in each cluster based on pairwise comparisons of log 2 fold change in expression between clusters. b Schematic of differential expression in purine metabolism in neoblast and neural cell types muscle (13/21 enzymes upregulated), supporting an increased need for energy production, and purine metabolism in neoblast (25/35 enzymes upregulated) and neural2 (9/35 enzymes upregulated) cell types (Additional file 6). The purine metabolites adenine and guanine can be synthesized in two distinct pathways: the de novo pathway from CO 2 , glycine, glutamine, aspartate, N 10 -formyltetrahydrofolate and ribose-5-phosphate, starting with phosphoribosyl pyrophosphate (PRPP) and ending in inosine monophosphate (IMP) synthesis; and the salvage pathway, which recycles purine bases by degradation of nucleic acids and nucleotides (Fig. 7b). The purine nucleotides adenosine monophosphate (AMP), guanosine monophosphate (GMP), and xanthosine monophosphate (XMP) are synthesized from IMP. The corresponding trinucleotides lead to generation of intracellular secondary messengers, such as cyclic AMP (cAMP) and cyclic GMP (cGMP). Conversely, the purine nucleotide monophosphates can also be generated by the salvage pathway, by attaching free purine bases to PRPP: via the hypoxanthine-guanine phosphoribosyltransferase (HGPRT) enzyme for IMP, XMP, and GMP synthesis and adenine phosphoribosyltransferase (APRT) for AMP synthesis. As expected, several enzymes of the de novo pathway are upregulated in neoblasts, along with HGPRT of the salvage pathway; however, synthesis of secondary messengers is downregulated. In contrast, there is a significant upregulation of enzymes producing cAMP and cGMP in cells of the neural2 cluster. It is worth noting that neoblasts, in addition to upregulated purine metabolism, are also enriched for pyrimidine metabolism (21/24 enzymes) and one carbon pool by folate (10/11 enzymes upregulated). The enriched synthesis of folate derivatives likely provides the carbon units powering the de novo synthesis of purines and pyrimidines.

Discussion
In this study, starting with an initial set of 83,469 transcripts, we used a hierarchical tiered approach based on protein prediction algorithms of varying stringency and genome assembly mapping to define a high-confidence set of 35,232 transcripts, with 33,487 transcripts (~95% of transcriptome) mapping to 20,483 loci associated with the recently published dd_Smes_g4 S. mediterranea genome [10]. The number of mapped loci is consistent with the number of gene models supported by RNA sequencing (RNAseq) data (n = 19,794) for the closely related regeneration-competent flatworm Macrostomum lignano [81], supporting the quality of the filtered transcriptome. The usage of a tiered approach, which differs from that used to generate other integrated transcriptomes, i.e., PlanMine [19], Oxford [14], and SmedGD [9], reveals that there are 5% unique transcripts in the Toronto transcriptomeof which 20% are supported by homology mapping and 74% by genome assembly mapping, adding to the existing S. mediterranea repertoire. Further, assessment of transcriptome completeness in terms of core eukaryotic and metazoan gene sets as defined by BUSCO v1 [38] reveals that, although the Toronto and PlanMine transcriptomes have the greatest coverage (81% of "core" eukaryotic genes, 78% of "core" metazoan genes), the Toronto dataset also comprises the fewest duplicates in comparison. However, we note that this could also be an artifact of transcript length, potential fusion products from mis-assembly, or spliceoforms, which we did not assess and may be superior in other datasets.
A systematic and comparative bioinformatics analysis of the Toronto transcriptome with the genomes of human, mouse, C. elegans, and close platyhelminth relatives reveals an abundance of transposase-related domains (270 transcripts; DNA transposons of type DDE_1 and DDE_Tnp_1_7), MATH domains (99 transcripts; matrix metalloproteases and TNF-receptor associated factors) and cadherins (100 transcripts) in the planarian. Although the presence of transposable elements is corroborated by previous studies in S. mediterranea [82][83][84] and the basal flatworm M. lignano [81], it is important to note that they are expressed at low RPKM and only a small percentage appear active. Of the 99 transcripts with MATH domains, most are likely to be homologs of TRAF proteins, involved in signal transduction, on the basis of their top homologs. In light of studies supporting the role of homologs of human TRAF-3 and TRAF-6 proteins in immune response in the closely related planarian Dugesia japonica [85], the repertoire of putative TRAF proteins identified in this study provides candidate transcripts that can be tested for their role in planarian immunity. Cadherins are involved in regulating cell-cell adhesion, morphogenesis, and cell recognition [51,52], with additional roles in cellular positioning and maintenance during and after development [86]. Phylogenetic analysis of putative cadherins obtained from humans, S. mediterranea, and other helminths predicts SmedASXL_013539 to be a calsyntenin-like protein, an ortholog of CASY-1 in C. elegans, which has been shown to be essential for learning [87], and CLSTN-1 and CLSTN-2 in humans, implicated in axonal anterograde transport and modulation of post-synaptic signals [88]. Functional characterization of these genes by RNAi may provide novel insights regarding immunity and learning, respectively, in planarians.
Our current understanding of S. mediterranea metabolism is limited [89]. Here we used an established enzyme prediction pipeline [90] to perform a metabolic reconstruction for S. mediterranea. Comparative analyses with other flatworms reveal that S. mediterranea encodes pathways for alternate sources of energy production, such as fatty acid metabolism and branched chain amino acid degradation. Our analyses also identified enzymes responsible for core 1 mucin-type O-glycosylation (notably absent in parasitic flatworms), which may be involved in the formation of the mucous coating, which is involved in locomotion, predation, innate immunity, and substrate adhesion [64].
Several studies have analyzed the role of transcription factors in S. mediterraneainvolving pigmentation [91], gametogenesis [92], epidermal lineage differentiation [93], regeneration [94], and glial cells [95]. Interestingly, Scimone et al. combined RNA sequencing of neoblasts from wounded planarians with expression screening to identify 33 transcription factors and proposed that cell fate for almost all cell types is decided by expression of distinct transcription factors in the neoblast cells [55]. In this study, we used a combination of profile-based approaches to predict 841 putative transcription factors in S. mediterranea. A comparative analysis of putative transcription factors with other eukaryotic species reveals that transcription factor classes belonging to zf-C2H2, Homeobox, zf-BED, bZIP, and HLH are well represented in most species. Several others, such as CSD, Ets, and bZIP-map, well represented in S. mediterranea and vertebrates, have not been studied in the planarian. Studying these transcription factors in S. mediterranea might provide insights into the understanding of the regeneration process.
Several whole-organism as well as tissue-specific bulk RNAseq analyses investigating gene expression differences between two or more treatment conditions have been undertaken in S. mediterranea. To date, 32 RNA-seq/transcriptome datasets are currently available through the NCBI Gene Expression Omnibus (GEO). These experiments provide insights into factors required for restricting injury responses in planarians [96], signaling in planarian glia [95], tissue embryogenesis, homeostasis, and regeneration [97], and transcriptional changes in neoblasts [98]. However, recent developments in scRNAseq technology [99] have provided a novel approach to more directly assess functional differences between different cell populations [100,101]. Recently, scRNAseq has been adopted by studies in S. mediterranea. A comprehensive study by Wurtzel et al. [26] https:// doi.org/10.1016/j.devcel.2015.11.004 using smart-seq2 scRNAseq technology on 619 cells predicted 13 distinct cell clusters and defined 1214 unique tissue markers. This landmark study showed that a generic wound response transcriptional program is activated in almost all cells irrespective of the injury, with most wound-induced genes expressed in muscle, epidermis, and stem cells [26]. A comparison of the cluster markers in our study with those from Wurtzel et al. [102] shows that, although the majority of the cluster markers are shared for muscle (109/122), neural (67/74), and neoblast (87/94) cells, several unique cluster markers are found from this study. Further, Cluster 11 shares 105/133 cluster markers with epithelial cell types, consistent with the presence of AGAT-1 + Smed-egr-5 + cells in this cluster (Additional file 7).
In this study, to better understand the dynamics of the transcriptome in a spatial context, we applied scRNAseq to~2000 cells, from which 25,168 transcripts were identified as expressed in at least one cell. Cluster analysis revealed 11 major clusters, with marker mapping identifying them to be associated with muscle, neural, neoblast, epithelial, and gut tissues, as well as a large cluster of cells likely representing transient transition states during neoblast differentiation (Cluster 1). Further, three novel clusters were identified: two cathepsin+ clusters consisting of four distinct mesenchymal cell types and a Smed-egr-5 + cluster involved in optic cup regeneration. Reassuringly, the cell types of four clustersmuscle, neural, neoblast, and epithelial cellswere recapitulated on the basis of GO term assignments from 1:1 orthologs of model organisms for the most differentially enriched transcripts in these clusters, demonstrating the ability to identify cell types solely on the basis of enrichment of GO terms if GO term assignments are available for differentially enriched transcripts. Differential expression analysis of transcription factors in these clusters identified several cluster-specific factors likely associated with driving the morphogenesis and maintenance of tissue-specific biochemical processes. Analyzing the differential expression of metabolic pathways in these clusters identified neoblast cells as the most metabolically active cell type in S. mediterranea, with highly upregulated purine and pyrimidine metabolism and folate interconversions for providing the key metabolic precursors for nucleotide production. Analysis of purine metabolism with respect to different cell types revealed additional cell-specific patterns of expression, including the upregulation of both de novo and salvage biosynthetic pathways in neoblast cells, as well as the upregulation of intracellular secondary messengers involved in neuronal signaling. Furthermore, our study revealed four cadherin and two MATH domain proteins to be significantly upregulated in neoblast cells, whereas one cadherin and four MATH domain proteins are significantly upregulated in neural cells, providing testable hypotheses for learning more about immunity and learning in planarians.
It should be noted that during the revision of this manuscript, two new studies describing single-cell sequencing in S. mediterranea were published [29,30]. Reassuringly, despite these new studies generating sequence data from 22,000 and~67,000 cells respectively, the results presented in both papers are consistent with our own findings. For example, our finding that pigment cells form a subcluster within the larger cathepsin+ cluster is consistent with the subclustering analysis performed in the Fincher study, in which pbgd-1 was found to mark a specific cathepsin+ subcluster [29]. Further, saposinB-2, which we found to be a specific marker for the cathepsin+ subcluster 2, is expressed in a cathepsin+ subcluster from the same study. This suggests that smaller scale datasets, such as the one presented here, are sufficient to recapitulate many of the conclusions of larger-scale studies and consequently represent a valuable experimental template to assay specific RNAi phenotypes with single-cell sequencing in the future.

Conclusions
Here we present a definitive set of transcripts for the freshwater planarian Schmidtea mediterranea. We further annotate all genes with identifiable homology and identify gene family expansions and losses. Interestingly, TRAF proteins have been disproportionately increased, while Myc and AF-4 transcription factors are absent. A genome-scale metabolic reconstruction was then performed to identify metabolic pathways conserved in platyhelminths, those that have been lost in parasitic flatworms and those that represent lineage-specific innovations in S. mediterranea. Sequencing transcripts associated with 2000 individual cells identified cell types by differential gene expression and further revealed additional genes and pathways specific to each cell type. These analyses also uncovered a novel cell type associated with a novel mesenchymal cell population. In summary, these analyses build a foundation of cell types and gene conservation profiles that will inform future gene function studies.

Methods
Culturing of S. mediterranea, in situ hybridization, and RNA interference Asexual individuals of S. mediterranea CIW4 strain were reared as previously described [103]. In situ hybridization was performed as previously described [18,104]. RNAi was performed as previously described [54], with either three or eight feeds as indicated in the text.

Generating a high-confidence S. mediterranea transcriptome
The initial transcriptome of 83,469 transcripts was an assembly collated from five separate experiments and more than 1 billion RNA-seq reads from whole animals, purified tissues, RNAi conditions, and irradiated whole animals [18,[31][32][33] (NCBI Bioproject PRJNA215411). The resulting transcriptome was filtered using various criteria in order to arrive at a high-confidence set of putative protein-coding transcripts (Fig. 1a). As a first step, likely contaminants were identified by a BLASTn (from BLAST+ 2.2.28) [40] search against the protein nucleotide (nt) database (2016) [105] to remove sequences matching other species at a sequence identity and query coverage cutoff of 95% (n = 237) as well as those matching vector sequences (n = 8). Next, likely mis-assembled transcripts were removed by identifying all transcripts with ≥ 25 unmapped bases to the transcriptome (n = 2387). Clustering approaches did not reduce the initial transcriptome to the expected range observed in regeneration-competent species such as M. lignano and D. japonica, suggesting the presence of contaminants, misassembled transcripts, split transcripts, alternative splice variants, and/or leaky transcripts. Therefore, the initial transcriptome was scrutinized via a multi-layered approach to identify potential protein-coding transcripts. The transcriptome was parsed through the prot4EST v3.1b [106] pipeline, an integrated approach which overcomes deficits in training data in order to convert transcripts into proteins. This multi-tiered program identifies coding transcripts in various stages. The first step identifies homologs of known RNA and protein sequences using the BLAST suite [40] -BLASTn (from BLAST 2.2.28) against the SILVA database (release 115) [107] at an E-value of 1e-65 for identifying RNA transcripts, BLASTx against the Mito-Miner database (v3.1) [35] at an E-value of 1e-08 and against the UniProt database [34] at an e value of 1e-05 for identifying mitochondrial and nuclear transcripts, respectively. From the remaining transcripts, the second step identifies likely protein-coding transcripts using ESTscan [v3.0.3] [39], a hidden Markov model (HMM)-based model trained to be error-tolerant, using a simulated S. mediterranea training set. Finally, the remaining transcripts are processed to identify the longest string of amino acids uninterrupted by stop codons from a six-frame translation of the sequence (LongestORFs). From the set of categorized transcripts, all transcripts with query coverage spanning two thirds of the reference sequence in RNA/mitochondrial/nuclear databases are retained. The rest of the transcripts are retained only if there is any support in terms of the following: (1) homology with respect to conserved eukaryotic gene sets (CEGMA v2.5 [37] and BUSCO v1.1 [38] using BLASTx at an E-value of 1e-08) and other helminth transcriptome EST datasets obtained from the NCBI (B. glabrata, C. sinensis, C. gigas, D. japonica, D. ryukyuensis, E. granulosus, E. multilocularis, H. robusta, H. medicinalis, H. microstoma, M. lignano, M. californianus, O. viverrini, S. japonicum, S. mansoni, T. solium) using BLASTn at an E-value of 1e-15; (2) annotation by InterPro [36] at an E-value of 1e-03; and (3) co-location of the draft S. mediterranea genome with ESTs from NCBI, transcripts from the Oxford dataset (v0.1) [14], or transcripts from SmedGD v2.0 using Spaln v2 [41] at a stringency filtering of F2 (corresponding to alignment length > 200 bp, sequence identity ≥ 93%, query coverage ≥ 93%).

Comparison with PlanMine genome and transcriptome
The Toronto transcriptome was mapped onto the Plan-Mine genome [10] using Spaln v2 [41] at stringency filtering cutoffs corresponding to F2 (sequence identity ≥ 93%, query coverage ≥ 93%) and F1 (sequence identity ≥ 75%, query coverage ≥ 75%) in order to identify the extent of overlap. Subsequently, the transcriptomes were compared using BLASTn [40] searches against each other using a relaxed word size (n = 7) in order to improve the stringency of the searches. BLASTn matches of the Toronto transcriptome to the PlanMine transcriptome were pruned based on the nearest bit score cutoff corresponding to the number of overlapping matches to the genome identified at F1 cutoff (corresponding to a bit score value ≥ 40). Based on this cutoff, matches were identified between the Toronto, PlanMine, Oxford, and SmedGD transcriptomes.

Functional annotation of the transcriptome
The predicted protein sequences generated from the high-confidence transcriptome were functionally annotated by (1) HMM searches against the curated Pfam-A database v31 using the PfamScan tool with hmmer-3.1b1 [44] at default cutoffs. Only those matches with an E-value cutoff of < 0.001 were considered for further analysis;

RPKM calculation
The expression levels of the transcripts were calculated by mapping the reads from 58 RNA-seq results (listed as the column headers under the RPKM section in Additional file 1) onto the initial transcriptome using Burrows-Wheeler Aligner (BWA) [110] and obtaining the number of reads mapped for each transcript. The normalized expression levels were quantified in RPKM units for each transcript for each RNA-seq experiment using the formula: RPKM = Number of Reads/(Transcript Length/1000 * Total Num Reads/1,000,000) where Total Num Reads consisted only of those transcripts with ≥ 10 reads mapped to them in a sample. Next, the mean, standard deviation, and median RPKM values for each transcript were calculated based on the number of RNA-seq experiments where the transcript was expressed. The mean values of all transcripts in the definitive transcriptome were used to derive a percentile distribution of RPKM values, which is used as a guide to derive the average level of expression of a transcript (low < 50th percentile, high > 20th percentile, medium ≤ 20th percentile and ≥ 50th percentile).

Enzyme annotation of the predicted proteome
For each of the predicted protein sequences, an initial set of enzyme commission (EC) predictions was obtained from several methods: (1) density estimation tool for enzyme classification (DETECT) v1.0 run using default parameters (here we retained hits with Integrated Likelihood Score (ILS) cutoff ≥ 0.9 from the top predictions file which also had ≥ 5 positive hits) [57]; (2) BLASTP (from BLAST+ 2.2.28) run against the Swiss-Prot database (release 2014-08) at an E-value cutoff of 1e-10; the enzyme annotations of top hits in the Swiss-Prot database were mapped to the query sequence [40]; and (3) PRIAM enzyme rel. Feb-2014 run using relaxed cutoffs specified for genome-wide annotations of organisms (minimum probability > 0.5, profile coverage > 70%, check catalytic -TRUE) [58]. From these assignments, a set of consolidated high-confidence predictions was derived using in-house scripts by retaining only those predictions identified by both PRIAM and BLASTP and combining them with the predictions from DETECT. Percent pathway conservation was calculated for the set of metabolic pathways as defined by KEGG v70 [80] using the following formula: (Number of predicted ECs in a KEGG pathway × 100)/Total number of ECs in the KEGG pathway.

Transcription factor prediction
The InterProScan v5.15.54.0 [108] outputs for all 35,235 high-confidence predicted protein sequences were scanned as follows in order to identify a set of putative transcription factors: (1) InterProScan hits with the description "transcription factor", (2) InterProScan hits to the Pfam families listed in the curated transcription factor database DNA-binding domain (DBD) v2.0 [117], (3) InterProScan hits to the Superfamily families listed in DBD v2.0. The hits from all of the above criteria were consolidated to arrive at the final predicted set of transcription factors for the organism.

Transposon analysis
RepeatMasker (2013) was used to predict repeats for the SmedAsxl genome v1.1. All transcripts assigned DDE transposase domains were mapped onto the masked SmedAsxl genome with the F2 cutoff of Spaln v2 [41] and searched for the presence of repetitive elements. For repetitive elements found within the mapped region, sequence regions flanking 1000 bp on either side of the repetitive element were extracted and its sequence divergence with the consensus of the repeat element calculated using the Needleman-Wunsch algorithm from the European Molecular Biology Open Software Suite (EMBOSS) package. A histogram of the extent of sequence divergence was analyzed in order to identify likely active elements, characterized by sequence divergence ≤5% from consensus element [118].

Generation of single-cell RNA-seq data
For single-cell RNA sequencing, a whole-animal cell suspension (in calcium-magnesium-free (CMF) + 10% glucose solution) was stained with the cell viability dye calcein (0.2 μg/ml), and calcein-positive cells were collected by fluorescence-activated cell sorting (FACS). Cells were then processed through a Drop-seq instrument and complementary DNA (cDNA) libraries were prepared as described in [66]. Libraries were sequenced on an Illumina NextSeq500 to a total depth of~480 million reads. The data are available at the NCBI GEO database under accession number GSE115280 (https:// www.ncbi.nlm.nih.gov/gds/?term=GSE115280) [119]. Reads were aligned to the S. mediterranea SmedASXL transcriptome assembly under NCBI BioProject PRJNA215411 using Bowtie2 with 15-bp 3′ trimming.

Identification of clusters and cluster markers using Seurat
To identify cell clusters enriched for transcriptionally co-expressed profiles, single-cell RNA-seq data were processed against the definitive Toronto transcriptome using the Seurat [67] pipeline while considering the standard default quality cutoffs optimized for a dataset of size~3000 cells, i.e., min.genes = 200, min.cells = 3, tot.expr = 1e4. The resolution parameter in the FindClusters function was varied from 0.4 to 4, and a resolution of 1 was chosen as it yielded the most visually distinct clustering pattern. In Seurat [67], cluster markers were identified using the FindAllMarkers function of the Seurat pipeline by considering transcripts that are expressed in at least 25% of the cells in the cluster, with an average expression ≥ 25% in comparison to their expression in all other clusters. The significance of the differential expression is calculated using the "bimod" likelihood-ratio test for single-cell gene expression [120] for all cells in one cluster vs all other cells and expressed as p values.

Differential expression of transcripts and identification of cluster markers in SCDE
Differential expression of transcripts between clusters was calculated using the SCDE R package, which employs a Bayesian approach to single-cell differential expression analysis [69], considering only those cells with a minimum library size of 500, and only those transcripts mapping to ≥ 10 reads and detected in ≥ 5 cells, since this yielded at least ten cells per cluster. Differential expression was calculated for all-vs-all pairwise combinations of clusters classified using Seurat, and the log 2 fold change and p values were noted. All transcripts that are significantly upregulated in 9/10 pairwise comparisons are considered as cluster markers.

Hypergeometric test for KEGG metabolic pathways
The enrichment of differentially expressed transcripts (both upregulated, corresponding to a log 2 fold change > 1; and downregulated, corresponding to a log 2 fold change < − 1, according to SCDE) was assessed using a hypergeometric test (using the phyper function in R) for all pairwise combinations of clusters classified using Seurat. All KEGG pathways with a p value < 0.05 were considered to be enriched.

Hypergeometric test for analyzing enrichment of Gene Ontology terms
Gene Ontology (GO) refers to a database providing a structured vocabulary for annotating genes [43]. The genes are annotated using specific biologically relevant terms corresponding to three main categories: Biological Process (BP), Molecular Function (MF), and Cellular Compartment (CC). Schmidtea transcripts were annotated with the GO terms from 1:1 orthologs from five model organisms: H. sapiens, M. musculus, D. rerio, C. elegans, and D. melanogaster, as identified by Inparanoid (annotations downloaded from GO website http://gene ontology.org/page/download-annotations). The annotations were transferred for GO terms designated by all methods other than Inference by Electronic Annotation (non-IEA) on the basis of Inparanoid mapping, using in-house scripts. The enrichment of significantly upregulated transcripts associated with the GO term (log 2 fold change > 1 calculated using SCDE) was assessed using a hypergeometric test (using the phyper function in R) for all pairwise combinations of clusters classified using