Mapping single-cell atlases throughout Metazoa unravels cell type evolution

Comparing single-cell transcriptomic atlases from diverse organisms can elucidate the origins of cellular diversity and assist the annotation of new cell atlases. Yet, comparison between distant relatives is hindered by complex gene histories and diversifications in expression programs. Previously, we introduced the self-assembling manifold (SAM) algorithm to robustly reconstruct manifolds from single-cell data (Tarashansky et al., 2019). Here, we build on SAM to map cell atlas manifolds across species. This new method, SAMap, identifies homologous cell types with shared expression programs across distant species within phyla, even in complex examples where homologous tissues emerge from distinct germ layers. SAMap also finds many genes with more similar expression to their paralogs than their orthologs, suggesting paralog substitution may be more common in evolution than previously appreciated. Lastly, comparing species across animal phyla, spanning sponge to mouse, reveals ancient contractile and stem cell families, which may have arisen early in animal evolution.


Introduction
There is much ongoing success in producing single-cell transcriptomic atlases to investigate the cell type diversity within individual organisms (Regev et al., 2017). With the growing diversity of cell atlases across the tree of life Cao et al., 2019;Fincher et al., 2018;Hu et al., 2020;Musser et al., 2019;Plass et al., 2018;Siebert et al., 2019;Wagner et al., 2018), a new frontier is emerging: the use of cross-species cell type comparisons to unravel the origins of cellular diversity and uncover species-specific cellular innovations (Arendt et al., 2019;Shafer, 2019). Further, these comparisons promise to accelerate cell type annotation and discovery by transferring knowledge from well-studied model organisms to under-characterized animals.
However, recent comparative single-cell analyses are mostly limited to species within the same phylum (Baron et al., 2016;Geirsdottir et al., 2019;Sebé-Pedró s et al., 2018;Tosches et al., 2018). Comparisons across longer evolutionary distances and across phyla are challenging for two major reasons. First, gene regulatory programs diversify during evolution, diminishing the similarities in cell-type-specific gene expression patterns. Second, complex gene evolutionary history causes distantly related organisms to share few one-to-one gene orthologs (Nehrt et al., 2011), which are often relied upon for comparative studies Shafer, 2019). This effect is compounded by the growing evidence suggesting that paralogs may be more functionally similar than orthologs across species, due to differential gain (neo-functionalization), loss (non-functionalization), or partitioning (sub-functionalization) events among paralogs (Nehrt et al., 2011;Prince and Pickett, 2002;Stamboulian et al., 2020;Studer and Robinson-Rechavi, 2009).
Here, we present the Self-Assembling Manifold mapping (SAMap) algorithm to enable mapping single-cell transcriptomes between phylogenetically remote species. SAMap relaxes the constraints imposed by sequence orthology, using expression similarity between mapped cells to infer the relative contributions of homologous genes, which in turn refines the cell type mapping. In addition, SAMap uses a graph-based data integration technique to identify reciprocally connected cell types across species with greater robustness than previous single-cell data integration methods (Haghverdi et al., 2018;Hie et al., 2019;Polański et al., 2019;Stuart et al., 2019).
Using SAMap, we compared seven whole-body cell atlases from species spanning animal phylogeny, which have divergent transcriptomes and complex molecular homologies ( Figure 1A-B and Supplementary file 1). We began with well-characterized cell types in developing frog and zebrafish embryos. We found broad concordance between transcriptomic signatures and ontogenetic relationships, which validated our mapping results, yet also detected striking examples of homologous cell types emerging from different germ layers. We next extended the comparison to animals from the same phylum but with highly divergent body plans, using a planarian flatworm and a parasitic blood fluke, and found one-to-one homologies even between cell subtypes. Comparing all seven species from sponge to mouse, we identified densely interconnected cell type families broadly shared across animals, including contractile and stem cells, along with their respective gene expression programs. Lastly, we noticed that homologous cell types often exhibit differential expression of orthologs and similar expression of paralogs, suggesting that the substitution and swapping of paralogs in cell types may be more common in evolution than previously appreciated. Overall, our study represents an important step toward analyzing the evolutionary origins of specialized cell types and their associated gene expression programs in animals.

Results
The SAMap algorithm SAMap iterates between two modules. The first module constructs a gene-gene bipartite graph with cross-species edges connecting homologous gene pairs, initially weighted by protein sequence similarity ( Figure 1C). In the second module, SAMap uses the gene-gene graph to project the two single-cell transcriptomic datasets into a joint, lower-dimensional manifold representation, from which each cell's mutual cross-species neighbors are linked to stitch the cell atlases together ( Figure 1D). Then, using the joint manifold, the expression correlations between homologous genes are computed and used to reweight the edges in the gene-gene homology graph in order to relax SAMap's initial dependence on sequence similarity. The new homology graph is used as input to the subsequent iteration of SAMap, and the algorithm continues until convergence, defined as when the cross-species mapping does not significantly change between iterations ( Figure 1E).
This algorithm overcomes several challenges inherent to mapping single-cell transcriptomes between distantly related species. First, complex gene evolutionary history often results in many-tomany homologies with convoluted functional relationships Nehrt et al., 2011). SAMap accounts for this by using the full homology graph to project each dataset into both its own and its partner's respective principal component (PC) spaces, constructed by the SAM algorithm, which we previously developed to robustly and sensitively identify cell types (Tarashansky et al., 2019). The resulting within-and cross-species projections are concatenated to form the joint space. For the cross-species projections, we translate each species' features into those of its partner, with the expression for individual genes imputed as the weighted average of their homologs specified in the gene-gene bipartite graph. Iteratively refining the homology graph to only include positively correlated gene pairs prunes the many-to-many homologies to only include genes that are expressed in the same mapped cell types.
Second, frequent gene losses and the acquisitions of new genes result in many cell type gene expression signatures being species-specific, limiting the amount of information that is comparable across species. Restricting the analysis of each dataset to only include genes that are shared across species would result in a decreased ability to resolve cell types and subtypes with many species-specific gene signatures. SAMap solves this problem by constructing the joint space through the concatenation of within-and cross-species projections, thus encoding all genes from both species.
Third, the evolution of expression programs gradually diminishes the similarity between homologous cell types. To account for this effect, SAMap links cell types across species while tolerating their differences. Cells are mapped by calculating each of their k mutual nearest cross-species neighbors   in the combined projection. To establish more robust mutual connectivity, we integrate information from each cell's local, within-species neighborhood ( Figure 1D), overcoming the inherent stochasticity of cross-species correlations. Two cells are thus defined as mutual nearest cross-species neighbors when their respective neighborhoods have mutual connectivity. It is important to note that the magnitude of connections is not directly calculated from their expression similarity, allowing cell types with diverged expression profiles to be tightly linked if they are among each other's closest cross-species neighbors. Lastly, SAMap is robust to technical batch effects between datasets that are collected through different platforms. For instance, we have succeeded in running SAMap on datasets containing hundreds of thousands of cells that were collected with different single-cell platforms, including 10X genomics, Drop-Seq, SmartSeq, and MARS-Seq. SAMap runtimes were typically less than an hour on an average desktop computer for the largest dataset we tested (Figure 1-figure supplement  1). Further, SAMap overcomes potential memory issues when running on large datasets by chunking its computationally intensive operations into smaller blocks, saturating the memory usage with respect to the number of cells (Figure 1-figure supplement 1B).

Homologous cell types emerging from distinct germ layers in frog and zebrafish
We first applied SAMap to the Xenopus and zebrafish atlases, which both encompass embryogenesis until early organogenesis Wagner et al., 2018). Previous analysis linked cell types between these two organisms by matching ontogeny, thereby providing a reference for comparison. SAMap produced a combined manifold with a high degree of cross-species alignment while maintaining high resolution for distinguishing cell types in each species (Figure 2A). We measured the mapping strength between cell types by calculating an alignment score (edge width in Figure 2B and color map in Figure 2C), defined as the average number of mutual nearest cross-species neighbors of each cell relative to the maximum possible number of neighbors.
SAMap revealed broad agreement between transcriptomic similarity and developmental ontogeny, linking 26 out of 27 expected pairs based on previous annotations ( Figure 2B and Supplementary file 2; Briggs et al., 2018). The only exception is the embryonic kidney (pronephric duct/mesenchyme), potentially indicating that their gene expression programs have significantly diverged. In addition, SAMap succeeded in drawing parallels between the development of homologous cell types and matched time points along several cell lineages ( Figure 2C). While the concordance was consistent across cell types, we noticed that the exact progression of developmental timing can vary, suggesting that SAMap can quantify heterochrony with cell type resolution.
SAMap also linked a group of secretory cell types that differ in their developmental origin, some even arising from different germ layers (highlighted edges in Figure 2B). Within ectoderm, frog cement gland cells map to zebrafish muc5ac + secretory epidermal cells, and frog small secretory cells (SSCs) map to zebrafish pvalb8 + mucous cells (Janicke et al., 2010). Across germ layers, SSCs also map weakly to zebrafish endodermal cells, and frog ectodermal hatching gland maps to zebrafish mesodermal hatching gland. These cell types are linked through a large set of genes, including proteins involved in vesicular protein trafficking and several conserved transcription factors (TFs) such as myb, foxa1, xbp1, and klf17 ( Figure 2D), which all have documented functions in controlling the differentiation of secretory cell types (Bennett et al., 2007;Dubaissi et al., 2014;Pan et al., 2014). For example, klf17 is expressed in zebrafish and frog hatching glands, and plays essential roles in regulating gland cell specification in both species (Kurauchi et al., 2010;Suzuki et al., 2019). Together, the conserved cell type specification programs (Erwin and Davidson, 2009) between developmentally distinct secretory cells support the notion that they may be transcriptionally and evolutionarily related despite having different developmental origins .
To benchmark the performance of SAMap, we used eggNOG (Huerta-Cepas et al., 2019) to define one-to-one vertebrate orthologs between zebrafish and frog and fed these gene pairs as input to several broadly used single-cell data integration methods, Seurat (Stuart et al., 2019), LIGER (Welch et al., 2019), Harmony (Korsunsky et al., 2019), Scanorama (Hie et al., 2019), and BBKNN (Polański et al., 2019). We found that they failed to map the two atlases, yielding minimal alignment between them ( Figure 2E and Figure 2-figure supplement 1). We also compared the results when restricting SAMap to using the one-to-one orthologs instead of the full homology graph. Even when removing the many-to-many gene homologies and the iterative refinement of the   homology graph, we identified similar, albeit weaker, cell type mappings. This suggests that, at least for the frog and zebrafish comparison, SAMap's performance is owed in large part to its robust, atlas stitching approach. Finally, to test if SAMap is robust to incomplete cell type atlases, we downsampled the frog and zebrafish data by systematically removing individual cell types. We found that cell types whose homologous partners were removed mapped weakly to closely related cell types, some of which were already present in the original mapping. For example, optic cells from both species were also connected to eye primordium, frog skeletal muscles to zebrafish presomitic mesoderm, and frog hindbrain to zebrafish forebrain/midbrain (Supplementary file 3). We observed several new mapping pairs, but their alignment scores were all barely above the detection threshold of SAMap. Moreover, most of these edges were mapped between cell types with similar developmental origins, with the only exception being the zebrafish neural crest mapped to the frog otic placode in the absence of frog neural crest cells. Examining the genes that support this mapping revealed that both cell types express sox9 and sox10, two TFs previously implicated to form a conserved gene regulatory circuit common to otic/neural crest cells (Betancur et al., 2011). Taken together, these results suggest that SAMap is more sensitive in linking homologous cell types compared to other existing methods, exhibits high robustness when applied to incomplete datasets, and yields mapping results that are well supported by conserved gene expression programs.

Paralog substitutions are prevalent between homologous cell types in frog and zebrafish
The key benefit of using the full homology graph is to enable the systematic identification of gene paralogs that exhibit greater similarity in expression across species than their corresponding orthologs. These events are expected to arise as the result of gene duplications followed by diversification of the resulting in-paralogs (Studer and Robinson-Rechavi, 2009). In an alternative scenario, genetic compensation by transcriptional adaptation, where loss-of-function mutations are balanced by upregulation of related genes with similar sequences (El-Brolosy et al., 2019), could also result in this signature.
In total, SAMap selected 8286 vertebrate orthologs and 7093 paralogs, as enumerated by egg-NOG, for manifold alignment. Paralogs were identified as non-orthologous genes that map to the same eggNOG orthology group ancestral to Vertebrata. Among these, 565 genes have markedly higher expression correlations (correlation difference >0.3) with their paralogs than their orthologs (see Figure 3A for examples), and 209 of them have orthologs that are either completely absent or lowly expressed with no cell-type specificity (Supplementary file 4). We term these events as 'paralog substitutions', as the orthologs may have lost or changed their functional roles at some point and were compensated for by their paralogs. Substituting paralogs were identified in most cells types with some (e.g. dorsal organizer) exhibiting higher rates than others ( Figure 3B), suggesting uneven diversification rates of paralogs across cell types. SAMap also linked an additional 297 homologous pairs previously unannotated by orthology or paralogy, but which exhibit sequence similarity and high expression correlations (>0.5 Pearson correlation). These likely represent unannotated orthologs/paralogs or isofunctional, distantly related homologs (Gabaldó n and Koonin, 2013).
We next asked whether paralog substitution rates depend on the evolutionary time since gene duplication. We categorized paralogs by the taxonomic level of their most recent shared orthology group and found that more recent paralogs substitute orthologs at higher rates than more ancient paralogs ( Figure 3C). This observation is consistent with the expectation that less diverged genes may be more capable of functionally compensating for each other. To rule out the possibility that these paralogs were linked spuriously during the homology refinement steps of SAMap, we repeated the paralog substitution analysis on an aligned manifold constructed using only one-to-one orthologs. We identified 70% of the paralog substitutions and observed similar patterns in evolutionary time and cell type dependencies ( . Failure to detect these substitutions was due to inaccurate imputation of gene expressions across species when restricting the mapping to one-toone orthologs, which resulted in weaker alignment with fewer cross-species edges. Altogether, these results illustrate the potential of SAMap in leveraging single-cell gene expression data for pruning the networks of homologous genes to identify evolutionary substitution of paralogs and, more generally, identify non-orthologous gene pairs that may perform similar functions in the cell types within which they are expressed.

Homologous cell types between two flatworm species with divergent body plans
To test if we can identify homologous cell types in animals with radically different body plans, we mapped the cell atlases of two flatworms, the planarian Schmidtea mediterranea (Fincher et al.,  2018), and the trematode Schistosoma mansoni, which we collected recently . They represent two distant lineages within the same phylum but have remarkably distinct body plans and autecology (Laumer et al., 2015;Littlewood and Waeschenbach, 2015). While planarians live in freshwater and are known for their ability to regenerate (Reddien, 2018), schistosomes live as parasites in humans. The degree to which cell types are conserved between them is unresolved, given the vast phenotypic differences caused by the transition from free-living to parasitic habits (Laumer et al., 2015).
SAMap revealed broad cell type homology between schistosomes and planarians. The schistosome had cells mapped to the planarian stem cells, called neoblasts, as well as most of the differentiated tissues: neural, muscle, intestine, epidermis, parenchymal, protonephridia, and cathepsin + cells, the latter of which consists of cryptic cell types that, until now, have only been found in planarians (Fincher et al., 2018; Figure 4A). These mappings are supported by both known cell-type-   Figure 4. SAMap transfers cell type information from a well-annotated organism (planarian S. mediterranea) to its less-studied cousin (schistosome S. mansoni) and identifies parallel stem cell compartments. (A) UMAP projection of the combined manifolds. Tissue type annotations are adopted from the S. mediterranea atlas (Fincher et al., 2018). The schistosome atlas was collected from juvenile worms, which we found to contain neoblasts with an abundance comparable to that of planarian neoblasts .  We next determined if cell type homologies exist at the subtype level. For this, we compared the stem cells, as planarian neoblasts are known to comprise populations of pluripotent cells and tissuespecific progenitors (Fincher et al., 2018;Zeng et al., 2018). By mapping the schistosome stem cells to a planarian neoblast atlas (Zeng et al., 2018), we found that the schistosome has a population of stem cells, e-cells , that cluster with the planarian's pluripotent neoblasts, both expressing a common set of TFs (e.g. soxp2, unc4, pax6a, gcm1) ( Figure 4C-D). The e-cells are closely associated with juvenile development and lost in adult schistosomes Nanes Sarfati et al., 2021), indicating pluripotent stem cells may be a transient population restricted to their early developmental stages. This is consistent with the fact that, whereas schistosomes can heal wounds, they have limited regenerative ability (Wendt and Collins, 2016). SAMap also linked other schistosome stem cell populations with planarian progenitors, including two populations of schistosome stem cells -denoted as m (Tarashansky et al., 2019) and m' -to planarian muscle progenitors, all of which express myoD, a canonical master regulator of myogenesis (Scimone et al., 2017). These likely represent early and late muscle progenitors, respectively, as m-cells do not yet express differentiated muscle markers such as troponin, whereas m'-cells do (Figure 4-figure supplement 2).

Cell type families spanning the animal tree of life
To compare cell types across broader taxonomic scales, we extended our analysis to include juvenile freshwater sponge (Spongilla lacustris) , adult Hydra (Hydra vulgaris) (Siebert et al., 2019), and mouse (Mus musculus) embryogenesis (Pijuan-Sala et al., 2019) atlases. In total, SAMap linked 1051 cross-species pairs of cell types, defined by the annotations used in each respective study. Of the cell type pairs, 95% are supported by at least 40 enriched gene pairs, and 87% are supported by more than 100 gene pairs, indicating that SAMap does not spuriously connect cell types with limited overlap in transcriptional profiles ( Figure 5-figure supplement 1A).
We next extended the notion of cell type pairs to cell type trios, as mapped cell types gain additional support if they share transitive relationships to other cell types through independent mappings, forming cell type triangles among species. The transitivity of a cell type pair (edge) or a cell type (node) can be quantified as the fraction of triads to which they belong that form triangles ( Figure 5A). The majority (81%) of cell type pairs have non-zero transitivity independent of alignment score and the number of enriched gene pairs ( Figure 5-figure supplements 1-2). Cell type pairs with fewer than 40 enriched gene pairs tend to have lower (<0.4) transitivity ( Figure 5-figure  supplement 1B). The transitivity measure can also be used to identify potentially spurious connections. 16% of mapped cell type pairs have zero edge transitivity but non-zero node transitivity. These cell types are connected to only a single member of an interconnected cell type group (motifs 2 and 3 in Figure 5B). Such links may be of lower confidence as they should connect to other members of the group and are thus excluded from downstream analysis.
Among the interconnected groups of cell types, we identified families of neural cells and contractile cells ( Figure 5C). Both cell type families are highly transitive compared to the overall graph transitivity (bootstrap p-value<1Â10 À5 ), meaning that their constituent cell types have more transitive edges within the group than outside the group ( Figure 5D). In addition, the dense, many-to-many connections within the contractile and neural families are each supported by at least 40 enriched gene pairs ( Figure 5E). Consistent with the nerve net hypothesis suggesting a unified origin of neural cell types (Tosches and Arendt, 2013), the neural family includes vertebrate brain tissues, both bilaterian and cnidarian neurons, cnidarian nematocytes that share the excitatory characteristics of neurons (Weir et al., 2020), and Spongilla choanocytes and apopylar cells, both of which are not considered as neurons but have been shown to express postsynaptic-like scaffolding machinery Wong et al., 2019). The contractile family includes myocytes in bilaterian animals, Hydra myoepithelial cells that are known to have contractile myofibrils (Buzgariu et al., 2015), and sponge pinacocytes and myopeptidocytes, both of which have been implicated to play roles in contractility Sebé-Pedró s et al., 2018). In contrast to the families encompassing all seven species, we also found a fully interconnected group that contains invertebrate multipotent stem cells, including planarian and schistosome neoblasts, Hydra interstitial cells, and sponge archeocytes (Alié et al., 2015). The lack of one-to-one connections across phyla is in keeping with recent hypotheses that ancestral cell types diversified into families of cell types after speciation events Arendt et al., 2019). Our findings thus suggest that these cell type families diversified early in animal evolution.

Transcriptomic signatures of cell type families
The high interconnectedness between cell types across broad taxonomic scales suggests that they should share ancestral transcriptional programs . SAMap identified broad transcriptomic similarity between bilaterian and non-bilaterian contractile cells that extends beyond the core contractile apparatus. It links a total of 23601 gene pairs, connecting 5471 unique genes, which are enriched in at least one contractile cell type pair. Performing functional enrichment analysis on these genes, we found cytoskeleton and signal transduction functions to be enriched (p-value<10 À3 ) based on the KOG functional classifications (Tatusov et al., 2003)    ( Figure 6A). These genes include orthology groups spanning diverse functional roles in contractile cells, including members of the adhesion complex that connects cells, actomyosin networks that drive contractility, and signaling pathways that stimulate contraction ( Figure 6B and Supplementary file 5). This observation suggests that contractile cells near the beginning of animal evolution already possessed the broad assemblage of gene modules associated with different functional aspects of derived muscle cell types in extant animals. We also identified several transcriptional regulators shared among contractile cells ( Figure 6B). Previously known core regulators involved in myocyte specification (Brunet et al., 2016) were enriched only in bilaterian (e.g. myod and tcf4/E12) or vertebrate contractile cells (e.g. mef2). In contrast, we found homologs of Muscle Lim Protein (Csrp) and Forkhead Box Group 1 (Larroux et al., 2008) enriched in contractile cells from all seven species. The Fox proteins included FoxC, which is known to regulate cardiac muscle identity in vertebrates (Brunet et al., 2016) and is contractile-specific in all species except schistosome and Spongilla. Notably, we also identified FoxG orthologs to

Signal transduction (T) Cytoskeleton (Z) Extracellular structures (W) Transcription (K) Inorganic ion metabolism (P) Carbohydrate metabolism (G) Secretion and transport (U) Defense mechanisms (V) Secondary metabolites (Q)
Translation ( be enriched in three of the four invertebrates ( Figure 6-figure supplement 1), suggesting that FoxG may play an underappreciated role in contractile cell specification outside vertebrates. Whether the most conserved regulators are positioned at the center of the cell type specification regulatory complex remains an important question to address in future studies.
For the family of invertebrate multipotent stem cells, we identified 3343 genes that are enriched in at least one cell type pair and observed significant enrichment (p-value<10 À3 ) of genes involved in translational regulation such as RNA processing, translation, and post-translational modification ( Figure 6C). We also observed cell cycle and DNA replication genes, as expected for actively dividing cells, though these were not the most enriched categories. Shared stem cell genes comprise 979 orthology groups, 17% of which are enriched in all cell types of this family (Supplementary file 5). Importantly, these genes have consistently higher expression in the multipotent stem cells (MSCs) compared to lineage-restricted stem cells (LSCs) and differentiated cells (DCs) for all four species ( Figure 6D), and may represent a large, deeply conserved gene module specifically associated with multipotency.
Next, we asked whether there is an MSC-specific transcriptional regulatory program. Notably, we identified a number of transcriptional regulators enriched in MSCs across all four invertebrates. This included TFs from sox and klf orthology groups, the transcriptional coactivator sub1, and several chromatin modifiers, including ssrp1, a subunit of the FACT complex, and kat7, a histone acetyltransferase ( Figure 6D and Supplementary file 6). sox and klf are well-known pluripotency factors (Bialkowska et al., 2017;Sarkar and Hochedlinger, 2013), whereas the others have been studied in various processes associated with stem cell functions (MacPherson et al., 2020;Sikder et al., 2019;Zeng et al., 2013). The conserved enrichment of these transcriptional regulators in invertebrate MSCs suggests that their functional role in stem cells may be more phylogenetically ancient than previously appreciated. Determining their roles in establishing and maintaining multipotency across diverse animal taxa is an important avenue of future research.

Discussion
Cell types evolve as their gene expression programs change either as integrated units or via evolutionary splitting that results in separate derived programs. While this notion of coupled cellular and molecular evolution has gained significant traction in the past years, systematically comparing cell type-specific gene expression programs across species has remained a challenging problem. Here, we map single-cell atlases between evolutionarily distant species in a manner that accounts for the complexity of gene evolution. SAMap aligns cell atlases in two mutually reinforcing directions, mapping both the genes and the cells, with each feeding back into the other. Although this algorithm scales well to the majority of presently available whole-organism cell atlases, the size of single-cell datasets will continue to increase. The most memory-intensive steps in SAMap are the neighborhood coarsening of cross-species edges and the cross-species imputation of gene expression to calculate gene-gene correlations, both of which could be intractable for datasets containing millions of cells. Our current solution is to chunk these operations into smaller blocks for large datasets to avoid memory limitations, but the runtime increases significantly as a result. An option for mapping massive datasets may be to downsample each atlas. So long as all cell types are preserved and remain separable during downsampling, we expect the mapping results to be the same.
SAMap allows us to identify one-to-one cell type concordance between animals in the same phylum, whereas between phyla, we observe interconnected cell types forming distinct families. These findings support the notion that cell types evolve via hierarchical diversification (Arendt et al., 2019), resulting in cell type families composed of evolutionarily related cell types sharing a regulatory gene expression program that originated in their common ancestor. One-to-one cell type homologies should exist only if no further cell type diversification has occurred since the speciation. To understand the genetic underpinnings of the observed cell type homologies, we have examined similarities in gene expression programs of several cell types in depth. Among various gene sets, we have focused on TFs, as they form regulatory networks that specify cell type identities and activate downstream differentiation gene batteries (Erwin and Davidson, 2009). Notably, many of the aligned cell types we identified share expression of transcription factors known to play important roles in cell type differentiation, suggesting SAMap alignments based on mutual connectivity reflect evolutionary homology, rather than convergent functional similarity. SAMap thus provides a roadmap for tracing cell type evolutionary history and identifying the molecular changes in transcription factor regulatory complexes that have driven cell type diversification.
In parallel, SAMap systematically identifies instances where paralogs exhibit greater expression similarity than orthologs across species. Paralog substitution can occur due to differential loss or retention of cell-type-specific expression patterns of genes that were duplicated in the common ancestor (Shafer et al., 2020;Studer and Robinson-Rechavi, 2009) or due to compensating upregulation of paralogs following a loss-of-function mutation acquired by an ortholog (El-Brolosy et al., 2019). Considering our observation that paralog substitutions occur at higher rates for more recent paralogs, which should be more capable of functionally compensating for each other, we expect the latter scenario to be more likely, at least between frog and zebrafish. Paralog substitutions may also play an important role in cell type diversification, enabling newly evolved sister cell types to subfunctionalize via the use of distinct paralogs. Whereas the analysis presented here focuses on comparisons between two species, incorporating multiple species into a single analysis that also accounts for their phylogenetic relatedness could enable determining the stability of paralog substitutions within clades and their associated cell type diversification events. However, this will require datasets that densely sample species within specific clades and at key branching points along the tree of life.
Besides applications in evolutionary biology, we anticipate SAMap can catalyze the annotation of new cell atlases from non-model organisms, which often represents a substantial bottleneck requiring extensive manual curation and prior knowledge. Its ability to use the existing atlases to inform the annotation of cell types in related species will keep improving as more datasets become available to better sample the diversity of cell types throughout the animal kingdom.

Data and code availability
The source code for SAMap is publicly available at Github (https://github.com/atarashansky/SAMap; copy archived at swh:1:rev:c696585f8fe41ec1599b0720df579f3cb14f935b; Tarashansky et al., 2021), along with the code to perform the analysis and generate the types of plots presented in the figures. We also provide a wrapper function to launch a graphical user interface provided by the SAM package to interactively explore both datasets in the combined manifold. The datasets analyzed in this study are detailed in Supplementary file 1 with their accessions and annotations provided.

The SAMap algorithm
The SAMap algorithm contains three major steps: preprocessing, mutual nearest neighborhood alignment, and gene-gene correlation initialization. The latter two are repeated for three iterations, by default, to balance alignment performance and computational runtime.

Preprocessing Generate gene homology graph via reciprocal BLAST
We first construct a gene-gene bipartite graph between two species by performing reciprocal BLAST of their respective transcriptomes using tblastx, or proteomes using blastp. tblastn and blastx are used for BLAST between proteome and transcriptome. When a pair of genes share multiple High Scoring Pairs (HSPs), which are local regions of matching sequences, we use the HSP with the highest bit score to measure homology. Only pairs with E-value <10 À6 are included in the graph.
Although we define similarity using BLAST, SAMap is compatible with other protein homology detection methods (e.g. HMMER [Eddy, 2008]) or orthology inference tools (e.g. OrthoClust [Yan et al., 2014] and eggNOG [Huerta-Cepas et al., 2019]). While each of these methods has known strengths and limitations, BLAST is chosen for its broad usage, technical convenience, and compatibility with low-quality transcriptomes.
We encode the BLAST results into two triangular adjacency matrices, A and B, each containing bit scores in one BLAST direction. We combine A and B to form a gene-gene adjacency matrix G. After symmetrizing G, we remove edges that only appear in one direction: where Recip only keeps reciprocal edges, and m 1 and m 2 are the number of genes of the two species, respectively. To filter out relatively weak homologies, we also remove edges where G ab <0:25 b max G ab ð Þ. Edge weights are then normalized by the maximum edge weight for each gene and transformed by a hyperbolic tangent function to increase discriminatory power between low and high edge weights, Construct manifolds for each cell atlas separately using the SAM algorithm The single-cell RNAseq datasets are normalized such that each cell has a total number of raw counts equal to the median size of single-cell libraries. Gene expressions are then log-normalized with the addition of a pseudocount of 1. Genes expressed (i.e. log 2 D þ 1 ð Þ>1) in greater than 96% of cells are filtered out. SAM is run using the following parameters: preprocessing = 'StandardScaler', weight_PCs = False, k = 20, and npcs = 150. A detailed description of parameters is provided previously (Tarashansky et al., 2019). SAM outputs N 1 and N 2 , which are directed adjacency matrices that encode k-nearest neighbor graphs for the two datasets, respectively.
SAM only includes the top 3000 genes ranked by SAM weights and the first 150 principal components (PCs) in the default mode to reduce computational complexity. However, downstream mapping requires PC loadings for all genes. Thus, in the final iteration of SAM, we run PCA on all genes and take the top 300 PCs. This step generates a loading matrix for each species i, L i 2 R 300Âmi .

Mutual nearest neighborhood alignment Transform feature spaces between species
For the gene expression matrices Z i 2 R niÂmi , where n and m are the number of cells and genes respectively, we first zero the expression of genes that do not have an edge in G and standardize the expression matrices such that each gene has zero mean and unit variance, yielding Z i . G represents a bipartite graph in the form of matrix and H is the biadjacency matrix. Letting H 1 ¼ H and H 2 ¼ H T encoding directed edges from species 1 to 2 and 2 to 1, respectively, we normalize the biadjacency matrix H i such that each row sums to 1: H i ¼ SumNorm H i ð Þ 2 R miÂmj , where the SumNorm function normalizes the rows to sum to 1. The feature spaces can be transformed between the two species via weighted averaging of gene

Project single-cell gene expressions into a joint PC space
We project the expression data from two species into a joint PC space (Barkas et al., 2019), We then horizontally concatenate the principal components P i and P ij to form P i 2 R niÂ600 .

Calculate k-nearest cross-species neighbors for all cells
Using the joint PCs, P i , we identify for each cell the k-nearest neighbors in the other dataset using cosine similarity (k ¼ 20 by default). Neighbors are identified using the hnswlib library, a fast approximate nearest-neighbor search algorithm (Malkov and Yashunin, 2020). This outputs two directed biadjacency matrices C i 2 R niÂnj for i; j ð Þ ¼ 1; 2 ð Þ or ð2; 1Þ with edge weights equal to the cosine similarity between the PCs.
Apply the graph-coarsening mapping kernel to identify cross-species mutual nearest neighborhoods To increase the stringency and confidence of mapping, we only rely on cells that are mutual nearest cross-species neighbors, which are typically defined as two cells reciprocally connected to one another (Haghverdi et al., 2018). However, due to the noise in cell-cell correlations and stochasticity in the kNN algorithms, cross-species neighbors are often randomly assigned from a pool of cells that appear equally similar, decreasing the likelihood of mutual connectivity between individual cells even if they have similar expression profiles. To overcome this limitation, we integrate information from each cell's local neighborhood to establish more robust mutual connectivity between cells across species. Two cells are thus defined as mutual nearest cross-species neighbors when their respective neighborhoods have mutual connectivity.
Specifically, the nearest neighbor graphs N i generated by SAM are used to calculate the neighbors of cells t i hops away along outgoing edges: where N i are adjacency matrices that contain the number of paths connecting two cells t i hops away, for i ¼ 1 or 2. t i determines the lengthscale over which we integrate incoming edges for species i. Its default value is 2 if the dataset size is less than 20,000 cells and 3 otherwise. However, cells within tight clusters may have spurious edges connecting to other parts of the manifold only a few hops away. To avoid integrating neighborhood information outside this local structure, we use the Leiden algorithm (Traag et al., 2019) to cluster the graph and identify a local neighborhood size for each cell (the resolution parameter is set to 3 by default). If cell a belongs to cluster c a , then its neighborhood size is l a ¼ jc a j. For each row a in N i we only keep the l a geodesically closest cells, letting the pruned graph update N i . Edges outgoing from cell a i in species i are encoded in the corresponding row in the adjacency matrix: C i;ai . We compute the fraction of the outgoing edges from each cell that target the local neighborhood of a cell in the other species: C i;aibj ¼ P c2Xj;b j C i;aic , where X j;bj is the set of cells in the neighborhood of cell b j in species j and C i;aibj is the fraction of outgoing edges from cell a i in species i targeting the neighborhood of cell b j in species j.
To reduce the density of C i so as to satisfy computational memory constraints, we remove edges with weight less than 0.1. Finally, we apply the mutual nearest neighborhood criterion by taking the element-wise, geometric mean of the two directed bipartite graphs: . This operation ensures that only bidirectional edges are preserved, as small edge weights in either direction results in small geometric means.
Assign the k-nearest cross-species neighborhoods for each cell Given the mutual nearest neighborhoods C 2 R n1Ân2 , we select the k nearest neighborhoods for each cell in both directions to update the directed biadjacency matrices C 1 and C 2 : C 1 ¼ KNN C ; k and

Stitch the manifolds
We use C 1 and C 2 to combine the manifolds N 1 and N 2 into a unified graph. We first weight the edges in N 1 and N 2 to account for the number of shared cross-species neighbors by computing the one-mode projections of C 1 and C 2 . In addition, for cells with strong cross-species alignment, we attenuate the weight of their within-species edges. For cells with little to no cross-species alignment, their within-species are kept the same to ensure that the local topological information around cells with no alignment is preserved. Specifically, we use N 1 and N 2 to mask the edges in the one-mode projections, where U E ð Þ sets all edge weights in graph E to 1 and Norm normalizes the outgoing edges from each cell to sum to 1. The minimum edge weight is set to be 0.3 to ensure that neighbors in the original manifolds with no shared cross-species neighbors still retain connectivity: Ñ 1;ij ¼ min 0:3; Ñ 1;ij and Ñ 2;ij ¼ min 0:3; Ñ 2;ij for all edges i; j ð Þ. We then scale the within-species edges from cell i by the total weight of its cross-species edges: Finally, the within-and cross-species graphs are stitched together to form the combined nearest neighbor graph N: The overall alignment score between species 1 and 2 is defined as S ¼ 1

Homology graph refinement Update edge weights in the gene-gene bipartite graph with expression correlations
To compute correlations between gene pairs, we first transfer expressions from one species to the other: Z i;nimj ¼ C i;ni Z j;mj , where Z i;nimj is the imputed expressions of gene m j from species j for cell n i in species i, and C i;ni is row n i of the biadjacency matrix encoding the cross-species neighbors of cell n i in species i, all for i; j ð Þ ¼ 1; 2 ð Þ and ð2; 1Þ. We similarly use the manifolds constructed by SAM to smooth the within-species gene expressions using kNN averaging: Z j;mj ¼ N j;mj Z j;mj , where N j is the nearest-neighbor graph for species j. We then concatenate the within-and cross-species gene expressions such that the expression of gene m j from species j in both species is For all gene pairs in the initial unpruned homology graph, G, we compute their correlations, where 0 ð Þ is a Heaviside step function centered at 0 to set negative correlations to zero. We then use the expression correlations to update the corresponding edge weights in G, which are again normalized through G ab ¼ 0:5 þ 0:5tanh 10G ab = b max G ab À 5 .

Annotation of cell atlases
To annotate the primary zebrafish and Xenopus cell types, the cell subtype annotations provided by the original publications Wagner et al., 2018) are coarsened using a combination of the manual matching and developmental hierarchies. For example, as 'involuting marginal zone' in Xenopus is manually matched to 'non-dorsal margin', 'dorsal margin' 'non-dorsal margin involuted', and 'dorsal margin involuted' in zebrafish, we label these cells as 'involuting marginal zone'. In cases where the matching is insufficient to coarsen the annotations, we use the provided developmental trees to name a group of terminal cell subtypes by their common ontogenic ancestor. Cell types that do not cluster well in the manifold reconstructed by SAM are excluded from the comparison. These include germline, heart, and olfactory placode cells, as they are mixed with other cell types in the Xenopus atlas. The germline cells are scattered across the reconstructed manifold and do not concentrate in a distinct cluster. The heart cells and olfactory placode cells are inextricably mixed with larger populations of intermediate mesoderm and placodal cells, respectively. Similarly, the iridoblast, epiphysis, nanog + , apoptotic-like, and forerunner cells are excluded because they do not cluster distinctly in the zebrafish atlas. The annotations provided by their respective studies are used to label the cells in the Spongilla, Hydra, planarian, and mouse atlases. To annotate the schistosome cells, we use known marker genes to annotate the main schistosome tissue types . Annotations for all single cells in all datasets are provided in Supplementary file 1.

Visualization
The combined manifold N is embedded into 2D projections using UMAP implemented in the scanpy package  by scanpy.tl.umap with the parameter min_dist = 0.1. The sankeyD3 package (https://rdrr.io/github/fbreitwieser/sankeyD3/man/sankeyD3-package.html) in R is used to generate the sankey plots. Edge thickness corresponds to the alignment score between mapped cell types. The alignment score between cell types a and b is defined as where c a and c b are the set of cells in cell types a and b, respectively. Cell type pairs with alignment score less than z are filtered out. By default, z is set to be 0.1. The network graphs in Figure 5C are generated using the networkx package (https://networkx. github.io) in python. To focus on densely connected cell type groups, we filter out cell type pairs with alignment score less than 0.05.

Identification of gene pairs that drive cell type mappings
We define g 1 and g 2 to contain SAMap-linked genes from species 1 and 2, respectively. Note that a gene may appear multiple times as SAMap allows for one-to-many homology. Let X a1b2 denote the set of all cells with cross species edges between cell types a 1 and b 2 . We calculate the average standardized expression of all cells from species i that are in X a1b2 : where Z i;x;gi 2 R jgij is the standardized expression of genes g i in cell x. The correlation between Y 1;g1 and Y 2;g2 can be written as Corr Y 1;g1 ; Y 2;g2 have zero mean and unit variance. We use the summand to identify gene pairs that contribute most positively to the correlation. We assign each gene pair a score: T Z ð Þ sets negative values in vector Z to zero in order to ignore lowly-expressed genes. To be inclusive, we begin with the top 1000 gene pairs according to h g and filter out gene pairs in which one or both of the genes are not differentially expressed in their respective cell types (p-value > 10 À2 ), have less than 0.2 SAM weight, or are expressed in fewer than 5% of the cells in the cluster. The differential expression of each gene in each cell type is calculated using the Wilcoxon rank-sum test implemented in the scanpy function scanpy.tl.rank_genes_groups.

Orthology group assignment
We use the eggNOG mapper (v5.0) (Huerta-Cepas et al., 2019) to assign each gene to an orthology group with default parameters. For the zebrafish-to-Xenopus mapping, genes are considered orthologs if they map to the same vertebrate orthology group. For the pan-species analysis, we group genes from all species with overlapping orthology assignments. In Figure 6B, each column corresponds to one of these groups. As each group may contain multiple genes from each species, we present the expression of the gene with the highest enrichment score per species. All gene names and corresponding orthology groups are reported in Supplementary file 5.

Paralog substitution analysis
SAMap outputs gene-gene correlations across the combined manifold for all pairs of genes in the homology graph. As determined by eggNOG, genes that map to the same orthology group for the two species' most recent common ancestor are considered orthologs, and those that map to the same orthology group more ancestral than Vertebrata are considered as paralogs. We note that as eggNOG does not provide an orthology group corresponding to the osteichtyan ancestor, our analysis does not include the paralogs that duplicated in between the osteichtyan and the vertebrate ancestors. If a gene has significantly higher correlation to one of its paralogs than its ortholog (>0.3 by default), we consider its ortholog to have been substituted. Paralog substitutions are identified using the samap.analysis.ParalogSubstitutions function provided by the SAMap package.
The evolutionary time period in which paralogs were duplicated can be inferred by identifying their most recent shared orthology group. We calculate the enrichment of paralog substitutions for each taxonomic level (i.e. Chordata, Bilateria, Metazoa, Opisthokonta, and Eukaryota) using the egg-NOG orthology group assignments. We normalize the number of substituting paralogs by the total number of paralogs at each level to calculate the rate of paralog substitution across evolutionary time.
To quantify the enrichment of substituting paralogs in each cell type, we define a cell type-specific substitution score. We first assign paralog substitution events to cell types if the paralogous gene pairs are enriched in any of their mappings. Each cell type k then has a set of substituting paralogs P k . The score S k for cell type k is calculated as S k ¼ P i2Pk 1Àni mk , where n i is the number of paralogs of ortholog i normalized by the maximum number of paralogs observed across all genes to accounts for the fact that genes with more paralogs are more likely to match with substituting paralogs by random chance, and m k is the number of differentially expressed genes in cell type k. Similarly, the denominator accounts for the fact that cell types with more differentially expressed genes are more likely to have paralog substitutions by random chance. The substitution scores for cell types with annotated homologs across species are averaged.

Phylogenetic reconstruction of gene trees
We generate gene trees to validate the identity of genes involved in putative examples of paralog substitution and of Fox and Csrp transcriptional regulators that are identified as enriched in contractile cells. For this, we first gather protein sequences from potential homologs using the eggnog version 5.0 orthology database (Huerta-Cepas et al., 2019). For the Fox and Csrp phylogenies, we include all Fox clade I (Larroux et al., 2008) and Csrp/Crip homologs, respectively, from the seven species included in our study.
Alignment of protein sequences is performed with Clustal Omega version 1.2.4 using default settings as implemented on the EMBL EBI web services platform (Madeira et al., 2019). Maximum likelihood tree reconstruction is performed using IQ-TREE version 1.6.12 (Nguyen et al., 2015) with the ModelFinder Plus option (Kalyaanamoorthy et al., 2017). For the Csrp tree, we perform 1000 nonparametric bootstrap replicates to assess node support. For Fox, we utilize the ultrafast bootstrap support option with 1000 replicates. For each gene tree we choose the model that minimizes the Bayesian Information Criterion (BIC) score in ModelFinder. This results in selection of the following models: DCMut+R4 (Csrp) and VT+F + R5 (Fox). The final consensus trees are visualized and rendered using the ETE3 v3.1.1 python toolkit (Huerta-Cepas et al., 2016) and the Interactive Tree of Life v4 .

KOG functional annotation and enrichment analysis
Using the eggNOG mapper, KOG functional annotations are transferred to individual transcripts from their assigned orthology group. For enrichment analysis, all genes enriched in the set of cell type pairs of interest are lumped to form the target set for each species. For example, the target set for Spongilla archaeocytes used in Figure 6C is composed of all genes enriched between Spongilla archaeocytes and other invertebrate stem cells. Note that this set includes genes from other species that are linked by SAMap to the Spongilla archeocyte genes. We include genes from other species in the target set to account for differences in KOG functional annotation coverage between species. As such, the annotated transcripts from all seven species are combined to form the background set. We use a hypergeometric statistical test (Eden et al., 2009) to measure the enrichment of the KOG terms in the target genes compared to the background genes.

Mapping zebrafish and Xenopus atlases using existing methods
For benchmarking, we use vertebrate orthologs as determined by eggNOG as input to Harmony (Korsunsky et al., 2019), LIGER (Welch et al., 2019), Seurat (Stuart et al., 2019), Scanorama (Hie et al., 2019), BBKNN (Polański et al., 2019), which are all run with default parameters. One-toone orthologs are selected from one-to-many and many-to-many orthologs by using the bipartite maximum weight matching algorithm implemented in networkx. When using the one-to-one orthologs as input for SAMap, we run for only one iteration. The resulting integrated lower-dimensional coordinates (PCs for Seurat, Harmony, and Scanorama and non-negative matrix factorization coordinates for LIGER) and stitched graphs (BBKNN and SAMap) are all projected into 2D with UMAP (Figure 2-figure supplement 1A). The integrated coordinates are used to generate a nearest neighbor graph using the correlation distance metric, which is then used to compute the alignment scores in Figure 2-figure supplement 1B. The alignment scores for SAMap and BBKNN are directly computed from their combined graphs.
In situ hybridization in schistosomes S. mansoni (strain: NMRI) juveniles are retrieved from infected female Swiss Webster mice (NR-21963) at~3 weeks post-infection by hepatic portal vein perfusion using 37˚C DMEM supplemented with 5% heat inactivated FBS. The infected mice are provided by the NIAID Schistosomiasis Resource Center for distribution through BEI Resources, NIH-NIAID Contract HHSN272201000005I. In adherence to the Animal Welfare Act and the Public Health Service Policy on Humane Care and Use of Laboratory Animals, all experiments with and care of mice are performed in accordance with protocols approved by the Institutional Animal Care and Use Committees (IACUC) of Stanford University (protocol approval number 30366). In situ hybridization experiments are performed as described previously (Tarashansky et al., 2019), using riboprobes synthesized from gene fragments