Tara Chari

Brandon Weissbourd

Jase Gehring

Anna Ferraioli

Lucas Leclère

Makenna Herl

Fan Gao

Sandra Chevalier

Richard R Copley
email: copley@obs-vlfr.fr

Evelyn Houliston
email: houliston@obs-vlfr.fr

David J Anderson

Lior Pachter
email: lpachter@caltech.edu

Whole Animal Multiplexed Single-Cell RNA-Seq Reveals Transcriptional Shifts Across Clytia Medusa Cell Types

We present an organism-wide, transcriptomic cell atlas of the hydrozoan medusa Clytia hemisphaerica and describe how its component cell types respond to perturbation. Utilizing multiplexed scRNA-seq, in which individual animals were indexed and pooled from control and perturbation conditions into a single sequencing run, we avoid artifacts from batch effects and are able to discern shifts in cell state in response to organismal perturbations. This work serves as a foundation for future studies of development, function, and regeneration in a genetically tractable jellyfish species. Moreover, we introduce a powerful workflow for high-resolution, whole animal, multiplexed single-cell genomics (WHAM-seq) that is readily adaptable to other traditional or non-traditional model organisms.

Introduction

Single-cell RNA sequencing (scRNA-seq) is enabling the survey of complete transcriptomes of thousands to millions of cells [START_REF] Hwang | Single-cell RNA sequencing technologies and bioinformatics pipelines[END_REF], resulting in the establishment of cell atlases across whole organisms [START_REF] Vergara | Whole-organism cellular gene-expression atlas reveals conserved cell types in the ventral nerve cord of Platynereis dumerilii[END_REF][START_REF] Sebé-Pedrós | Cnidarian Cell Type Diversity and Regulation Revealed by Whole-Organism Single-Cell RNA-Seq[END_REF][START_REF] Cao | Comprehensive single-cell transcriptional profiling of a multicellular organism[END_REF][START_REF] Han | Mapping the Mouse Cell Atlas by Microwell-Seq[END_REF][START_REF] Wagner | Singlecell mapping of gene expression landscapes and lineage in the zebrafish embryo[END_REF], exploration of the diversity of cell types throughout the animal kingdom [START_REF] Sebé-Pedrós | Cnidarian Cell Type Diversity and Regulation Revealed by Whole-Organism Single-Cell RNA-Seq[END_REF][START_REF] Fincher | Cell type transcriptome atlas for the planarian Schmidtea mediterranea[END_REF][START_REF] Siebert | Stem cell differentiation trajectories in Hydra resolved at single-cell resolution[END_REF][START_REF] Gerber | Single-cell analysis uncovers convergence of cell identities during axolotl limb regeneration[END_REF], and investigation of transcriptomic changes under perturbation [START_REF] Dixit | Perturb-Seq: Dissecting Molecular Circuits with Scalable Single-Cell RNA Profiling of Pooled Genetic Screens[END_REF][START_REF] Mcfarland | Multiplexed single-cell transcriptional response profiling to define cancer vulnerabilities and therapeutic mechanism of action[END_REF]. However, scRNA-seq studies involving multiple samples can be costly, and may be confounded by batch effects resulting from multiple distinct library preparations [START_REF] Butler | Integrating single-cell transcriptomic data across different conditions, technologies, and species[END_REF][START_REF] Tran | A benchmark of batch-effect correction methods for single-cell RNA sequencing data[END_REF]. Recent developments in scRNA-seq multiplexing technology expand the number of samples, individuals, or perturbations that can be incorporated within runs, facilitating well-controlled scRNA-seq experiments [START_REF] Mcfarland | Multiplexed single-cell transcriptional response profiling to define cancer vulnerabilities and therapeutic mechanism of action[END_REF][START_REF] Gehring | Highly multiplexed single-cell RNA-seq by DNA oligonucleotide tagging of cellular proteins[END_REF][START_REF] Mcginnis | MULTI-seq: sample multiplexing for single-cell RNA sequencing using lipid-tagged indices[END_REF][START_REF] Guo | CellTag Indexing: genetic barcode-based sample multiplexing for single-cell genomics[END_REF][START_REF] Srivatsan | Massively multiplex chemical transcriptomics at single-cell resolution[END_REF][START_REF] Stoeckius | CITE-seq. Protocol Exchange[END_REF]. These advances have created an opportunity to explore systems biology of whole organisms at single-cell resolution, merging the concepts of cell atlas surveys with multiplexed single-cell experimentation.

Here we apply this powerful experimental paradigm to a planktonic model organism. We examine the medusa (free-swimming jellyfish) stage of the hydrozoan Clytia hemisphaerica, with dual motivations. Firstly, Clytia is a powerful, emerging model system spanning multiple fields, from evolutionary and developmental biology to regeneration and neuroscience [START_REF] Sinigaglia | Pattern regulation in a regenerating jellyfish[END_REF][START_REF] Bosch | Back to the Basics: Cnidarians Start to Fire[END_REF][START_REF] Leclère | The genome of the jellyfish Clytia hemisphaerica and the evolution of the cnidarian life-cycle[END_REF][START_REF] Leclère | Hydrozoan insights in animal development and evolution[END_REF][START_REF] Steinmetz | Independent evolution of striated muscles in cnidarians and bilaterians[END_REF][START_REF] Kamran | In vivo imaging of epithelial wound healing in the cnidarian Clytia hemisphaerica demonstrates early evolution of purse string and cell crawling closure mechanisms[END_REF]. While previous work has characterized a number of cell types in the Clytia medusa [START_REF] Leclère | The genome of the jellyfish Clytia hemisphaerica and the evolution of the cnidarian life-cycle[END_REF], a wholeorganism atlas of transcriptomic cell types has been lacking. Such an atlas is a critical resource for the Clytia community, and an important addition to the study of cell types across animal phylogeny.

Secondly, emerging multiplexing techniques present new opportunities for systems-level studies of cell types and their changing states at unprecedented resolution in whole organisms. The Clytia medusa offers an appealing platform for pioneering such studies. It is small, transparent, and has simple tissues and organs, stem cell populations actively replenishing many cell types in mature animals, and remarkable regenerative capacity [START_REF] Sinigaglia | Pattern regulation in a regenerating jellyfish[END_REF][START_REF] Leclère | Hydrozoan insights in animal development and evolution[END_REF][START_REF] Kamran | In vivo imaging of epithelial wound healing in the cnidarian Clytia hemisphaerica demonstrates early evolution of purse string and cell crawling closure mechanisms[END_REF][START_REF] Galliot | Cnidarians as a model system for understanding evolution and regeneration[END_REF][START_REF] Amiel | Clytia hemisphaerica: a cnidarian model for studying oogenesis. Oogenesis: the universal process[END_REF][START_REF] Leclère | Diversity of Cnidarian Muscles: Function, Anatomy, Development and Regeneration[END_REF]. Furthermore, the 1cm-diameter adult medusae used in this study contain on the order of 10 5 cells, making it possible to sample cells comprehensively across a whole animal in a cost-effective manner using current scRNA-seq technology (Table S1,S2; Fig. S1). In this study, we generate a cell atlas for the Clytia medusa while simultaneously performing a whole organism perturbation study, providing the first medusa single-cell dataset and an examination of changing cell states across the organism. Our work also provides a proof-of-principle for perturbation studies in nontraditional model organisms, using multiplexing technology and a reproducible workflow with lessened reliance on functional annotation, from the experimental implementation to the data processing and analysis.

Results

We compared control versus starved animals, as this strong, naturalistic stimulus was likely to cause significant, interpretable changes in transcription across multiple cell types. Laboratory-raised, young adult, female medusae were split into two groups of five animals, one deprived of food for four days, and the second fed daily (see Methods). We observed numerous phenotypic changes in starved animals, including a dramatic size reduction reflecting two-to three-fold fewer cells (28) (Fig. 1; Fig. S2; see Methods), and a striking reduction in gonad size. Correspondingly, the number of eggs released per day decreased (29) (Fig. S3).

For scRNA-seq, single-cell suspensions were prepared from each whole medusa and individually labeled with unique ClickTag barcodes ( 14) using a sea-water compatible workflow (see Methods; Table S2,S3; Supplementary Methods; Fig. S4). All labeled suspensions were pooled and processed with the 10X Genomics V2.0 workflow and Illumina sequencing, allowing construction of a combined dataset across organisms and treatments, without requiring batch correction (Fig. 1; Table S1; Fig. S4A-D). A total of 13,673 single-cell profiles derived from 10 individuals (5 control, 5 starved) passed quality control, with high concordance in cell type abundance and gene expression among animals in the same treatment condition (see Methods; Fig. S5). From this gene expression matrix, we 1) derived a Clytia medusa cell atlas, and 2) generated a high-resolution resource of the transcriptional impact of starvation across all observed cell types.

In order to validate the cell atlas, and assess technical variability across and within multiplexed experiments, we performed an additional, independent round of sequencing from twelve individuals. We found that cell types were highly concordant between experiments, and confirmed a reduction of batch effect-induced variability within multiplexed experiments (see below). During this second sequencing run, we took advantage of our multiplexing approach to perform an experiment designed both to search for transcripts with "immediate early gene (IEG)"like behavior in Clytia and to test its sensitivity for detecting more rapid or subtle gene changes than those of the extreme starvation perturbation. For this we exposed Clytia medusae to multiple transient, ionic stimuli and dissociated ~1h later. This paradigm allowed us to identify candidate genes with IEG-like properties across many cell types, including neurons (30) (Fig. S6-S8; Table S4). IEGs are valuable tools in neuroscience, to identify neurons that are active following a specific stimulus or behavior [START_REF] Sheng | The regulation and function of c-fos and other immediate early genes in the nervous system[END_REF]. This methodology is thus able to detect transcriptional responses across diverse stimulus-response paradigms (Table S5).

A Clytia Cell Atlas

To generate the cell atlas, we clustered the cells using the gene expression matrix, extracting 36 cell types and their corresponding marker genes (Fig. 2a,b; Fig. S9-S11; Table S5; see Methods). Each of the cell types was present in each of the individual animals sequenced (Fig. S12). We then generated a low-dimensional representation [START_REF] Wolf | PAGA: graph abstraction reconciles clustering with trajectory inference through a topology preserving map of single cells[END_REF][START_REF] Mcinnes | UMAP: Uniform Manifold Approximation and Projection for Dimension Reduction[END_REF] of these cell types (Fig. 2A). We could group the cell types into seven broad classes (Fig. 2A) which correspond to the outer epidermis, the inner gastrodermis, and to likely derivatives of the multipotent interstitial stem cell population (i-cells). I-cells are a specific feature of hydrozoans, and are particularly well characterized in Hydra, where they generate neural cells, gland cells, and stinging cells (nematocytes), as well as germ cells [START_REF] Siebert | Stem cell differentiation trajectories in Hydra resolved at single-cell resolution[END_REF][START_REF] Bosch | Back to the Basics: Cnidarians Start to Fire[END_REF][START_REF] Hemmrich | Molecular signatures of the three stem cell lineages in hydra and the emergence of stem cell function at the base of multicellularity[END_REF]. Our dataset was derived from female medusae so it lacks male germ cells, and late stage oocytes are expected to be too large for capture by the dissociation procedure.

The 36 cell types (see Methods; Fig. 2b-d) were concordant between the two separate multiplexed experiments ('Starvation' and 'Stimulation'), and robust to different transcriptome annotations (Fig. S6,S13). For some of them, cell type identity could be assigned on the basis of published information on gene expression in Clytia and/or of homologous genes in other animals, while for the others we performed in situ hybridization for selected marker genes (Fig. 2C; Fig. S11,S14; Table S3). Previously known cell types apparent in our data included i-cells [START_REF] Leclère | Maternally localized germ plasm mRNAs and germ cell/stem cell formation in the cnidarian Clytia[END_REF] and nematocytes at successive stages of differentiation [START_REF] Condamine | Molecular characterisation of a cellular conveyor belt in Clytia medusae[END_REF][START_REF] Denker | Ordered progression of nematogenesis from stem cells through differentiation stages in the tentacle bulb of Clytia hemisphaerica (Hydrozoa, Cnidaria)[END_REF][START_REF] Sunagar | Cell type-specific expression profiling unravels the development and evolution of stinging cells in sea anemone[END_REF], as well as oocytes [START_REF] Takeda | Identification of jellyfish neuropeptides that act directly as oocyte maturation-inducing hormones[END_REF], gonad epidermis, manubrium epidermis, and bioluminescent cells in the tentacles that each express specific endogenous green fluorescent proteins (GFPs) [START_REF] Fourrage | An endogenous green fluorescent protein--photoprotein pair in Clytia hemisphaerica eggs shows co-targeting to mitochondria and efficient bioluminescence energy transfer[END_REF].

In situ hybridization for a selection of diagnostic muscle cell type genes allowed us to describe cell types making up the smooth and striated muscles, for instance, distinguishing the striated muscle cells lining the bell (subumbrella) and velum (Fig. 2C,D; Fig. S14) [START_REF] Steinmetz | Independent evolution of striated muscles in cnidarians and bilaterians[END_REF][START_REF] Leclère | Diversity of Cnidarian Muscles: Function, Anatomy, Development and Regeneration[END_REF]. Within known cell types, clustering revealed an unappreciated degree of cell heterogeneity, yielding novel subtypes. For example, eight cell types could be distinguished within the gastrodermis, six of which were designated gastro-digestive (GD A-F) on the basis of a largely shared set of marker genes (Fig. 2B), including enzymes associated with intracellular digestion, such as CathepsinL [START_REF] Steinmetz | A non-bilaterian perspective on the development and evolution of animal digestive systems[END_REF]. Unlike most other clusters the GD clusters differ primarily in their relative levels of gene expression, rather than by unique marker genes. They thus likely represent variations on a similar digestive-absorptive epithelial cell type with different functional specializations, distributed across the main digestive compartments of the gastrodermis -the manubrium, gonad, and tentacle bulb -as well as the gastrovascular canals that link them (Fig. S11,S14). Comparison of gene modules discriminating these gastrodermal clusters (Fig. S15) indicates that GD-B may have a particular role in TGF-β signaling, likely involving the ligands BMP2/4 and BMP5/8. GD-D is enriched for a module associated with cell-cell junctions, while GD-F shows relative depletion, suggesting poorer integration into the gastrodermal epithelium (Fig. S15) and possible involvement in GD cell mobilization during starvation and regeneration [START_REF] Sinigaglia | Pattern regulation in a regenerating jellyfish[END_REF]. GD-C cells, localized closest to the endodermal plate, are enriched for transcripts associated with extracellular matrix and mesoglea (jelly) production (Fig. 2C,D ; Fig. S11,S14). Expression of these and other genes implicated in mesoglea production, such as fibrillar collagens, is also a characteristic of endodermal plate cells (cluster 33) and proximal tentacle-bulb endoderm cells (cluster 16).

Digestive gland cells fell into five types expressing different mixtures of enzymes for extracellular digestion. These showed overlapping distributions in the mouth and stomach regions of the manubrium. Two subtypes of gland cells (type C and E) were also present within the gonad gastroderm. Four broad clusters corresponding to neural cells each appeared to represent mixed populations, and could be subdivided by further analyses to define 14 likely subpopulations of neurons (see below). Seven major clusters could be assigned identities as nematocytes at different developmental stages, comprising two groups with highly distinct transcriptional signatures. Four of these we designate "nematoblasts" on the basis of high levels of transcripts related to formation of the nematocyst (stinging capsule) [START_REF] Condamine | Molecular characterisation of a cellular conveyor belt in Clytia medusae[END_REF][START_REF] Denker | Ordered progression of nematogenesis from stem cells through differentiation stages in the tentacle bulb of Clytia hemisphaerica (Hydrozoa, Cnidaria)[END_REF][START_REF] Denker | Horizontal gene transfer and the evolution of cnidarian stinging cells[END_REF]. The other three, designated as differentiating and mature "nematocytes", show no enrichment of these nematocyst transcripts, but strongly express highly conserved proteins of the actin-rich "stereovilli" of vertebrate hair cells, including Whirlin, Harmonin, and Sans/USH-IG. This is consistent with observations of similar actin-based protrusions surrounding a central cilium in many of the mechanosensory cell types described in other cnidarian species [START_REF] Bezares-Calderón | Diversity of cilia-based mechanosensory systems and their functions in marine animal behaviour[END_REF]. Related but more elaborate actin structures are associated with the cnidocil of mature nematocytes, but it had not previously been known to share functional hair-cell components [START_REF] Bezares-Calderón | Diversity of cilia-based mechanosensory systems and their functions in marine animal behaviour[END_REF][START_REF] Mcpherson | Sensory Hair Cells: An Introduction to Structure and Physiology[END_REF]. Nematocilin, a hydrozoan-specific component of the nematocil (ciliary trigger for nematocyte discharge) [START_REF] Hwang | Cilium evolution: identification of a novel protein, nematocilin, in the mechanosensory cilium of Hydra nematocytes[END_REF], is also expressed in these clusters (Table S5; Fig. S14; see below). In situ hybridizations revealed marker expression in morphologically distinguishable nematocytes, notably including two lines along the oral face of each tentacle (Fig. 3E; Fig. S14), a striking arrangement overlooked in previous studies.

A remarkable feature of the Clytia medusa is that it constantly generates many cell types, notably neural cells and nematocytes from prominent i-cell pools in the tentacle bulb epidermis [START_REF] Denker | Ordered progression of nematogenesis from stem cells through differentiation stages in the tentacle bulb of Clytia hemisphaerica (Hydrozoa, Cnidaria)[END_REF], as well as at other sites [START_REF] Leclère | Maternally localized germ plasm mRNAs and germ cell/stem cell formation in the cnidarian Clytia[END_REF]. Within our dataset, we thus expected to be able to capture dynamic information relating to the development of i-cell derived cell types, similar to that extracted from Hydra polyp single-cell transcriptome data [START_REF] Siebert | Stem cell differentiation trajectories in Hydra resolved at single-cell resolution[END_REF]. As in Hydra, our cell atlas revealed clear connections between the neuronal and nematocyte populations and the i-cell population (3, 8) (Fig. 2A; Fig. S11,S14), likely corresponding to differentiation trajectories [START_REF] Condamine | Molecular characterisation of a cellular conveyor belt in Clytia medusae[END_REF][START_REF] Denker | Ordered progression of nematogenesis from stem cells through differentiation stages in the tentacle bulb of Clytia hemisphaerica (Hydrozoa, Cnidaria)[END_REF]. In contrast, we found no clear developmental connection between i-cells and gland cells, and little to no expression of markers of the common neuronal-gland cell precursors identified in Hydra (8) (Fig. S16). In Hydra, gland cells are generated not only from the i-cell lineage but also by processes of self-renewal and position-dependent transdifferentiation [START_REF] Siebert | Stem cell differentiation trajectories in Hydra resolved at single-cell resolution[END_REF]. In the Clytia medusa, digestive gland cells show widespread distribution across distinct regions of the manubrium and gonad compartments of the gastrovascular system, spatially separated from i-cell populations positioned proximally in both these organs (Fig. S14). It is possible that these alternative pathways may dominate over direct differentiation from the i-cell lineage in this system.

To address the developmental relationships between the different neural and nematocyte clusters and identify developmental markers, we assigned pseudotime values to the cells and ranked genes in each trajectory (Fig. 3A, see Methods). This revealed trajectories consistent with these cell types both deriving from i-cells (8) (Fig. 3A,B). Examination of the nematocyte trajectory revealed the early expression of genes previously not associated with this process, including Znf845 and Mos3 (45) (Fig. 3C,D). Nematocyst-related genes such as Minicollagens, Polyglutamate synthases, Dkk3, and NOWA, were then expressed during a first major phase of nematogenesis, consistent with previous reports [START_REF] Condamine | Molecular characterisation of a cellular conveyor belt in Clytia medusae[END_REF][START_REF] Denker | Ordered progression of nematogenesis from stem cells through differentiation stages in the tentacle bulb of Clytia hemisphaerica (Hydrozoa, Cnidaria)[END_REF][START_REF] Sunagar | Cell type-specific expression profiling unravels the development and evolution of stinging cells in sea anemone[END_REF][START_REF] Denker | Horizontal gene transfer and the evolution of cnidarian stinging cells[END_REF]) (Fig. 3C,D; corresponding expression domains in Fig. 3E; Fig. S14). The trajectory analysis confirmed continuity between the 'nematoblast' clusters and the distinct and underappreciated nematocyte differentiation phase, characterized by expression of putative nematocil structural proteins and Nematocilin expression at the end of the trajectory (Fig. 3C,D; Fig. S11,S14) (see above). The two phases of nematogenesis were linked by the expression of rare specific marker genes for cluster 17 (e.g. M14 peptidase in Fig. 3E; Fig. S14). Consistent with this linking of the nematoblast and differentiation phases revealed in trajectory analysis, in situ markers showed distinct expression territories in the tentacle bulb and tentacle respectively (Fig. 3E). Furthermore, we found that markers of both phases and their respective orthologs, including the 'hair cell' gene set, were appropriately distributed amongst transcriptomes derived from dissected Clytia bulb and tentacle regions [START_REF] Condamine | Molecular characterisation of a cellular conveyor belt in Clytia medusae[END_REF], and between developing and mature nematocyte scRNA-seq clusters in Hydra [START_REF] Siebert | Stem cell differentiation trajectories in Hydra resolved at single-cell resolution[END_REF].

Cnidarian nervous systems represent both valuable points of phylogenetic comparison and tractable platforms for systems neuroscience [START_REF] Sebé-Pedrós | Cnidarian Cell Type Diversity and Regulation Revealed by Whole-Organism Single-Cell RNA-Seq[END_REF][START_REF] Siebert | Stem cell differentiation trajectories in Hydra resolved at single-cell resolution[END_REF][START_REF] Bosch | Back to the Basics: Cnidarians Start to Fire[END_REF]. However, the molecular heterogeneity of neural cell types and their developmental progression remains largely unexplored, particularly in the more complex medusa forms. We therefore extracted genes expressed during neural development which included those encoding bHLH, Sox, and other transcription factors with potential roles in neurogenesis or fate specification, and numerous other genes of interest in neuronal development, such as cell adhesion molecules [START_REF] Denker | Ordered progression of nematogenesis from stem cells through differentiation stages in the tentacle bulb of Clytia hemisphaerica (Hydrozoa, Cnidaria)[END_REF][START_REF] Massari | Helix-loop-helix proteins: regulators of transcription in eucaryotic organisms[END_REF][START_REF] Kurki | Expression of proliferating cell nuclear antigen (PCNA)/cyclin during the cell cycle[END_REF][START_REF] Sanes | Synaptic Specificity, Recognition Molecules, and Assembly of Neural Circuits[END_REF] (Fig. 4A-C; Table S5). In mature neurons, neuropeptides are thought to be the dominant neurotransmitters in cnidarians [START_REF] Jékely | The chemical brain hypothesis for the origin of nervous systems[END_REF][START_REF] Grimmelikhuijzen | Mini-review: the evolution of neuropeptide signaling[END_REF] but are challenging to identify due to rapid sequence evolution [START_REF] Jekely | Global view of the evolution and diversity of metazoan neuropeptide signaling[END_REF][START_REF] Nielsen | De novo transcriptome assembly of the cubomedusa Tripedalia cystophora, including the analysis of a set of genes involved in peptidergic neurotransmission[END_REF]. In conjunction with sequence-based analysis, we were able to identify ten new likely neuropeptides on the basis of their inclusion as marker genes for the four basic neural clusters (6, 9, 26 and 31 in Fig. 2B,D), increasing the number of predicted Clytia neuropeptides to twenty-one (Table S3). Our pseudotime ranking revealed that many of these predicted neuropeptides mark the later stages of neural cluster trajectories, likely defining distinct, mature neural subpopulations (Fig. 4D).

We extracted and re-clustered the neural supergroup ('Neural', Fig. 2A,B) to characterize neural subtypes. This distinguished 14 subpopulations of neurons as well as a progenitor population, expressing cell cycle and conserved neurodevelopmental genes including the bHLH transcription factor Neurogenin (sub-cluster 0, Fig. 4D). Strikingly, the neuronal subpopulations show combinatorial neuropeptide precursor expression, often with a distinct and identifying neuropeptide (Fig. 4D). Expression of putative neuropeptide processing enzymes was detected across the subpopulations, as well as in the nematocytes and gland cells (Fig. S17). The UMAP expression embedding for neuropeptide precursors suggest that further complexity remains to be discovered within these subpopulations, for instance, subpopulation 4 includes cells either expressing pp1 or pp11. In situ hybridization for a set of neuropeptide precursors indicated that some were broadly distributed across the animal while others had highly specific spatial locations, suggestive of distinct functions in regulating physiology and behaviors relating to swimming, feeding, and orientation (Fig. 4E; Fig. S18). For example, pp5+ (GRFamide precursor) cells were widely detected across the tentacles, bulbs, nerve rings, subumbrella, and mouth, while pp11+ (GLWamide precursor) cells were detected predominantly in the manubrium and nerve ring. Pp7+ cells were located both around the rim of the mouth and in patches of the nerve ring, close to the statocysts (vestibular organs). Subpopulations of neurons also occupied different expression domains within the tentacles: pp25, which generates distinct RFamide family neuropeptides (Table S3; Fig. S18), distinguishes a subpopulation of the pp5+ neurons positioned on the aboral side of the tentacle, pp17 labels distinct cells in the same region, and pp20 labels a small group of cells at the base of each tentacle (Fig. 4E,F; Fig. S18).

Unlike cells using neuropeptides, the transcripts for which are directly assayable in scRNA-seq data, cells using classical chemical neurotransmitters are identified by enzymatic or transporter proxies. Cells using glutamate as a neurotransmitter are usually inferred via the presence of vesicular glutamate transporter markers (vGluts -in human SLC17A6/7/8) [START_REF] Serrano-Saiz | Modular control of glutamatergic neuronal identity in C. elegans by distinct homeodomain proteins[END_REF]. Of the closest Clytia homologs of human vGluts, one was expressed in neurons, as well as nonneural cell types, and several were marker genes for nematocytes. However, in common with nearly all cnidarian genes annotated in silico as vesicular glutamate transporters, Clytia sequences lack an arginine residue, conserved in all bilaterian vGluts, and recently shown to be required for vGlut function [START_REF] Li | Ion transport and regulation in a synaptic vesicle glutamate transporter[END_REF]. Further, we detected no neuronal expression of the closest homologs of glutamate decarboxylase (GAD), a marker for GABAergic neurons in Bilateria. GAD homologs were instead detected in some gastrodermal subtypes. There were a number of other interesting neuron subtype-specific genes potentially involved in 'chemical' neurotransmission, for example, a possible nitric oxide synthase, a choline (SLC5A7-like) transporter, a gene encoding a taurine dioxygenase-like domain, and a member of the Slc6 family of transporters (Table S5).

Cell State Shifts in Response to Starvation Across the Cell Atlas

To assess the transcriptional impact of starvation, we mapped individual cells to their corresponding control or starved labels. As there are around 60% fewer cells in a starved animal (Fig. S2), we first asked whether there were significantly different numbers of cells per cluster between control and starved conditions. We found that only one cluster had a significant difference (cluster 11, early nematoblasts in Fig. S5B), suggesting a nearly uniform reduction across cell types in the starved condition. In contrast, the distribution of cells from control versus starved animals across the atlas embedding showed dramatic shifts in the local density of cells from the two conditions within most clusters (Fig. 5A).

Given the cell type resolution of the atlas, as determined operationally by clustering, we then asked how drastic the transcriptional changes incurred by perturbation were in comparison to the transcriptional differences defining the cell types i.e. are the perturbation-induced changes encompassed within these cell type designations or are they larger in magnitude. We thus compared distances between control and starved cells within clusters to the distances between clusters. As a metric we used the L1 distance (see Methods), the sum of the absolute differences between centroid coordinates in PCA-reduced space. We found that the L1 distances between control and starved cells within a cell type, versus between cell types (regardless of condition), formed nearly non-overlapping distributions (Fig. 5B) (see Methods). This suggests that, overall, in Clytia the transcriptional responses to starvation are defined by cell state shifts, and their cell type repertoire is well represented by the original clusters. However, the impact of starvation was variable across cell types, as reflected by the range of internal (state) distances (Fig. 5B). Starvation produced the largest perturbations in cells of the gastrovascular system, causing control-vs-starved distances large enough to overlap with the smallest inter-type distance, i.e. that between the stem cells and nematocyte precursors (Fig. 5B). This distinction between state shifts and type was also clearly visible in the lack of overlap between the distributions of inter-and intra-cluster distances within the second multiplexed experiment (Fig. S6E). Though classification and distinction of cell state and type is a complex task [START_REF] Trapnell | Transcript assembly and abundance estimation from RNA-Seq reveals thousands of new transcripts and switching among isoforms[END_REF], this analysis, based on relative distance in transcriptional space, provides a quantitative basis for delineation of type/state effects that may be useful in other contexts. We additionally validated the ability of this method to recapitulate the magnitude of state shifts in response to graded stimuli as well as state versus type distinctions, on two other published, multi-perturbation datasets (Fig. S19) (see Methods).

To characterize gene-level responses underlying these starvation-induced shifts, we then asked if responses are shared or unique across the cell types and compared the extent of the responses, in terms of gene quantity and expression level, across the atlas. For each cell type, we collected genes that were differentially expressed under starvation ("perturbed genes" Fig. 6A) (see Methods). For a high-level view of the general functions and processes affected by starvation, and their cell type specificity, we clustered perturbed genes into apparent "gene modules" (56) by their pattern of co-expression across cells (Fig. 6A) (see Methods). We assigned putative functions to these gene modules through GO term enrichment, giving a global view of affected processes (Fig. S20), and examined the distribution of cell types across modules by asking in how many cell types is a given gene a perturbed gene (Fig. 6B). We found that certain gene modules were broadly shared across cell types, while others were almost entirely cell typespecific (Fig. S20). Striking examples include gene module 5, which is enriched in proteolytic genes (Fig. 6C; Fig. S20) and has shared expression across multiple GD cell types (Fig. 6C). Notably, there is also divergent gene expression between GD types (Fig. 6C). In comparison, gene module 3 is largely composed of early oocyte gene expression (~70%), and is enriched in cell cycle and developmental genes, which are commonly enriched in growing oocytes (Fig. 6D). Changes in expression of these genes likely reflect the processes of oocyte phagocytosis activated in the gonads of starving animals (see below). Overall, these modules thus provide an overview of which processes affected by starvation are shared across cell types and reveal divergent expression potentially reflecting different motility [START_REF] Sinigaglia | Pattern regulation in a regenerating jellyfish[END_REF], locations, or as-yet-undescribed functional differences between the cell types (Fig. S21).

To examine how individual perturbed genes are distributed across cell types, we visualized, for each cell type, how many perturbed genes it had, and how many of these genes are unique versus shared with other cell types (Fig. 6E). We found a large number of perturbed genes (~72%) were cell type specific (Fig. 6E). For the most perturbed cell types, we examined whether the state shifts that we had observed were due to changes in a large or small number of genes, and how highly these genes were expressed. Consistent with the marked shrinkage of the gonads during starvation treatment (Fig. 1), early oocytes contained the highest number of perturbed genes, which were spread across many gene modules (Fig. 6D). In contrast, the GD cell types had fewer perturbed genes that were expressed at higher levels and localized to more specific gene modules (Fig. 6C), highlighting the diverse logic employed by cell types under starvation (Fig. S22).

In accordance with these distinct responses in GD cells and oocytes, comparison of the cellular organization of gonads from control and starved medusae revealed major reorganization of both the gastrodermis and the oocyte populations (Fig. 7). Most strikingly, the population of mid-sized, growing oocytes, which progress daily through vitellogenesis in conditions of normal feeding [START_REF] Amiel | Three distinct RNA localization mechanisms contribute to oocyte polarity establishment in the cnidarian Clytia hemisphaerica[END_REF], was largely depleted following starvation, leaving a majority of pre-vitellogenic oocytes (Fig. 7A). A sparse population of large oocytes in starved gonads likely results from growth of a minor subpopulation of oocytes fueled by recycling of somatic tissue and oocytes (disintegration and phagocytosis of smaller oocytes visible in Fig. 7A, asterisks). Consistently, GD cells in many parts of the gonad lost their regular epithelial organization and, despite the absence of any external food supply, showed evidence of active phagocytosis involving variablysized vesicles (arrows in Fig. 7A). Changes in organization and activity of the gonad gastrodermis were also evident from in situ hybridization images for the GD cell marker CathepsinL, while reduced expression was confirmed for a protease (ShKT-TrypA) expressed in gland cell types A and B positioned within the manubrium gastroderm, which is down-regulated during the starvation treatment (Fig. 7B). Shifts between gonad gastrodermis organization and transcriptional profiles induced by starvation thus accompany activation of tissue autodigestion programs, and likely also the mobilization of GD cells (termed MGD for Mobilizing Gastro Digestive cells [START_REF] Sinigaglia | Pattern regulation in a regenerating jellyfish[END_REF]) from the gonad through the gastrovascular canal system, which has been observed both under conditions of starvation and during regeneration of the feeding organ [START_REF] Sinigaglia | Pattern regulation in a regenerating jellyfish[END_REF].

Discussion

The Clytia medusa single-cell atlas presented here is an important addition to the growing number of single-cell atlases across the animal tree of life. It is available for easy interrogation by the community through the UCSC Genome Browser (see Data Availability). This provides the first cell-level transcriptomic characterization of a pelagic medusa stage, the most complex of the life cycle forms within the large and diverse phylum Cnidaria. Reflecting this complexity, we found greater cell type diversity in the Clytia medusa than in its polyp-only hydrozoan cousin Hydra [START_REF] Siebert | Stem cell differentiation trajectories in Hydra resolved at single-cell resolution[END_REF]. The outer, epidermal body layer could be sub-divided into seven clusters encompassing all of the described Clytia muscle types, including two types of fast-contracting striated swimming muscle [START_REF] Steinmetz | Independent evolution of striated muscles in cnidarians and bilaterians[END_REF][START_REF] Leclère | Diversity of Cnidarian Muscles: Function, Anatomy, Development and Regeneration[END_REF]. Rich diversity was also uncovered in the inner gastroderm layer, which is elaborated in the medusa into distinct digestive compartments (mouth, stomach, gonad, and tentacle bulb) and also generates the thick mesoglea (jelly) characteristic of the medusa form. Of the eight gastroderm cell clusters, four could be mapped to distinct sites by marker gene in situ hybridization including three likely involved in mesoglea/ECM production and thus in modulating the medusa structure. In contrast, the six 'GastroDigestive' clusters, characterized by a largely shared transcriptomic profile (Fig. 2B; Table S5), are distributed more uniformly across the digestive compartments (Fig. S14), suggesting a dominant role in intracellular digestion but with functional specializations that remain to be fully understood. Our starvation experiment analyses revealed that these clusters were maintained operationally as distinct 'cell types' rather than 'cell states' between the two extreme conditions tested, but we cannot rule out that responses to other environmental or physiological perturbations may reveal plasticity between these clusters, for instance transdifferentiation between muscle and nerve cell types is well documented in hydrozoan medusae (overview in [START_REF] Leclère | Diversity of Cnidarian Muscles: Function, Anatomy, Development and Regeneration[END_REF]).

In addition to the epithelial cell types of the epidermis and gastroderm, our single-cell atlas confirms the presence of an interstitial stem cell (i-cell) population in Clytia providing a similar set of somatic cell types to that described in Hydra, as well as the germ cells [START_REF] Siebert | Stem cell differentiation trajectories in Hydra resolved at single-cell resolution[END_REF][START_REF] Leclère | Maternally localized germ plasm mRNAs and germ cell/stem cell formation in the cnidarian Clytia[END_REF]. In this medusa data we do not find strong evidence for direct progression from i-cells to gland cells, or for the shared neural-gland cell progenitors described in [START_REF] Siebert | Stem cell differentiation trajectories in Hydra resolved at single-cell resolution[END_REF]. In contrast, our pseudotime analyses provide transcriptional signatures of the progressive stages of nematogenesis and neurogenesis from i-cells that will guide future studies of their developmental regulation. The large representation of nematogenic stages in this Clytia medusa scRNA dataset allowed us to link two distinct phases of nematocyte formation with extremely different transcriptional profiles. The initial phase covering nematocyst formation has been the focus of many studies [START_REF] Condamine | Molecular characterisation of a cellular conveyor belt in Clytia medusae[END_REF][START_REF] Denker | Ordered progression of nematogenesis from stem cells through differentiation stages in the tentacle bulb of Clytia hemisphaerica (Hydrozoa, Cnidaria)[END_REF][START_REF] Sunagar | Cell type-specific expression profiling unravels the development and evolution of stinging cells in sea anemone[END_REF][START_REF] Denker | Horizontal gene transfer and the evolution of cnidarian stinging cells[END_REF], but the terminal phase has been largely overlooked in previous transcriptomics studies, likely due to the relatively low mRNA content [START_REF] Siebert | Stem cell differentiation trajectories in Hydra resolved at single-cell resolution[END_REF][START_REF] Condamine | Molecular characterisation of a cellular conveyor belt in Clytia medusae[END_REF] and the extremely abrupt degradation of nematocystrelated mRNAs before the terminal phase [START_REF] Sunagar | Cell type-specific expression profiling unravels the development and evolution of stinging cells in sea anemone[END_REF]. We uncovered 14 mature neuronal subtypes in Clytia, which is similar to the number reported in Hydra and Nematostella [START_REF] Sebé-Pedrós | Cnidarian Cell Type Diversity and Regulation Revealed by Whole-Organism Single-Cell RNA-Seq[END_REF][START_REF] Siebert | Stem cell differentiation trajectories in Hydra resolved at single-cell resolution[END_REF]. It is likely that further heterogeneity exists within these 14 subpopulations. Spatial expression analysis of neuropeptides that contributed to the signatures of one or more revealed a wide variety of neuronal populations either associated with specific anatomical structures, such as the tentacles, nerve rings, and manubrium, or distributed across the medusa (Fig. 4D,E; Fig. S18). How molecular cell type maps to function both within and across body parts, the roles of these peptides as primary transmitters and/or neuromodulators, and the uses, if any, of classical, small-molecule neurotransmission, remain unknown. Moving forward, with this cell atlas as the foundation, the ability to perform whole-organism, multiplexed scRNA-seq, in combination with emerging genetic tools and advantageous life history traits, makes Clytia a powerful, tractable platform for high-resolution systems biology. This work further serves as a case study in using multiplexed single-cell transcriptomics to assess cellular responses to whole organism perturbations, and provides a guide for deployment in other organisms. We anticipate WHAM-seq will benefit researchers studying various biologies, from developing embryos to organoids to non-model organisms. The techniques for multiplexed experimentation that underlie this study are also well suited to large-scale perturbation studies (such as temperature, pH, or other environmental disturbances) in other marine organisms given the sea-water compatible workflow. Though the inclusion of multiple animals and conditions may currently limit the detection of very rare cell populations (Fig. S1), as sequencing costs drop and cell throughput in scRNA-seq grows, WHAM-seq should become tractable for larger, more complex systems. The lack of library-induced batch effects demonstrates how large-scale experiments can be conducted without introduction (or minimizing introduction) of confounding factors from multiple experiments, which can be highly non-linear and difficult to account for [START_REF] Tran | A benchmark of batch-effect correction methods for single-cell RNA sequencing data[END_REF]. The second perturbation dataset also demonstrates both how batch effects are reduced within multiplexed experiments (Fig. S11) (see Methods) and the utility of this multiplexed approach in discerning cell type-specific activity markers (Fig. S11) (see Methods).

The fully reproducible framework we have presented, which includes code that can be run on a laptop or for free in the cloud (see Code Availability), will further assist in extending this expression-based analysis to other organisms. By relying on expression, our strategy reduces the reliance on prior gene functional annotation, using specificity of expression to identify genes of interest, allowing for targeted annotation. This includes determination of strong diagnostic markers for cell type definition, cell type specific and shared transcriptional responses to starvation, and 'modules' of co-expressed genes underlying these responses. The extent of these expression-based changes additionally highlights areas of the organism's biology that are strongly or uniquely affected by a perturbation. By applying simple and interpretable quantitative analyses to the various cell-type specific perturbation responses, we revealed the large-scale downregulation of gene expression in two GastroDigestive cell types and severe disruption of oocyte development under starvation. Together, this approach dramatically lowers the barriers for working with non-traditional models, and affords opportunities to match uniquely-suited organisms to specific questions. Moving forward, the combination of scRNA-seq and other sequencing-based genomics techniques with multiplexing and annotation-agnostic analyses, could foster comprehensive high-resolution molecular studies of diverse organisms and their responses to numerous environmental perturbations.

Materials and Methods

Animal culture and experimental setup Starvation

Culture of the Clytia life cycle was carried out as previously described [START_REF] Lechable | An improved whole life cycle culture protocol for the hydrozoan genetic model Clytia hemisphaerica[END_REF], with some modifications to the tank design. The system used in this study to culture polyps uses zebrafish tanks (Pentair), with polyp slides held in glass slide racks (Fisher, cat#02-912-615). Medusae used in this experiment were raised in 4L beakers with a circular current generated by stirring with a constant speed 5rpm DC motor (Uxcell), attached to the lid of a multiwell tissue culture plate. Artificial sea water for culture and experiments was made using Red Sea Salts (Bulk Reef Supply, cat# 207077) diluted into building deionized water to 36ppt. Experiments used ~1cm female medusa of the Z4B strain.

For the experiment, baby medusae were collected overnight and then cultured together until they reached ~1cm (about 2 weeks). Animals were fed once per day using 2-4 day old brine shrimp. Prior to the experiment, animals were split into two beakers. Feeding continued as before for one beaker, while the other was starved for 4 days. The "control" group was not fed on the day of the experiment. The 4 day time point was chosen as animals show strong phenotypic changes (Fig. 1) but it is far from their survival limit following starvation, as Clytia medusae can survive for more than 3 weeks with no food.

Stimulation

Rearing of the medusae prior to the experiment was performed as described above. On the day before the experiment, each Clytia medusa (3-5 weeks old) was placed in a separate container (~150 mL seawater), which was covered with foil and moved to the experimentation area to acclimate overnight in the dark. The morning following overnight acclimation, lights were turned on for ~2.5 hours to allow for spawning to complete [START_REF] Quiroga Artigas | A gonad-expressed opsin mediates light-induced spawning in the jellyfish Clytia[END_REF]. Each animal was then given repeated bouts of stimulation over a period of 30 mins, with each stimulus administered every 2 minutes. 100 uL of each stimulant (150 mM KCl, DI water, or seawater (SW) as a control) was gently added just below (or just above for KCl) each medusa by pipette (Fig. S7A). Stimuli were chosen based on their ability to reliably induce crumpling behavior, a protective response in which the bell is drawn in towards the mouth using the radial muscle [START_REF] Hyman | OBSERVATIONS AND EXPERIMENTS ON THE PHYSIOLOGY OF MEDUSAE[END_REF]. Animals were dissociated 30 mins following the last stimulation.

Single-cell suspension and multiplexing

For the starvation experiment, animals were washed with hypertonic PBS (500 mM NaCl, 2.7 mM KCl, 8 mM Na2HPO4, and 2 mM KH2PO4, pH 7.4) by serial transfer from seawater through three successive containers each with 150 mL hypertonic PBS, to prepare cells for fixation and to avoid precipitation of seawater salts in methanol. The animals were then thoroughly homogenized with a dounce. After homogenization, cells were collected by centrifugation at 500 x g for 5 minutes and resuspended in 100 uL hypertonic PBS. Cells were then fixed by addition of 400 uL ice-cold methanol and stored at -80C until sample indexing and library preparation. Each sample was labeled according to the ClickTag labeling procedure described previously [START_REF] Gehring | Highly multiplexed single-cell RNA-seq by DNA oligonucleotide tagging of cellular proteins[END_REF]. ClickTags used for each animal are outlined in Table S3. Each sample was labeled with two distinct and unique ClickTags, and samples were then pooled after addition of the 'blocking' oligo. Control and starved samples were counted on a Countess, both to estimate cell numbers per animal and to determine concentrations for 10X loadings. 200,000 cells/mL were counted for the starved sample and 1 M cells/mL for the control. Around 120-150,000 cells of the starved and 200,000 cells of the control were then pooled. A loading of 20,000 pooled cells was then used as input into two lanes of the 10X Chromium Controller with v2 chemistry. Sample tag libraries were separated and processed after an SPRI size-selection step as previously described [START_REF] Gehring | Highly multiplexed single-cell RNA-seq by DNA oligonucleotide tagging of cellular proteins[END_REF]. cDNA samples were run on 2 lanes of HiSeq 4000 (two HiSeq 3000/4000 SBS 300 cycle kits) and tag libraries were run on 2 lanes of MiSeq (using MiSeq v3 150 cycle kits).

For cell counting (Fig. S2), the starvation experiment was repeated with 4 animals per condition with animals dissociated and resuspended in 500 uL seawater, and then cells counted on a hemocytometer (InCyto DHC-B02) with a 10X objective. Two 16 square grids were counted per sample. The same protocol was followed for fixing and labeling cells from the stimulation experiment. Each animal in this case was assigned one unique ClickTag and one ClickTag per condition (Table S3). After the separation of the cDNA and ClickTag samples, ClickTags were added at a 3% final concentration to the cDNA samples sequenced on the HiSeq, in addition to separate sequencing on the MiSeq. The full protocol is described in the Supplementary Methods.

In Situ Hybridization

Colorimetric in situ hybridization (ISH, Fig. 2,6b) was performed as previously [START_REF] Sinigaglia | A safer, urea-based in situ hybridization method improves detection of gene expression in diverse animal species[END_REF] with minor modifications. In brief, 2-week-old medusae (Z4B strain) were relaxed in 0.4 mM menthol in sea water and tentacles were trimmed prior to fixation in a pre-chilled solution of 3.7% Formaldehyde, 0.2% Glutaraldehyde in PBS on ice for 40 min. Specimens were then washed thoroughly with PBST (PBS + 0.1% Tween20), dehydrated in methanol stepwise and stored in 100% methanol at -20°C. Hybridisation (at 62°C for 72h ) and washing steps were performed in a robot (Intavis AG, Bioanalytical Instruments) using 20X SSC pH adjusted to 4.7 throughout. Acetylation steps using 0.1M triethanolamine in PBST (2x5min) then 0.25% acetic anhydride in 0.1M triethanolamine (2x5min) followed by PBST washes (3x10min) were included before prehybridization to reduce probe non-specific binding. Incubation with Anti-DIG AP, 1:2000 in 1X blocking solution was performed for 3h before washing and the NCB-BCIP colour reaction at pH 9.5. Following post-fixation, washing and equilibration of samples in 50% glycerol/PBS, images were acquired using a Zeiss AxioImager A2.

Probes were generated by PCR from cDNA clones corresponding to our EST collection [START_REF] Chevalier | Polarised expression of FoxB and FoxQ2 genes during development of the hydrozoan Clytia hemisphaerica[END_REF] or from medusa cDNA; the Elav probe was synthesized as a gBlock by Integrated DNA Technologies (details in Table S3). For probes against the Elav, the T3 polymerase recognition site (AATTAACCCTCACTAAAGGG) was added to the 3'-end of the PCR product, or gBlock, respectively. Products were TOPO cloned (ThermoFisher cat#K280020) and sequence verified. All probes were labelled with DIG RNA labelling mix (Sigma-Aldrich 11277073910), and purified with ProbeQuant G-50 Micro Columns (GE HealthCare Life Sciences Cat# 28-9034-08).

Initial Cell Ranger Demultiplexing and Clustering for cDNA

Initial processing of starvation cDNA libraries was performed with the 10X Cell Ranger pipeline, using Cell Ranger 3.0 count and aggr functions to align and quantify the HiSeq reads, and combining the counts from the two lanes using the denoted sample IDs. This was followed by filtering cells for the high-quality cells chosen during ClickTag analysis, in addition to filtering cells by thresholding the rank-UMI vs cell barcodes plot. Values were log1p normalized, mean-centered, and scaled for downstream dimensionality reduction and visualization using Scanpy [START_REF] Wolf | SCANPY: large-scale single-cell gene expression data analysis[END_REF]. [Code] We then conducted Louvain clustering [START_REF] De Meo | 11th International Conference on Intelligent Systems Design and Applications (ieeexplore.ieee.org[END_REF] on the data mapped to a lower dimensional space, by applying PCA to the expression data filtered for highly variable genes, initially using Scanpy's filter_genes_dispersion on only the log-normalized data. This resulted in the identification of 36 clusters (Fig. 2B; Fig. S9-S11), which we also refer to as cell types. The marker genes were selected by analyzing the top 100 markers extracted by Scanpy's rank_genes_groups using default settings (p-values adjusted with the Benjamini-Hochberg method for multiple testing). The clusters were annotated and validated with marker genes previously identified in the literature, and manually categorized into the seven classes in Fig. 2A based on the marker gene patterns and their functional annotations. [Code]

kallisto bustools for Demultiplexing and Clustering: Standardization of Workflow

In order to integrate and update the analysis using a platform with streamlined ClickTag demultiplexing and count matrix generation workflows, we used the kITE demultiplexing protocol (link), which is based on the kallisto-bustools workflow and is described in the ClickTag demultiplexing protocol [START_REF] Gehring | Highly multiplexed single-cell RNA-seq by DNA oligonucleotide tagging of cellular proteins[END_REF]. Briefly, MiSeq reads are aligned to possible tag sequences (Hamming distance 1 away from designed oligo sequence whitelist) by building a kallisto index and pseudoaligning reads to this index. Counts for these sequences were then collapsed into counts for their respective ClickTags, creating a cell-by-tag count matrix. We used Louvain clustering of cell barcodes based on the observed ClickTags, to filter for clearly delineated cells (clusters strongly marked by the individual's two corresponding tags) and exclude sample doublets [START_REF] Gehring | Highly multiplexed single-cell RNA-seq by DNA oligonucleotide tagging of cellular proteins[END_REF]. We also followed similar pre-processing to standard cell-by-gene workflows, using the inflection point in rank-UMI vs cell barcodes knee plots to filter cell barcodes based on their tag UMI counts. We do find, similar to findings in the original ClickTag methodology [START_REF] Gehring | Highly multiplexed single-cell RNA-seq by DNA oligonucleotide tagging of cellular proteins[END_REF], that the number of ClickTags per fed cell is higher than for the starved counterparts, possibly following the previous observation of ClickTag number reflecting cell size (Fig. S4E

, shown for oocytes). [Code]

For the stimulation experiment we concatenated sequencing data from the MiSeq and HiSeq as input to the previously described kallisto-kITE workflow. With the same clustering procedure we selected cell barcodes in clusters with strong overlapping expression of both individual and condition ClickTags. [Code] To standardize the cDNA analysis workflow, we re-processed the starvation data using the kallisto-bustools workflow to generate gene count matrices for each lane, which were then concatenated. Cells were also filtered based on the ClickTag analysis. Values were log1p normalized, mean-centered and scaled, and filtered for highly variable genes using the same procedure described above for downstream dimensionality reduction and visualization (e.g. PCA) using Scanpy. We found that with the kallisto-processed data (with the same Cell Ranger clustering applied to the cells), the top 100 markers for each of the 36 clusters determined with Scanpy's rank_genes_groups function (Table S5) (using the nonparametric Wilcoxon test) overlap with markers in the Cell Ranger expression data, verifying that the cluster labels were concordant (Fig. S10). To ensure that we additionally detect low expression marker genes we extracted low expression markers in Table S5, which consist of genes with at least 10 counts over all cells, and with 90% of those counts deriving from the same cell type. We used the PAGA graph abstraction method [START_REF] Wolf | PAGA: graph abstraction reconciles clustering with trajectory inference through a topology preserving map of single cells[END_REF] to generate an underlying graph representation of the connectivity between cells (determining connectivity by the number of inter-edges between cell groups compared to the number of inter-edges under random assignment). We then generated a 2D UMAP (32) embedding initialized with the PAGA graph structure for cell atlas visualization (Fig. 2A) using Scanpy. For Fig. 2B, 100 cells were randomly subsampled from each cell type to generate the heatmap. [Code] The stimulation cDNA data was processed with the same kallisto-bustools workflow and commands as the starvation experiment data. [Code] We initially used Louvain clustering to also filter low UMI count clusters that were then removed from downstream analysis.

[Code], [Code]

Distance-based Comparative Analysis of Clusters [Code]

We first used L1 distances between starved and control cells (within each cell type) to assess how comprehensive our cell type designations were (Fig. 5). Centroids for a given cell type were calculated for starved and control cells separately, in PCA-reduced space (60PC coordinates for each cell as opposed to the raw gene expression matrix). The centroid vectors are represented as cs and cc for the Starved and Control cells respectively. The L1 distance (d) between them was calculated as the sum of absolute difference between the centroid coordinates:

𝑑 = ∑ $% &'% | 𝑐 + , -𝑐 . , |. (1)
These intra-cluster distances were then compared to the pairwise L1 distances between the cell types, with centroids calculated for all cells in a given type (c1 and c2 for cell types 1 and 2) in the same manner as described above. The L1 distances were then calculated for all possible pairs of these cell type centroids using [START_REF] Hwang | Single-cell RNA sequencing technologies and bioinformatics pipelines[END_REF]. The distributions of the inter-and intra-cluster distances are shown in (Fig. 5B). We chose to use the L1 distance metric as it tends to better retain relative distances in high dimensions, particularly in comparison to the commonly used Euclidean distance or other higher L-norms [START_REF] Ntranos | Fast and accurate single-cell RNA-seq analysis by clustering of transcript-compatibility counts[END_REF][START_REF] Aggarwal | On the Surprising Behavior of Distance Metrics in High Dimensional Space[END_REF].

We validated this method of transcriptional distance measurement on the multiplexed perturbation of mouse NSCs in [START_REF] Gehring | Highly multiplexed single-cell RNA-seq by DNA oligonucleotide tagging of cellular proteins[END_REF] and the multiple immunomodulatory drug treatments across heterogeneous cell populations in human PBMCs in [START_REF] Chen | Dissecting heterogeneous cell populations across drug and disease conditions with PopAlign[END_REF]. For [START_REF] Gehring | Highly multiplexed single-cell RNA-seq by DNA oligonucleotide tagging of cellular proteins[END_REF], we calculated centroid distances between perturbed and control populations, in PCA space [START_REF] Gehring | Highly multiplexed single-cell RNA-seq by DNA oligonucleotide tagging of cellular proteins[END_REF], as well as pairwise distances between individual perturbed and control cells in [START_REF] Gehring | Highly multiplexed single-cell RNA-seq by DNA oligonucleotide tagging of cellular proteins[END_REF] to demonstrate the ability of the L1 distance to recapitulate the graded response to a perturbant, with larger distances representing greater perturbant impact (with respect to the control cells) (Fig. 19A,B). We additionally compared all pairwise L1 distances between perturbed and control cells within the Monocyte population or T-cell PBMC populations in (71) using the oNMF space. We then compared these measurements against the magnitude of cell type distances (between control Monocyte and T-cell populations) to highlight 'state'-versus 'type'-level transcriptional differences (Fig. 19C). To then extract the most perturbed Monocyte populations we calculated centroid distances between the centroid of each perturbation condition and the CTRL1 (control) condition to rank populations by perturbation 'distance' (Fig. 19

D). [Code Gehring et al.] [Code Chen at al.]

We then used the stimulation dataset to assess the validity of the clusters/cell types generated with the starvation data. We created a joint representation of the two datasets by using a concatenated cell-by-gene matrix including only the genes highly variable in both the control and starved datasets, and used 70% of the starvation data to train a k-nearest neighbor (KNN) classifier (using sklearn's KNNClassifier with k = 15) to assign cluster labels to the remaining starvation cells and to the stimulation dataset. This showed that the stimulation cells' labels from their neighbors in starvation data were assigned at the same accuracy as the test starvation data, meaning that clusters from the starvation data are applicable to the stimulation dataset and capture the main features of the stimulation dataset to the same extent (Fig. S6A).

We also examined batch effects in the stimulation experiment. As with the starvation experiment, the L1 metric was utilized to visualize the magnitude of batch effect within the multiplexed experiments compared to between experiments (Fig. S8). We used the merged representation (used for the KNN assignment) between both experiments to find the average pairwise distances between cell types of control condition individuals within the starvation experiment, within the stimulation experiment, and across both experiments. We found that cell type distances between organisms were reduced within multiplexed experiments compared to distances across experiments.

The merged atlas was also used for determination of strong in situ markers, mainly for gastrodermal cell types, since these highly related types share many marker genes.

RNA-seq Analysis & Clustering with MARIMBA Annotation

To make our dataset more easily searchable with the MARIMBA v.1 (created on 05/30/2016) transcriptome annotation [START_REF] Leclère | The genome of the jellyfish Clytia hemisphaerica and the evolution of the cnidarian life-cycle[END_REF], we generated single-cell gene count matrices with respect to transcript sequences distributed via the MARIMBA website (link). The gene count matrices were generated using the same kallisto-bustools workflow, and compared the application of the previous clustering/cell type assignment (including the overlap in differentially expressed genes delineating these clusters) to validate the quantification derived from the Trinity/Cuffcompare assembled transcriptome (Fig. S13). [Code] We also produced a notebook for visualization of gene expression in this dataset to facilitate its use in future studies. This establishes a code base for rapidly and transparently processing and comparing single-cell datasets with future transcriptome annotations. [Code] Neural Analysis [Code] We clustered all cells within the broad class of 'Neural' with Louvain clustering to obtain distinct subpopulations (labeled in Fig. 2A). Markers were determined with Scanpy's rank_genes_groups function (using the Wilcoxon test) for each subpopulation (Table S3).

Marker genes from neural clusters, with predicted signal peptides [START_REF] Almagro Armenteros | SignalP 5.0 improves signal peptide predictions using deep neural networks[END_REF], were screened for candidate neuropeptide cleavages sites (regular expression G[KR][KRED]). In cases where a sequence had more than one match to this motif, the 6 residues immediately N-terminal to the motif were inspected for similarity to each other -when similarity was present, the protein was considered a neuropeptide candidate. Predicted sequences can be found in Table S3.

Pseudotime Analysis [Code]

We selected cells from cell types of interest (i-cells, nematoblasts, nematocytes, and neural cells) and used diffusion mapping [START_REF] Haghverdi | Diffusion pseudotime robustly reconstructs lineage branching[END_REF] to create a reduced dimension representation of cells, along with Scanpy's dpt function which uses geodesic distance along the graph of cells (in the determined 'diffusion component' space) to estimate pseudotime. We then computed a PAGAbased embedding to visualize the cells in the context of the different trajectories with a ForceAtlas2 layout. To determine which genes constituted the important features in a given pseudotime trajectory, we implemented a method based on the random forest method used in the dynverse R package for extracting 'important' genes [START_REF] Saelens | A comparison of single-cell trajectory inference methods[END_REF]. A random forest regression model implemented with sklearn's random_forest_regressor was used to identify genes that were good predictors of the generated pseudotime values (grouped into quantiles). This was run for each of the two inferred trajectories separately (stem cells to nematocytes, and stem cells to neurons). The training set consisted of 80% of the genes' expression data, within which 80% was used for optimizing the model, and 20% was used to evaluate the mean squared error (MSE). Both trajectory models had an 𝑅 1 of 0.85 or greater. The remaining 20% of data from the full dataset was used to calculate the gene-wise permutation importance scores, providing a ranking of each feature (gene) in terms of its contribution to the model's predictive capabilities using sklearn's permutation_importance. In the ranking, positive scores indicate importance, with the importance 𝑖 3 for gene j being the difference between the original score s (MSE) calculated from the original model, and the average score across K=5 random permutations of the feature columns in the validation dataset: Genes with non-zero (positive) scores in the permutation test were retained and ranked (Table S5).

𝑖 3 = 𝑠 -
To assess the relationship of the gland cells to the neural populations, in investigation of a neural-gland progenitor population, we applied the URD tree construction pipeline [START_REF] Farrell | Single-cell reconstruction of developmental trajectories during zebrafish embryogenesis[END_REF] to create a pseudotime-biased dendrogram of the neural and gland cell types beginning at the defined i-cell population. This included the i-cell, gland cell, and neuron clusters. We additionally plotted expression along the generated tree of Myb and Myc3 orthologs corresponding to the neural-gland progenitor markers in [START_REF] Siebert | Stem cell differentiation trajectories in Hydra resolved at single-cell resolution[END_REF] as well as other progenitor/developmental markers from [START_REF] Siebert | Stem cell differentiation trajectories in Hydra resolved at single-cell resolution[END_REF] for the gland and neural cell types (Fig. S16).

Perturbation Response Analysis Extracting DE (perturbed) genes

We used a likelihood ratio test (LRT) with negative-binomial based models of each gene's expression in DESeq2 [START_REF] Love | Differential analysis of count data--the DESeq2 package[END_REF][START_REF] Wang | Comparative analysis of differential gene expression analysis tools for single-cell RNA sequencing data[END_REF], with the reduced model not including the condition label (control or starved) and single-cells treated as individual replicates. There are other approaches that can be used [START_REF] Crowell | muscat detects subpopulation-specific state transitions from multi-sample multi-condition single-cell transcriptomics data[END_REF][START_REF] Luecken | Current best practices in single-cell RNA-seq analysis: a tutorial[END_REF]; we found that this method yielded many biologically relevant genes, which we could validate via the in situ experiments. All nonzero-expression genes were used for analysis. From the LRT we obtained p-values (corrected with the Benjamini-Hochberg multiple testing correction method across the number of genes tested) to determine if the gene's expression was significantly affected by the condition, and thus the perturbation. Genes with alpha < 0.05 and |log2FC| > 1 were selected as significant. We used the parameters sfType="poscounts", minmu=1e-6, minReplicatesForReplace=Inf in the DESeq2 model. Clusters with greater than 100 cells in each condition were subsampled evenly from fed and starved cells, to reduce the effect of uneven cluster sizes on differential expression analysis (Fig. S5). Clusters with less than 10 cells in any condition were not used for this analysis. To create the UpSet Plot (80) in Fig. 2B, p-values from the LRT analysis were corrected for multiple testing (Bonferroni correction with n = number of clusters, since the test is assessing the intersection of genes across all cell types). All genes are included in Table S5.

[Code for Starvation Responses] [Code for Stimulation Responses]

We additionally include low expression DE genes detected by pseudo-bulk analysis in Table S5. We applied the same methodology as above, treating each animal as a replicate (5 control and 5 starved replicates) and summing the counts per gene for each organism. For the DESeq2 model we used the default parameters, as we were simulating bulk results.

[Code] qPCR Validation for Upregulated, Stimulation-Responsive Genes

To further validate the gene candidates from the stimulation multiplexed experiment, which displayed immediate early-gene like expression under KCl or DI stimulation, we chose broadly expressed candidates i.e. DE genes in multiple cell types, for whole animal, bulk qPCR. Here we looked for upregulation of these genes as compared to a 'housekeeping' collagen marker, XLOC_008048. Two organisms were housed together in 150ml of seawater, with four animals total per each condition (SW, KCl, DI) treated as described previously in the Methods. The two animals were homogenized together in 1ml of Trizol on ice with a syringe. 0.2ml of chloroform were added per 1ml Trizol and shaken for 15s. After spinning at 12000xg for 15min at 4C the top phase is removed and in a new tube mixed with an equal volume of 60% EtOH. The RNEasy Mini kit (Qiagen Catalog # 74004) was then used for bulk rna extraction from the sample. Probes for the candidates and housekeeping gene were produced as previously described, and the LightCycler® 480 SYBR Green I Master Mix (2X) was used for 20ul qPCR reactions with 5ul of RNA template per reaction, using the standard cycling procedure with a 62C annealing temperature. For each rna sample (n = 6), three qPCR replicate reactions were done, resulting in a total of six replicates for each perturbation condition. Fold-change (log2FC) values presented in Fig. S7 were calculated by first normalizing Ct values for each replicate to the average housekeeping gene Ct value in the same condition (ΔCtNorm = CtRep -CtHouse_Avg). The log2FC was then calculated as log(2 ΔΔCtNorm ) where ΔΔCtNorm = ΔCtNorm -avg(ΔCtNorm_SW ), is the difference of each normalized Ct from the average normalized Ct for the SW (control) condition. [Code]

De novo perturbed gene clusters

To cluster genes affected by perturbation (Table S5) and to obtain information on their co-expression and possible functional similarity, we transposed the cell-by-gene expression matrix (obtaining a gene-by-cell matrix) for only the aggregated perturbed genes with padj > 0.05 (adjusted across genes and clusters) from the DeSeq2 analysis. We used Louvain clustering on the gene expression matrix, identifying both co-expressed genes and cell type specific genes, similar to Monocle's procedure for detecting genemodules [START_REF] Trapnell | The dynamics and regulators of cell fate decisions are revealed by pseudotemporal ordering of single cells[END_REF]. We then used the aggregated information from the 'modules' to determine putative functions/response types. The topGO weight algorithm (81) was used to determine gene ontology (GO) terms that were significantly enriched in each gene module compared to GO terms in all other groups, with significance threshold alpha < 0.05. The p-values were also adjusted for multiple testing over the number of different gene modules using Bonferroni correction, and only significant GO terms were used to label the response types among the modules (Fig. S20; Table S5). [Code] This analysis was replicated to determine gene modules distinguishing the GD cell types by clustering groups of co-expressed marker genes across the subtypes.

FiguresFig. 1 .Fig. 2 .

 12 Figures

Fig. 3 .

 3 Fig. 3. Developmental and Hierarchical Analysis of Nematocytes. (A) Subset of i-cell, neural, and nematocyte populations for trajectory analysis with respective clusters denoted. (B) Pseudotime trajectories for differentiating nematocytes. (C) Diagram indicating the two main phases of nematogenesis with the corresponding 'nematoblast' and 'nematocyte' clusters and representative genes. (D) Heatmap of expression of dynamic genes in nematocyte development across the maturing populations. * denotes published nematogenesis genes. (E) Trajectory embeddings for selected nematogenesis markers and corresponding in situs. Scale bar is 100µm. [Code]

Fig. 4 .

 4 Fig. 4. Developmental and Hierarchical Analysis of Neural Cell Types. (A) Pseudotime trajectories for differentiating neurons. (B) Heatmap of expression of dynamic genes in neural development over pseudotime. * denotes previously reported genes. (C) Embeddings show gene expression of markers for neural development across the neurons. (D) Diagram of the 14 neural subpopulations and progenitor cells along with a heat map of neuropeptide markers and the bHLH TFs Neurogenin and Hlh6 across these subpopulations. Hlh6, a member of the Atonal superfamily, corresponds to Nematostella Ath-like[START_REF] Richards | Regulation of Nematostella neural progenitors by SoxB, Notch and bHLH genes[END_REF]. + denotes previously unpublished neuropeptide candidates. XLOC identifiers and accession numbers in TableS3. (E) Subcluster embedding for neural subpopulation neuropeptide markers and corresponding in situs. Scale bar is 100µm. (F) Comparison of subpopulation marker expression between medusa territories.

Fig. 5 .Fig. 6 .

 56 Fig. 5. Starvation-Induced Shifts in Cell Types. (A) Density map of fed and starved cells across cell types, with each cell colored by the number of its 15 nearest neighbors that are control cells, to highlight the clusters most affected by starvation i.e. most clearly separated from their control neighbors. Clusters 14 and 19, with particularly dense and separated control/starved regions are circled. [Code] (B) Histogram of L1 distances between centroids of control and starved cells within cell types, versus pairwise L1 distances between centroids of the 36 cell types. Clusters 14 and 19, with largest internal distances are highlighted, and clusters with small pairwise distances also noted for comparison. Inset illustrates inter-and intra-L1 distance calculations. [Code]

Fig. 7 .

 7 Fig. 7. Perturbations to gastroderm cell types in response to starvation. (A) Confocal sections through gonads from control and starved medusae with cell morphology revealed by phalloidin staining of cell boundaries (magenta/grey). The first panel of each row shows co-staining of nuclei with Hoechst (blue) and endogenous GFP4 (green) in the outer epidermis (epi). Vitellogenic oocytes (vo) are largely absent during starvation, leaving a majority of previtellogenic oocytes (po). The gastroderm (g) is heavily reorganized, with evidence of active phagocytosis (vesicles arrowed) and disintegrating oocytes (asterisks). The third panel in each row is a higher magnification of the boxed area in the second panel, and the fourth panel shows a second example gonad for each condition. Scale bars all represent 50µm. (B) In situ hybridization visualization of example cell-types (left) and genes (right) affected by starvation. The general GD cell marker CathepsinL confirms extensive reorganization of the digestive gastroderm in starved animals, especially its collapse into central regions of the strongly reduced, oocyte-depleted, gonad, and fragmentation in the circular canal; Expression of the gland cell A/B marker "ShKT and trypsin domain protein A" (ShKT-TrypA) in the manubrium consistent with its strong downregulated following starvation. Gonads (go), manubrium (m), tentacle bulbs (tb), gastroderm (g), circular canal (cc). Scale bars all represent 200µm.

Acknowledgements

We thank Xiaolin Da and Xiao Wang for their technical assistance, Tsuyoshi Momose for his assistance with the single-cell experimentation, the Caltech Single-Cell Profiling and Engineering Center for use of their single-cell and sequencing tools, and the Caltech Bioinformatics Resource Center for transcriptome assembly and annotation analysis. We additionally thank the Caltech Center for Evolutionary Science for the bioinformatics resources to create a local UCSC genome browser. We thank A. Sina Booeshaghi for help with kallisto, bustools, the kITE demultiplexing of the ClickTag reads, and for rescuing the stimulation experiment sequencing data. We thank Jocelyn Malamy for helping to establish Clytia work at Caltech. We thank Sophie Peron for initial characterization of some of the cell type marker genes, Pascal Lapébie for identification of novel neuropeptide sequences, Muriel Jager for valuable advice on the in situ protocol and Julia Ramon Mateu for kindly providing pp11 probe.

Funding

JG, MH and LP were supported in part by a seed grant from the Chen Institute at the California Institute of Technology.

TC, JG and LP were supported in part by NIH U19MH114830 and NIH RF1AG062324A.

We thank the Marine Resources Centre (CRBM and PIV imaging platform) of Institut de la Mer de Villefranche (IMEV), supported by EMBRC-France. The French state funds of EMBRC-France are managed by the ANR within the investments of the Future program. LL was supported by the Agence Nationale de la Recherche (ANR-19-CE13-0003) AF, RRC, and EH, were supported by the H2020 / Marie Skłodowska-Curie ITN "EvoCell" Grant agreement N° 766053.

BW was supported in part by a Howard Hughes Medical Institute Fellowship of the Life Sciences Research Foundation and in part by NIH K99NS119749. This work was in part supported by the Whitman Center of the Marine Biological Laboratory in Woods Hole, MA, and a visiting grant from EMBRC-France.

DJA is an Investigator of the Howard Hughes Medical Institute.

Data Availability

All raw sequencing and processed data files used for analysis are available from CaltechData with links provided via the notebooks in the code repository. The sequencing read alignments as well as an interactive browser for gene expression visualization are publicly hosted on a UCSC Genome Browser at the Caltech Bioinformatics Resource Center.

Code Availability

All the code used to perform the analyses and generate the results and figures is available in Google Colab notebooks at https://github.com/pachterlab/CWGFLHGCCHAP_2021. The notebooks, which include the complete pre-processing of the raw data, provide a transparent implementation of the methods, and can be run for free in the Google cloud.

Confocal Microscopy

Visualization of cell morphology within the gonads of control and starved young adult female medusae (Z4B strain) by confocal microscopy was performed as previously [START_REF] Quiroga Artigas | A gonad-expressed opsin mediates light-induced spawning in the jellyfish Clytia[END_REF]. Fixation used 4% EM-grade paraformaldehyde in 0.1M Hepes pH 6.9/50mM EGTA/10mM MgSO4/80mM Maltose/0.2% Triton X-100 for 2h at room temperature. Specimens were washed 3x15min with PBS/ 0.02% Triton and 3x5 min in PBS. Cell boundaries and nuclei staining was performed by overnight incubation in 1/50 rhodamine-phalloidin (Molecular Probes 1mg/ml) and 1/5000 Hoechst 33258 (Img/ml stock, Sigma Aldrich) in PBS. Samples were washed 3x15min with PBS /0.02% Triton and 3x5min with PBS and equilibrated in 50% PBS/ Citifluor (Citifluor AF1) before imaging using a Leica SP5 confocal microscope. Control medusae were fixed 24h after the last feeding while starved ones were fixed 4 days after the last feeding.

Generation of Reference Transcriptome

Assignments to PANTHER database entries (version 11) were made using the 'pantherScore2.0.pl' script available from the database website (www.pantherdb.org) [START_REF] Mi | Protocol Update for large-scale genome and gene function analysis with the PANTHER classification system (v.14.0)[END_REF]. Human -Clytia orthologs were assigned using the OMA program (63) as described in [START_REF] Leclère | The genome of the jellyfish Clytia hemisphaerica and the evolution of the cnidarian life-cycle[END_REF], and taken from the pairwise human/Clytia orthologs output, rather than the orthologous groups.

To ensure highly sensitive transcriptome alignment, a new transcriptome assembly for Clytia hemisphaerica was generated from bulk RNA-seq data produced from Clytia medusae (organisms at the same life stage as in the single-cell experiments) (link). We used the Trinity (64) de novo assembler, with default parameters, to generate a transcriptome (link) and the Cufflinks Cuffcompare utility [START_REF] Trapnell | Transcript assembly and abundance estimation from RNA-Seq reveals thousands of new transcripts and switching among isoforms[END_REF] to merge the Trinity assembled transcripts with any XLOC annotations from the MARIMBA v.1 (created on 05/30/2016, link) transcriptome assembly [START_REF] Leclère | The genome of the jellyfish Clytia hemisphaerica and the evolution of the cnidarian life-cycle[END_REF]. Then with the CD-HIT (link) clustering method, assembled sequences with at least 95% were clustered and only one representative sequence was kept for each cluster. With the published MARIMBA v.1 genome sequence as a reference (created on 5/30/2016, link) [START_REF] Leclère | The genome of the jellyfish Clytia hemisphaerica and the evolution of the cnidarian life-cycle[END_REF], the GMAP (link) aligner converted the collapsed Trinity fasta records to gff3 coordinate file (link). Most differentially expressed genes found in single-cell RNA-seq data from this study were previously identified and annotated (Fig. S13). This annotation was used for the pre-processing and quantification of single-cell RNA-seq described below. Protein sequences were then obtained by running Transdecoder [START_REF] Haas | Others, TransDecoder (find coding regions within transcripts[END_REF] with default settings for the Trinity transcriptome (link).

Pre-Processing and Clustering of Sequencing Data Initial Cell Ranger Demultiplexing for ClickTags

Initial demultiplexing of ClickTag libraries was done using output from the 10X Cell Ranger pipeline, using Cell Ranger 3.0 count and aggr functions with the MiSeq ClickTag fastqs as input, and combining the counts from the two lanes using the denoted sample IDs. A ClickTag count matrix (ClickTag x cell) was generated by counting ClickTag barcodes that had high sequence similarity to the designed sequences using the Python fuzzywuzzy package to identify targets within Levenshtein distance 1. [Code] We additionally quantified gene expression with the kallisto-bustools workflow [START_REF] Melsted | Modular and efficient pre-processing of single-cell RNA-seq[END_REF], which reproduced concordant results, described below.

Author Contributions

Conceived

Competing Interests statement

The authors declare no competing interests.

Software used

Cell Ranger 3.0.