Whole-animal multiplexed single-cell RNA-seq reveals transcriptional shifts across Clytia medusa cell types

Chari

et al

Sci. Adv.

7

, eabh1683 (2021) 26 November 2021

SCIENCE ADVANCES

|

RESEARCH ARTICLE

1 of 17

GENETICS

Whole-animal multiplexed single-cell RNA-seq reveals

transcriptional shifts across

Clytia

medusa cell types

Tara Chari

†, Brandon

Weissbourd

1,2,3

†, Jase

Gehring

†, Anna

Ferraioli

†, Lucas

Leclère

Makenna Herl

, Fan

Gao

, Sandra

Chevalier

, Richard R.

Copley

*, Evelyn

Houliston

David J.

Anderson

1,2,3

*, Lior

Pachter

1,8

We present an organism-wide, transcriptomic cell atlas of the hydrozoan medusa

Clytia hemisphaerica

and describe

how its component cell types respond to perturbation. Using multiplexed single-cell RNA sequencing, 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 per-

turbations. 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 that is readily adaptable to other traditional or nontraditional

model organisms.

INTRODUCTION

Single-cell RNA sequencing (scRNA-seq) is enabling the survey of

complete transcriptomes of thousands to millions of cells (

), re-

sulting in the establishment of cell atlases across whole organisms

(

–

), exploration of the diversity of cell types throughout the ani-

mal kingdom (

–

), and investigation of transcriptomic changes

under perturbation (

). However, scRNA-seq studies involving

multiple samples can be costly and may be confounded by batch

effects resulting from multiple distinct library preparations (

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 (

–

). These advances have created an opportunity

to explore systems biology of whole organisms at single-cell resolu-

tion, 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.

First,

Clytia

is a powerful, emerging model system spanning multiple

fields, from evolutionary and developmental biology to regenera-

tion and neuroscience (

–

). While previous work has characterized

a number of cell types in the

Clytia

medusa (

), a whole-organism

atlas of transcriptomic cell types has been lacking. Such an atlas is a

critical resource for the

Clytia

community and an important addi-

tion to the study of cell types across animal phylogeny.

Second, emerging multiplexing techniques present new oppor-

tunities for system-level studies of cell types and their changing

states at unprecedented resolution in whole organisms. The

Clytia

medusa offers an appealing platform for pioneering these 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 (

–

). Furthermore,

the 1-cm-diameter adult medusae used in this study contain on the

order of 10

cells, making it possible to sample cells comprehensively

across a whole animal in a cost-effective manner using current

scRNA-seq technology (fig. S1 and tables S1 and S2). 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 imple-

mentation to the data processing and analysis.

RESULTS

We compared control versus starved animals, as this strong, natural

istic stimulus was likely to cause notable, 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 4 days, and the second fed daily (see Mate-

rials and Methods). We observed numerous phenotypic changes in

starved animals, including a marked size reduction reflecting two-

to threefold fewer cells (Fig. 1, fig. S2, and see Materials and Methods)

(

), and a notable reduction in gonad size. Correspondingly, the

number of eggs released per day decreased (fig. S3) (

For scRNA-seq, single-cell suspensions were prepared from each

whole medusa and individually labeled with unique ClickTag bar-

codes (

) using a seawater (SW) compatible workflow (see Materials

and Methods, tables S2 and S3, Supplementary Methods, and fig. S4).

All labeled suspensions were pooled and processed with the 10X

Genomics V2.0 workflow and Illumina sequencing, allowing con-

struction of a combined dataset across organisms and treatments,

Division of Biology and Biological Engineering, California Institute of Technology,

Pasadena, CA 91125, USA.

Tianqiao and Chrissy Chen Institute for Neuroscience,

Pasadena, CA 91125, USA.

Howard Hughes Medical Institute, California Institute of

Technology, Pasadena, CA 91125, USA.

Department of Genome Sciences, University

of Washington, Seattle, WA 98195, USA.

Sorbonne Université, CNRS, Laboratoire

de Biologie du Développement de Villefranche-sur-mer (LBDV), 06230, France.

University of New Hampshire School of Law, Concord, NH 03301, USA.

Caltech

Bioinformatics Resource Center, California Institute of Technology, Pasadena, CA

91125, USA.

Department of Computing and Mathematical Sciences, California

Institute of Technology, Pasadena, CA 91125, USA.

*Corresponding author. Email: copley@obs-vlfr.fr (R.R.C.); houliston@obs-vlfr.fr (E.H.);

wuwei@caltech.edu (D.J.A.); lpachter@caltech.edu (L.P.)

†These authors contributed equally to this work.

The Authors, some

rights reserved;

exclusive licensee

American Association

for the Advancement

of Science. No claim to

original U.S. Government

Works. Distributed

under a Creative

Commons Attribution

NonCommercial

License 4.0 (CC BY-NC).

Downloaded from https://www.science.org at California Institute of Technology on December 07, 2021

Chari

et al

Sci. Adv.

7

, eabh1683 (2021) 26 November 2021

SCIENCE ADVANCES

|

RESEARCH ARTICLE

2 of 17

without requiring batch correction (Fig. 1, fig. S4, A to D, and table

S1). A total of 13,673 single-cell profiles derived from 10 individuals

(5 control and 5 starved) passed quality control, with high concor-

dance in cell type abundance and gene expression among animals in

the same treatment condition (see Materials and Methods and fig.

S5). From this gene expression matrix, we (i) derived a

Clytia

me-

dusa cell atlas and (ii) generated a high-resolution resource of the

transcriptional impact of starvation across all observed cell types.

To validate the cell atlas and assess technical variability across

and within multiplexed experiments, we performed an additional,

independent round of sequencing from 12 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 ~1 hour later. This

paradigm allowed us to identify candidate genes with IEG-like proper

ties across many cell types, including neurons (figs. S6 to S8 and

table S4) (

). IEGs are valuable tools in neuroscience, to identify

neurons that are active following a specific stimulus or behavior (

This methodology is thus able to detect transcriptional responses

across diverse stimulus-response paradigms (table S5).

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 (see Materials and Methods; Fig. 2, A and B; figs. S9 to

S11; and table S5). Each of the cell types was present in each of the

individual animals sequenced (fig. S12). We then generated a low-

dimensional representation (

) of these cell types (Fig. 2A). We

could group the cell types into seven broad classes (Fig. 2A) that

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 particu-

larly well characterized in

Hydra

, where they generate neural cells,

gland cells, and stinging cells (nematocytes), as well as germ cells

(

). 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 Materials and Methods and Fig. 2, B to D)

were concordant between the two separate multiplexed experiments

(see “Starvation” and “Stimulation” sections in Materials and Methods)

and robust to different transcriptome annotations (figs. S6 and 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 homol-

ogous genes in other animals, while for the others we performed in situ

hybridization for selected marker genes (Fig. 2C, figs. S11 and S14,

and table S3). Previously known cell types apparent in our data in-

cluded i-cells (

) and nematocytes at successive stages of differen-

tiation (

–

), as well as oocytes (

), gonad epidermis, manubrium

epidermis, and bioluminescent cells in the tentacles that each express

specific endogenous green fluorescent proteins (GFPs) (

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. 2, C and D, and fig.

S14) (

). Within known cell types, clustering revealed an unap-

preciated degree of cell heterogeneity, yielding novel subtypes. For

example, eight cell types could be distinguished within the gastro-

dermis, 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

(

). 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, dis

tributed across the main digestive compartments of the gastrodermis—

the manubrium, gonad, and tentacle bulb—and the gastrovascular

canals that link them (figs. S11 and S14). Comparison of gene modules

Tags

cDNA

Perturbation

Dissociation

and

ClickTag

multiplexing

Pooling

and

library

generation

Sequencing

Multi-

organism

data

analysis

Control

Starved (4 days)

Animal

Cells

Animals

Genes

Fig. 1. Overview of whole-organism multiplexed experimentation.

Experimen-

tal design of the starvation experiment showing (1) images of control versus 4-day

starved animals (scale bars, 0.5 cm), (2) dissociation of individual medusa and

chemical tagging of cells with ClickTags to enable multiplexed scRNA-seq, (3) pooling

of cells and library generation from lysed cells to generate (4) sequencing libraries

for the multiplexed cDNA and ClickTag data and create (5) single-cell resolved

gene expression count matrices from all animals (see Materials and Methods).

Downloaded from https://www.science.org at California Institute of Technology on December 07, 2021

Chari

et al

Sci. Adv.

7

, eabh1683 (2021) 26 November 2021

SCIENCE ADVANCES

|

RESEARCH ARTICLE

3 of 17

Gastroderm

Epidermal/muscle

Gland cell

Stem cell/

germ cell

Nematocyte

Neural

Bioluminescent

cells

Top marker genes

AB

C

Broad cell type categories

D

( )

0. 0

1. 0

-scored expression

Radial

canal

Circular

canal

Velu

Subumbrella

Exumbrella

Tentacle

Manubrium

Gonad

Mesoglea

Tentacle

bulb

3 GastroDigestiv

e-A

7 GastroDigestiv

e-B

14 GastroDigestiv

15 GastroDigestiv

e-D

19 GastroDigestiv

e-E

24 GastroDigestiv

e-F

ntacle bulb distal gastroder

33 Endodermal plate

ntacle GFP cell

1 Exumbrella epidermis

4 Manubrium epidermi

8 Striated muscle of subumbrella

18 Ra

dial smooth muscle

ntacle epidermis

29 Gonad epidermis

30 Striated muscle of velu

6 Neural cell early stages

9 Neural cells A (incl. GL

, MIH cells)

26 Neural cells B (incl. RF

amide cells)

31 Neural cells C (incl.

amide cells

)

22 Gland cells-A

27 Gland cells-B

25 Gland cells

-C

32 Gland cells-D

34 Gland cells-E

0 i-

cells

ry early oo

tes

13 Small ooct

2 Medium oo

tes

12 Nemato

te precursors

11 Early nematoblasts

23 M

id nematoblasts

17 Late nematoblasts

21 Mature nematoc

ytes

ll types

within

each class

Whole medusa

Manubrium

Gonad

Tentacle bulb

FibColl-B

TPM-B

TPM-A

CathepsinL

Chitinase

ShKT-

ypA

ST MyHCb

FibCdom-

CathepsinL

ST MyHCb

ST MyHCa

FibColl-B

Znf845

amide

BP10-lik

Tentacle

Mout

ST MyHC

Nanos1

Mos3

Cathepsin

Nematocilin

Fig. 2. The

Clytia

Medusa Cell Atlas.

(

) Two-dimensional UMAP embedding of cells labeled by seven cell type classes. Class colors are retained in (

) to (

). (B) Heatmap

of top marker genes from the sequencing data with 36 Louvain clusters comprising the seven cell type classes. (C) In situ hybridization patterns for a selection of cluster

marker genes providing spatial location on the animal (comprehensive set in fig. S14). The label GD denotes general markers for GastroDigestive cell types. Scale bars, 100



(D) Schematics of

Clytia

medusa, manubrium, gonad, and tentacle bulb showing the main cell types. Abbreviations of cell class names: GD, GastroDigestive; BC, bioluminescent

cells; EM, epidermal/muscle; GC, gland cells; SC, stem cell/germ cell; NY, nematocytes; NE, neural cells.

Downloaded from https://www.science.org at California Institute of Technology on December 07, 2021

Chari

et al

Sci. Adv.

7

, eabh1683 (2021) 26 November 2021

SCIENCE ADVANCES

|

RESEARCH ARTICLE

4 of 17

discriminating these gastrodermal clusters (fig. S15) indicates that

GD-B may have a particular role in transforming growth factor–



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 (

). GD-C

cells, localized closest to the endodermal plate, are enriched for

transcripts associated with extracellular matrix and mesoglea (jelly)

production (Fig. 2, C and D, and figs. S11 and 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 mix-

tures 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 identi-

ties 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)

(

). 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 (

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 (

). Nematocilin, a

hydrozoan-specific component of the nematocil (ciliary trigger for

nematocyte discharge) (

), is also expressed in these clusters (fig.

S14, table S5, and 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 and

fig. S14), a notable 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 (

) and at

other sites (

). 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 (

). As in

Hydra

, our cell atlas revealed clear

connections between the neuronal and nematocyte populations and

the i-cell population (Fig. 2A and figs. S11 and S14) (

), likely

corresponding to differentiation trajectories (

). 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

(

) (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 transdifferentia-

tion (

). 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 pseudo-time values to the cells and ranked genes in each

trajectory (see Materials and Methods and Fig. 3A). This revealed

trajectories consistent with these cell types both deriving from

i-cells (Fig. 3, A and B) (

). Examination of the nematocyte trajectory

revealed the early expression of genes previously not associated

with this process, including

Znf845

and

Mos3

(Fig. 3, C and 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

(Fig. 3, C and D; corresponding expression domains in Fig. 3E and

fig. S14) (

–

). The trajectory analysis confirmed continuity

between the “nematoblast” clusters and the distinct and under-

appreciated nematocyte differentiation phase, characterized by ex-

pression of putative nematocil structural proteins and nematocilin

expression at the end of the trajectory (Fig. 3, C and D, and figs. S11

and 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 and 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

among transcriptomes derived from dissected

Clytia

bulb and ten-

tacle regions (

) and between developing and mature nematocyte

scRNA-seq clusters in

Hydra

(

Cnidarian nervous systems represent both valuable points of

phylogenetic comparison and tractable platforms for systems neuro

science (

). 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

that 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 (Fig. 4, A to C, and table S5) (

–

). In

mature neurons, neuropeptides are thought to be the dominant

neurotransmitters in cnidarians (

) but are challenging to

identify because of rapid sequence evolution (

). In conjunc-

tion with sequence-based analysis, we were able to identify 10 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. 2, B and D),

increasing the number of predicted

Clytia

neuropeptides to 21

(table S3). Our pseudo-time 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 reclustered the neural supergroup (“Neural;”

Fig. 2, A and B) to characterize neural subtypes. This distinguished

14 subpopulations of neurons and a progenitor population, ex-

pressing cell cycle and conserved neurodevelopmental genes including

the

bHLH

transcription factor Neurogenin (subcluster 0; Fig. 4D).

Notably, the neuronal subpopulations show combinatorial neuro-

peptide precursor expression, often with a distinct and identifying

Downloaded from https://www.science.org at California Institute of Technology on December 07, 2021