Molecular characterization of cell types in the squid Loligo vulgaris

  1. Jules Duruz
  2. Marta Sprecher
  3. Jenifer C Kaldun
  4. Al-Sayed Al-Soudy
  5. Heidi EL Lischer
  6. Geert van Geest
  7. Pamela Nicholson
  8. Rémy Bruggmann
  9. Simon G Sprecher  Is a corresponding author
  1. Department of Biology, Institute of Zoology, University of Fribourg, Switzerland
  2. Interfaculty Bioinformatics Unit and Swiss Institute of Bioinformatics, University of Bern, Switzerland
  3. Institute of Genetics, University of Bern, Switzerland

Abstract

Cephalopods are set apart from other mollusks by their advanced behavioral abilities and the complexity of their nervous systems. Because of the great evolutionary distance that separates vertebrates from cephalopods, it is evident that higher cognitive features have evolved separately in these clades despite the similarities that they share. Alongside their complex behavioral abilities, cephalopods have evolved specialized cells and tissues, such as the chromatophores for camouflage or suckers to grasp prey. Despite significant progress in genome and transcriptome sequencing, the molecular identities of cell types in cephalopods remain largely unknown. We here combine single-cell transcriptomics with in situ gene expression analysis to uncover cell type diversity in the European squid Loligo vulgaris. We describe cell types that are conserved with other phyla such as neurons, muscles, or connective tissues but also cephalopod-specific cells, such as chromatophores or sucker cells. Moreover, we investigate major components of the squid nervous system including progenitor and developing cells, differentiated cells of the brain and optic lobes, as well as sensory systems of the head. Our study provides a molecular assessment for conserved and novel cell types in cephalopods and a framework for mapping the nervous system of L. vulgaris.

Editor's evaluation

This article describes cell types in the head of the squid, Loligo vulgaris, through expression patterns of key genes identified in single-cell transcriptomics. This topic is generally of great comparative interest. The data presented here are convincing, and these valuable findings will contribute to a better understanding of the cephalopod nervous and sensory systems, providing a basis for future comparative and evolutionary research.

https://doi.org/10.7554/eLife.80670.sa0

Introduction

Cephalopods are part of a long tradition in evolutionary neuroscience and are often highlighted as an ideal example of convergent evolution of nervous systems (Young, 1965; Young, 1985; Abbott et al., 1995). Their large brains, relative to body size, and behavioral repertoires offer opportunities for comparative analysis of the evolution of learning and memory mechanisms. Cephalopods possess a sophisticated nervous system that resembles that of vertebrates in relative size and cognitive abilities (Young, 1991; Hochner et al., 2006). These complex behavioral abilities include various forms of learning, problem-solving, and body color and pattern changes used either as predatory or anti-predatory behaviors or for social communication (Fiorito et al., 1990; Langridge et al., 2007; Staudinger et al., 2013; Darmaillacq et al., 2014; Schnell et al., 2015; Richter et al., 2016; Lin and Chiao, 2017; Marini et al., 2017; Mather and Dickel, 2017; Hanlon and Messenger, 2018). The fossil record suggests that cephalopods were more prominent in marine ecosystems over 300 million years ago with a higher biodiversity than today. It is assumed that the rise of bony fishes and resulting competition for resources in that ecological niche played a major role in the extinction of many clades of cephalopods, including all the shelled cephalopods except few members of the subclass Nautiloidea. Today, the class Cephalopoda has over 800 current living species with extremely diversified lifestyles and body shapes (Kröger et al., 2011). The vast majority of current cephalopod species belong to the soft-bodied subclass Coleoidea (cuttlefishes, octopuses, and squids), and a small handful belong to the subclass Nautiloidea (nautiluses). Coleoid cephalopods, which have internalized or lost their shell over evolutionary time, were able to swim more easily and could compete with other predators in open waters. The continuous competition with other animals to access resources has been proposed to have been the main trigger for the selection of cephalopods that could develop specialized strategies for accessing resources unavailable to competing species. In coleoid cephalopods, this evolutionary process resulted in larger brains, camera-type eyes, as well as abilities to voluntarily change their skin pattern and texture for camouflage and social communication (Young, 1973; Hanlon and Messenger, 2018; Amodio et al., 2019).

The European squid Loligo vulgaris is well known for its economic role in the food and fishing industry but is much less studied for its fascinating biology. This species is known for its striking ability for camouflage typically in open-water environments by the use of its specialized iridescent cells to reflect light in different manners depending on its surroundings. The neuroanatomy and nervous system morphology of L. vulgaris have been described in detail in histological studies (Young, 1976; Wild et al., 2015). Detailed staging of the embryonic development of the related species Loligo pealii was established and can be used as a reference for L. vulgaris (Arnold, 1965). While the gross anatomy of squids is well known from histology and morphological studies, we are currently lacking resolution at the molecular level to identify genes that specify the many cell types that make up the unique cephalopod body-plan and nervous system. Since the publication of the first complete genome sequence of Octopus bimaculoides, there has been increasing interest in cephalopod biology and neuroscience. In addition, the availability of genomic sequences of other cephalopods, such as Euprymna scolopes (Sepiida), Architeuthis dux (Oegopsida), and Nautilus pompilius (Nautilida), opens the possibility for molecular comparisons between these distant clades of cephalopods with drastically different morphologies and lifestyles. These comparisons constitute a major step for understanding the molecular basis of cephalopod evolution (Albertin et al., 2015; Belcaid et al., 2019; da Fonseca et al., 2020; Zhang et al., 2021).

Single-cell transcriptomics experiments in different cephalopods have shed light on the cell type diversity of the O. bimaculoides visual system, the brain of Octopus vulgaris, and the nervous system of the squid Euprymna berryi (Gavriouchkina et al., 2022; Songco-Casey et al., 2022; Styfhals et al., 2022).

To gain further insights into the cellular diversity and cell type-specific gene expression in squids, we performed single-cell transcriptomics on whole L. vulgaris heads (stages 28–30 based on Arnold, 1965, pre-hatchling stage). Additionally, we used hybridization chain reaction in situ hybridization (HCR ISH) to visualize gene expression in situ and correlate single-cell transcriptomes with anatomical features. By combining these methods, we were able to identify and assess common cell types found in most bilaterians, such as different types of neurons, muscles, and connective tissues. Interestingly, even in these canonical tissues, we found certain clusters that highly express cephalopod-specific genes, as exemplified by some subtypes of connective tissues. Thus, while cells can be assigned a certain identity based on well-characterized marker genes, such as collagens, these cell types have likely further diverged in cephalopods. Neural tissues are prominent in the head and include the brain, optic ganglia, eyes, and other sensory organs. We identified and analyzed new markers of different cells in the nervous system, including widely expressed pan-neuronal genes, genes marking distinct stages of neural differentiation, markers for neural subsets, as well as sensory cell types. In addition, we also characterized cephalopod-specific cell types, such as the chromatophores or cells of the suckers. Our study provides new insights into the transcriptomic signatures of cell types in L. vulgaris. By analyzing the expression patterns of many cell type-specific genes, we assessed the conservation of features previously described in other animal clades but also previously unknown features and modalities of cell types that appear specific for cephalopods.

Results

Single-cell transcriptome of the L. vulgaris head

To generate a reference transcriptome of L. vulgaris, embryos were collected at late embryonic stages for RNA extraction and subsequent cDNA synthesis. The transcriptome was sequenced, assembled de novo, and annotated (for details, see 'Materials and methods').

Single-cell RNA sequencing was performed from the dissociated cells of the whole heads of L. vulgaris pre-hatchlings. Cells were captured and sequenced in two different experiments and the resulting barcoded sequences were mapped to the reference transcriptome. The mapping rate of the sequences was of 29.3% and 29.8% for samples 1 and 2, respectively (Figure 1—figure supplement 1A). This low mapping rate can be explained by L. vulgaris genes not being mapped to the transcriptome because the genes are missing from the transcriptome and by the presence of cells and ambient RNA from other species in the single-cell suspension. The assessment of the transcriptome assembly using benchmark universal single-copy orthologs (BUSCO) showed a score of 54.4% (49.8% single-copy, 4.6% duplicated, 1.9% fragmented), which indicates that a number of genes are likely missing from our transcriptome that in turn may impact downstream analysis. Cells that were identified and sequenced were clustered and analyzed. 20,000 cells were analyzed in our study. Cells were filtered to keep only the ones with a gene-per-cell count comprised between 200 and 2000 (median = 700 genes/cell). This filtering step resulted in a total of 19,974 cells. Cells were plotted to identify the presence of outliers for gene-per-cell count, UMI-per-cell count, or percentage of mitochondrial genes (Figure 1—figure supplement 1B). Clustering was performed with Seurat using 25 principal components with a resolution value of 2. This resulted in a total of 34 cell clusters numbered from 0 to 33. Clusters were assigned to presumed cell type identities based on the genes differentially expressed in each cluster (Cl). Clusters were manually annotated and fitted into cell type categories based on the expression of specific markers (Figure 1B, Figure 1—figure supplement 2). These defined cell type categories include muscle, connective tissue, stem cells, epidermis, papilin+ cells, neurons, foxD1+ cells, reflectin+ cells, ASc+ cells, cilia-related cells, photoreceptors, sucker, and lens.

Figure 1 with 2 supplements see all
Single-cell transcriptome of the L. vulgaris head.

(A) Ventral view of a recently hatched L. vulgaris with relevant external anatomical structures labeled. The dotted line represents the plane of amputation. The region anterior to this line is what was used in this study and is referred to as the head region. (B) Two-dimensional UMAP representing the 34 cell clusters and their cell type characterization based on marker gene expression. These defined cell type categories include muscle, connective tissue, stem cells, epidermis, papilin+ cells, neurons, foxD1+ cells, reflectin+ cells, ASc+ cells, cilia-related cells, photoreceptors, sucker, and lens. (C) Successive confocal micrographs of the ventral–dorsal view of L. vulgaris hatchling heads stained with DAPI (nuclei). The entire ventral–dorsal surface of the head of L. vulgaris is labeled with relevant anatomical structures: (a) labels a1–4 = arms 1–4, (bl) brachial lobe, (dbl) dorsal basal lobe, (de) dorsal epidermis, (fl) frontal lobe, (fu) funnel, (ibl) inferior buccal lobe, (L1 ,L2) epidermal lines, (Oo) olfactory organ, (ol) optic lobe, (pr) proximal segment of retinas, (s) suckers, (sbl) superior buccal lobe, (svl) subvertical lobe, (t) tentacles, and (ve) ventral epidermis. Scale bars = 500 μM.

The identification of muscle cells is based on the expression of various genes involved in the formation of contractile fibers such as myosin heavy chain, troponin, myosin regulatory light chain, paramyosin, tropomysin (Cl00). A second muscle cluster (Cl25) expresses markers including myosin catalytic light chain and paramyosin but is also distinguished by the expression of the transcription factor soxE. Smooth muscle-like cells (Cl25) are characterized by the expression of the homeobox transcription factor nkx2-5, which plays a crucial role in cardiac development in human embryos and angiopoietin-1 receptor, involved in angiogenesis in many species (Elliott et al., 2006). Because of the presence of these two markers, we hypothesize that these cells could be a part of the squid circulatory system.

Connective tissues are defined by the very broad expression of many different types of collagens, which are one of the main components of connective tissues (Cl01, Cl06, Cl23, Cl28). The presence of enriched laminin and fibroblast growth factor receptor further suggests that these cell clusters include precursors or differentiated cells of connective tissues. These markers that are typical of connective tissues are, however, co-expressed with several uncharacterized genes that have known orthologs in other cephalopod species. These previously unidentified markers could provide new avenues for the study of cephalopod connective tissues and provide a molecular complement to the existing morphological data (Thompson and Kier, 2001; Kier and Stella, 2007).

Cells of the epidermis are characterized based on the expression of an egf-like epidermal growth factor (Cl08, Cl09). Interestingly, Cl08 comprises chromatophores defined by the co-expression of tyrosinase and tryptophan 2,3-dioxygenase, two members of the melanin synthesis pathway that, when knocked out, have been shown to result in depigmentation in Doryteuthis pealeii (Crawford et al., 2020). Cl31 was characterized by a different epithelial marker epithelial splicing regulatory protein (esrp).

Three clusters (Cl03, Cl04, Cl12) highly express papilin, a gene which was suggested to play a role in the formation of extracellular matrix in Drosophila melanogaster (Fessler et al., 2004; Apte, 2020). The function of this gene in other clades is, however not known.

Cells of the suckers (Cl29) could be identified by the enriched expression of suckerin genes that have been shown to be expressed in the sucker ring teeth of squids (Kumar et al., 2016).

Neurons were identified based on widely used markers, including two genes belonging to the ELAV RNA-binding protein family elav-like 1 and elav-like 2 (Cl05, Cl07, Cl11, Cl13, Cl15). Known markers for the presynaptic machinery such as synaptobrevin, synaptogyrin, synaptophysin, and synapse-associated protein 1 are also broadly expressed in these clusters as well as amyloid beta-binding-like, tetraspanin-8, and rab3.

Cl27 is characterized by the expression of different cilia-related and flagella-related markers, together with the neuronal gene rab3 and synaptogyrin-1. The modality of these cells remains unknown, but the presence of cilia/flagella-related genes together with neuronal markers could suggest a mechanosensory function (see Figure 4).

Photoreceptors (Cl32) could be identified based on the expression of various members of the canonical phototransduction cascade, including trp-like, visual arrestin, retinal binding protein, phospholipase C, and rhodopsin.

Cl16, Cl19, and Cl30 are very clearly defined by the expression of the forkhead box transcription factor foxD1. In other phyla, this gene is involved in various developmental processes. However, because the presence of this gene is not sufficient to affirm a specific cell type, we named these cells foxD+ cells. A large proportion of these cells co-express the transcription factor ets4. In Cl30, foxD1 was co-expressed with a GATA zinc-finger protein (gata-2) and a few blood-related markers (ferritin, hemocyte, apolipophorin).

Cl20 and Cl22 are characterized by the expression of different types of reflectin genes, a family of proteins used in cephalopods to control light reflection by modulating the iridescence of the skin (Crookes et al., 2004). Interestingly, the transcriptional repressor scratch-2, which is highly expressed in several neuronal clusters, is also highly expressed in reflectin+ cells. We also observe expression of reflectins in neurons in agreement with this notion.

Different types of neural progenitors (Cl21, Cl26) are characterized by the high expression of a soxB1 and an achaete-scute (asc). Genes of the achaete-scute complex have been broadly studied for their implication in the specification of the neuronal lineage, particularly in the fruit fly (Campuzano and Modolell, 1992; Skeath and Doe, 1996; Zhao et al., 2007; Arefin et al., 2019; Soares et al., 2022). Both clusters also express histone-related genes and markers of proliferation, suggesting that they may be replicating undifferentiated cells. The presence of a soxB1 in these cells indicates that they may be proliferating neuronal progenitors, based on previous evidence of expression of soxB1 genes during neurogenesis in cephalopods (Miyagi et al., 2009; Focareta and Cole, 2016).

The lens is a specialized part of the eye for which we identified a cluster (Cl33) based on the expression of s-crystallin, an important structural component of the tissue of the lens (Tomarev and Piatigorsky, 1996).

To next move beyond the in silico assessment of cell cluster identity and describe where distinct cell types are localized, we analyze the expression of an array of marker genes specifically enriched in cell clusters of particular interest using HCR ISH. We further assess the newly identified cell types and describe their localization anatomically (Figure 1C).

Molecular mapping of the L. vulgaris nervous system

The central nervous system of L. vulgaris could be anatomically characterized using DAPI staining and were reinforced by HCR ISH experiments using broad neuronal markers. The anatomical characterization of brain structures is usually relying on the identification of neuropils, which are regions that have large numbers of axons and synapses but few cell bodies. Since the methods used here to describe gene expression are detecting RNA molecules that are mostly located in the cell bodies, usually surrounding nuclei, we briefly describe here the anatomy of the L. vulgaris nervous system in a nonexhaustive manner to facilitate the understanding of the results. The largest structure in the L. vulgaris nervous system are the optic lobes. The cortex of the optic lobes is composed of an outer and inner plexiform layers and outer/inner granular layers. The rest of the optic lobe is referred to as the medulla. The suboesophageal mass and supraoesophagal mass make up the central part of the brain and contain several different lobes. Discussed in this article from anterior to posterior are the inferior buccal lobe, superior buccal lobe, frontal lobe, anterior basal lobe, precommissural lobe, vertical lobe, medial basal lobe, and dorsal basal lobe. The smaller peduncle lobe and olfactory lobe are posterior to the optic lobes. Note that the nervous system of L. vulgaris has several posterior lobes (posterior to the basal lobes) that are not shown on the confocal micrographs of this study as they have been sometimes lost during sample preparation after amputation of the head. A rough schematic of the arrangement of the different lobes (not at scale) can be found in Figure 2—figure supplement 1.

To search for neuronal subtypes in the nervous system, we looked specifically at the cell clusters that express neuronal makers (Figure 2A). We then assessed the expression patterns of neuronal markers using HCR ISH to gain better insights into the molecular identities of neuronal populations. Probes were generated for the genes elav-like 1, elav-like 2, amyloid beta-binding-like, and tetraspanin-8, which were all found to be expressed broadly across all neuronal clusters (Figure 2B). HCR ISH experiments for these probes confirmed their expression throughout the nervous system of L. vulgaris, but different genes showed specificities in their expression domains and expression levels. Some regions of the nervous system show co-expression of elav-like 1 and elav-like 2 (Figure 2C). Both genes are highly expressed in the optic lobes and in the neurons of the supraesophageal and suboesophageal mass. However, differences in the expression patterns of elav-like 1 and elav-like 2 can be observed: elav-like 1 is strongly expressed in the medial part of the supraesophageal mass, including neurons of the dorsal basal lobe, vertical lobe, and subvertical lobe, as well as in the optic lobes in the medial portion of the medulla (Figure 2C, Figure 2—figure supplement 2). Elav-like 2, on the other hand, is expressed broadly in the lateral edges of the supraesophageal mass near the optic lobes and in the optic lobes. Interestingly, photoreceptors of the retina do not express elav-like 1 but only elav-like 2 (Figure 2C, Figure 2—figure supplement 2). The different expression domains are reflected by the scRNAseq data that shows, for instance, that elav-like 2 expression level is higher than that of elav-like 1 in Cl05, whereas the opposite can be observed in Cl07 (Figure 2B, Figure 2—figure supplement 2). Differences in the expression of these two ELAV homologs suggest that these genes could either serve different functions specific to the neuronal subpopulations in which they are expressed or that the different expression patterns reflect distinct developmental or differentiation states. The data shows that elav-like 2-expressing cells are present in the central and peripheral nervous system as well as in the lateral lips, which are shown to be connected to stem cells and neuronal differentiation later in the article (see Figure 3).

Figure 2 with 8 supplements see all
Molecular mapping of the L. vulgaris nervous system.

(A) UMAP with clustered identified as neuronal or neuron-related highlighted. (B) Heatmap showing the expression of neuronal genes across each cell of the dataset. Cells are grouped by clusters. (C) Confocal micrographs showing mRNA expression detected by hybridization chain reaction (HCR) for broad neuronal markers identified in the scRNAseq data. The right-hand pictures are feature plots (UMAP) showing the expression of genes across the subclustered neuronal cells. Cartoons on the top right corner of each micrograph indicate the position of the image on a frontal view (top) and the plane of acquisition on sagittal view (bottom). Scale bars = 250 μM. (D) Confocal micrographs showing mRNA expression detected by HCR for additional neuronal markers; some of these markers were not differentially expressed in specific clusters. Note that DAPI staining in some cases has not penetrated the tissue, therefore explaining its absence from the optic lobes. Cartoons on the top-right corner of each picture indicate the position of the image on a frontal view (top) and the plane of acquisition on sagittal view (bottom). Scale bars = 250μM. (abl) anterior basal lobe, (dbl) dorsal basal lobe, (fl) frontal lobe, (ibl) inferior buccal lobe, (ll) lateral lip, (mbl) medial basal lobe, (ol) optic lobe, (Oo) olfactory organ, (ogl) outer granular layer, (pcl) precommissural lobe, (pdl) peduncle lobe, and (sbl) superior buccal lobe.

Figure 3 with 5 supplements see all
Stem cells and neuronal differentiation.

(A) UMAP with clusters identified as either stem cell, Asc+, or FoxD+ highlighted. (B) Dot plot showing the average expression of marker genes and proportion of cells in each cluster that express these genes. (C) Confocal micrographs showing mRNA expression detected by hybridization chain reaction (HCR) for genes involved in stem cells and neuronal differentiation. Cartoons on the top-right corner of each picture indicate the position of the image on a frontal view (top) and the plane of acquisition on sagittal view (bottom). Red circles highlight the location of the lateral lips. Scale bars = 250μM. (co) cortex of the optic lobe, (eye) eye, (ll) lateral lip, (me) medulla of the optic lobe, and (ol) optic lobe. (D) Graphical representation summarizing the gene expression patterns of the genes assessed by HCR in situ hybridization (ISH) in the head of L. vulgaris with dorsal view (top) and transversal view (bottom). P/A, posterior/anterior; D/V, dorsal/ventral. (E) Plots showing the expression patterns of indicated genes across neurons and stem cells. The clusters and the trajectories are determined using monocle.

Expression of the gene tetraspanin-8 was assessed and observed to be present in all cells throughout the brain, providing tetraspanin-8 as a strongly expressed pan-neuronal marker. This high neural expression pattern suggests that tetraspanin-8 is for the function of differentiated neurons in L. vulgaris providing an interesting candidate for the study of neuronal function in cephalopods. The pan-neuronal expression of tetraspanin-8 allowed precise description of neuroanatomical regions of the L. vulgaris central nervous system (Figure 2—figure supplement 3A–C). Note that in many pictures shown in this study, DAPI appears absent in nervous tissue. This is the result of a lack of penetrability of DAPI inside the tissue but does not in any case signify the absence of nuclei in these structures.

Similarly, the expression of amyloid beta-binding-like is observed broadly throughout the entire brain, retina, and other components of the nervous system (Figure 2C). However, at higher magnification, it becomes evident that expression levels differ between cells. This is specifically apparent in the cortex layers of the optic lobe, where the cells of the inner granular layer express lower levels of amyloid beta-binding-like compared to tetraspanin-8 (Figure 2C, Figure 2—figure supplement 3D–G). Low expression of amyloid beta-binding-like can also be observed in the olfactory organ (Figure 2—figure supplement 3F). Additionally, the scRNAseq data shows that the expression of amyloid beta-binding-like protein is lower in the reflectin+ clusters (Cl20, Cl22, Figure 1—figure supplement 2).

The transcriptional repressor scratch-2 was found to be highly expressed in two neuronal clusters (Cl05, Cl15) as well as in a reflectin+ cluster (Cl22). HCR ISH expression analysis showed that this gene is mostly expressed in the medulla of the optic lobes and in the supraesophageal mass in two domains lateral–ventral to the vertical lobe and anterior to the basal lobe that extends anteriorly toward the frontal lobes, as well as a stripe of cells connecting the two domains posterior to the vertical lobe (Figure 2C, Figure 2—figure supplement 4). The analysis of the expression of scratch-2 in the scRNAseq dataset highlights that the cells highly expressing this gene are distinct from neurons expressing elav-like 1, tetraspanin-8, and amyloid beta-binding-like but show more transcriptional similarities with elav-like 2-expressing cells (Figure 2C).

Fmrf was shown to be expressed in some cells within the nervous system in agreement with previous anatomical studies (Wollesen et al., 2010; Wollesen et al., 2012). It is expressed in a region located between the dorsal basal lobe and the anterior basal lobe, surrounding the precommissural lobes, the peduncle lobe, the olfactory lobes, on the posterior part of the buccal mass and in the olfactory organ. Some cells expressing fmrf are also scattered throughout the optic lobes (Figure 2D, Figure 2—figure supplement 5).

The homeobox transcription factor LIM homeobox is expressed in the vertical lobe and subvertical lobe in a domain located between the dorsal basal lobes and the frontal lobes. It is also expressed around the peduncle lobe and in some cells at the surface of arms and tentacles (Figure 2D, Figure 2—figure supplement 6).

The expression of additional genes that were not identified as cell type specific in the scRNAseq but traditionally characterize neuronal populations based on their neurochemical modalities were assessed with HCR ISH. We identified a serotonin transporter gene that is expressed in the nervous system and more particularly in the medulla of the optic lobes. Low expression of this gene can be observed throughout the nervous system and indicates that the expression of the gene in the nervous system is not specific to an anatomical structure (Figure 2D). This transporter was, however, very highly expressed in the cornea of eyes (Figure 2—figure supplement 6D–F).

Two receptors for glutamate (glutr3) and GABA (gaba-r) were also tested with HCR ISH. The results show that a very large proportion of neurons express glutr3 but gaba-r expression is mostly restricted to the medial part of the brain and a few cells of the optic lobes (Figure 2D). More specifically, gaba-r is expressed in the region posterior to the dorsal basal lobes, between the anterior basal lobes and the frontal lobes, in the peduncle lobes and the superior buccal lobes (Figure 2—figure supplement 7). Gaba-r appears to also be highly expressed in both neuronal and non-neural parts of the eyes but is absent from the optic lobes. Expression of glutr3, on the other hand, is observed all over the animals in neuronal tissues but also in epithelial and muscular tissues. This extremely broad expression pattern for glutr3 indicates that this gene is likely to have other functions unrelated to neuronal communication. The genes involved in neurotransmission normally used to differentiate neuronal populations were expressed only in few cells in the single-cell dataset (Figure 2—figure supplement 8). This could be caused by issues with the quality of the cells or a lack of mapping to the reference transcriptome. Possible reasons for these issues are discussed later in the article.

Stem cells and neuronal differentiation

Single-cell sequencing makes it possible to identify cells not only based on their functional cell type identity but also based on their developmental stage. By identifying transcription factors that have been shown to be involved in cell type specification in previous studies and by following the expression of these genes across different cells, single-cell transcriptomics can provide insights into the different developmental lineages that specify cell types in L. vulgaris.

Several marker genes expressed in specific clusters suggest an involvement in proliferation and possess some conserved markers of neuronal differentiation. Two clusters (Cl02, Cl14, Figure 3A) do not express specific transcription factors that can be associated with developmental processes but show enrichment of genes involved in cellular proliferation and protein synthesis (Figure 1—figure supplement 2). In Cl10 and Cl17, the transcription factor e2f3 is expressed together with other indicators of cellular proliferation, cellular growth, and nucleic acid synthesis such as histone proteins (H3), ribosomal proteins (40s ribosomal protein S13, 40s ribosomal protein s19), pcna, and DNA replication licensing factors (mcm3, mcm4, mcm5). The role of e2f3 to promote cellular proliferation has been shown to be conserved between D. melanogaster, Caenorhabditis elegans. and mammalians (Dynlacht et al., 1994; Attwooll et al., 2004; DeGregori and Johnson, 2006). Because these clusters highly express proliferation markers but lack any distinctive features of differentiated cell types, we propose that they may be stem cells. To visualize the expression pattern of this gene, we performed HCR ISH experiments for e2f3. The results show that this gene is specifically expressed in the anatomical structure called the lateral lips, which have been described in O. vulgaris (Figure 3C and D, Figure 3—figure supplement 1, Deryckere et al., 2021). This suggests the presence and proliferation of stem cells or progenitors in the lateral lips of squids, consistently with the proposed neurogenic regions of O. vulgaris.

In Cl17, Cl21, and Cl26, the transcription factor achaete-scute (asc), which has a well-studied role in the regulation of nervous system development in D. melanogaster, is expressed together with similar proliferation markers as indicated previously (Cabrera et al., 1987; Romani et al., 1987; Skeath and Carroll, 1992). This is particularly interesting to assess the possible conservation of this gene’s function in nervous system development. The expression pattern revealed by HCR ISH showed a strong similarity with the expression of e2f3 in the lateral lips (Figure 3C and D, Figure 3—figure supplement 2). This suggests that neuronal progenitors and stem cells are co-located in the lateral lip where they proliferate. Frequent co-expression with soxB1 in these clusters further supports their neuronal fate because of the conserved involvement of soxB1 genes in neuronal development in very distantly related species such as C. elegans and mammals as well as in cephalopods (Arsic et al., 1998; Miyagi et al., 2009; Alqadah et al., 2015; Focareta and Cole, 2016).

We next analyzed the expression of soxB1, which revealed enriched expression in the lateral lips, similarly to e2f3 and asc. Interestingly, the pattern of expression of soxB1 is broader than asc and is also present in most lobes of the brain (Figure 3—figure supplement 3). Together, the overlapping expression patterns of e2f3, asc, and soxB1 in domains associated with cephalopod neurogenesis strongly suggest that stem cells and neuronal progenitors are present the lateral lips of L. vulgaris, suggesting that these cells are likely migrating to the nervous system at later differentiation stages. In agreement with a pro-neural function, these cells do not yet express markers of terminally differentiated neurons. It was previously shown that soxB1 is expressed throughout late embryogenesis in the developing eye of D. pealeii (Napoli et al., 2022). We indeed found comparable expression of soxB1 in the developing eye in L. vulgaris pre-hatchling and further detected expression in the developing olfactory organ, suggesting a broader role of soxB1 in sensory organ development or function (Figure 3—figure supplement 3).

Cl16, Cl19, and Cl30 are characterized by the expression of the forkhead box transcription factor foxD1. FoxD genes have been documented to be expressed in the nervous systems of Mus musculus, Xenopus laevi, C. elegans, and D. melanogaster and to play a role in the development of neuronal cell types although not exclusively (Herrera et al., 2004; Polevoy et al., 2017; Newman et al., 2018; Janssen et al., 2022). HCR ISH assessment of the expression of foxD1 showed that these cells are strongly expressed in the cortex of the optic lobes and in the medial basal and anterior basal lobes (Figure 3C). The expression of foxD1 is very strong in the cortex of the optic lobes and in the medial and anterior basal lobes but is also sparsely located in the medulla (Figure 3—figure supplement 2). The expression pattern of foxD1 indicated that these cell types are present in the nervous system and could therefore be types of neurons or neuronal progenitors. Cl30 also expressed a GATA zinc finger transcription factor (gata-2) together with some markers typically associated with blood cells (hemocyte protein, immunoglobulin, soma-ferritin, apolipophorin). The expression pattern of gata-2 is very similar to foxD1 and is also strongly expressed in the cortex of the optic lobes (Figure 3C, Figure 3—figure supplement 3). Gata-2 expression can also be found at the surface of the arms and in the suckers. The HCR ISH experiments with these different genes showed that the expression of e2f3 and asc is limited to the proliferating cells located in the lateral lips, whereas the expression of foxD1 and gata-2 is not found in the lateral lips but inside nervous tissue. The expression pattern of soxB1 revealed that it is present both in the lateral lips and the nervous tissue of the optic lobes (Figure 3D).

To assess the possible developmental lineages that could link these different clusters together, we performed trajectory analysis on a subset of cells. We selected all the cells from the stem cell clusters (Cl02, Cl10, Cl14, Cl17), the neuronal clusters including photoreceptors and presumed sensory cells (Cl05, Cl07, Cl11, Cl13, Cl27, Cl32), the foxD+ clusters (Cl16, Cl19, Cl30), and the asc+ clusters (Cl21, Cl26). These cells were used for subclustering and subsequent trajectory analysis using Monocle3 (Trapnell et al., 2014; Qiu et al., 2017b; Qiu et al., 2017a). As monocle operates with its own clustering algorithm, we fist verified that the markers tested previously were still clustered together and then proceed to calculate trajectories. We selected the e2f3 expressing cluster as a starting point for the trajectories and the results indicated that following the proliferative stage where cells express e2F3 and asc, cell trajectories go through a stage where foxD1 and gata-2 are strongly expressed (Figure 3E). After this stage, we see increased expression of the markers of differentiated neurons amyloid beta-binding-like (Figure 3E). The expression of soxB1 along the trajectory axis appears to be present in the early stages as well as in the differentiated neurons, consistently with the expression patterns previously observed by HCR ISH (Figure 3C and E). The expression level of soxB1 is, however, very low in the foxD1-expressing stage. Additional markers could be plotted along the trajectory axis: the proliferation markers pcna and histone H3 and the transcription factor elf-2 show similar expression as e2f3 and asc in the early stages of differentiation (Figure 3—figure supplement 4). The hypothetical intermediate stage characterized by the expression of foxD1 and gata-2 could also be correlated with the expression of the transcription factors ets4 and hes-4 as well as chemotrypsin (Figure 3—figure supplement 4). Elav-like 1, elav-like 2, and synaptogyrin-1 are highly expressed in the differentiated, and although they are present across all clusters, the average expression increases along the differentiation axis (Figure 3—figure supplement 4).

To support the results shown by the monocle 3 trajectory inference, we performed a similar analysis with slingshot (Street et al., 2018). The same cells as for the monocle analysis were selected and re-clustered. The results obtained with slingshot were very similar to those obtained with monocle and showed a trajectory that supports foxD1+ cells as a possible intermediate stage between the proliferative stage where e2f3 and asc are expressed and the differentiated neurons expressing amyloid beta-binding-like and higher levels of elav-like 1 and elav-like 2 (Figure 3—figure supplement 5). The calculation of trajectories using these two separate methods provides insights into the developmental relationship between these different identified cell types and provides new hypotheses for the study of neurogenesis of cephalopods.

Sensory cells

Vision is a crucial sense for various behaviors of the squids as it is needed for navigating efficiently in open waters in the wild. The eyes of squids are very large and the optic lobes, where visual information is processed and integrated, take up a significant part of the overall volume of the brain. The phototransduction cascade that allows light stimuli to be transformed into an electrical signal in the neurons is well conserved among cephalopods, and we were able to clearly identify photoreceptors in the scRNAseq data (Davies et al., 1996; Monk et al., 1996; Arendt, 2003). Many important components of the phototransduction cascade are highly expressed in Cl32 such as trp-like, retinochrome, arrestin, phospholipase C, and rhodopsin. To verify the expression of members of the phototransduction cascade in the photoreceptors, HCR ISH probes for trp-like and rhodopsin were selected and tested. Both genes showed striking specificity and are expressed in the photoreceptors of the retina (Figure 4A, Figure 4—figure supplement 1). The olfactory organs of L. vulgaris are located on the ventral surface of the head (see Figure 1C). The genes involved in the function of these sensory organs remain unknown and thus no cell type-specific markers were readily available to identify chemosensory cells of the olfactory organ. However, through screening of some cluster markers of interest, we could observe that the olfactory organ strongly expresses FMRF-related neuropeptides (fmrf). Interestingly, the pattern of expression of fmrf in the olfactory organ (Figure 4B, Figure 2—figure supplement 4) strongly resembles that of reflectin-8 (Figure 4B, Figure 4—figure supplement 2), a gene typically involved in the function of iridophores, which are a specialized subclass of chromatophores. FMRFamide has been previously shown to be involved in the control of chromatophores in squids (Loi and Tublitz, 2000; Sweedler et al., 2000; Satpute-Krishnan et al., 2006; Mobley et al., 2008). Additionally, we observed expression of elav-like 1 in neurons of the epidermal lines, stripes of cells located on the dorsal epidermis of cephalopods that have been proposed to have a function similar to the lateral lines of fish (Figure 4B, Figure 1—figure supplement 1; Budelmann and Bleckmann, 1988). Epidermal lines were suggested to be mechanosensory cells, used to sense water currents in the cuttlefish Sepia officinalis and in the squid Lolliguncula brevis (Komak et al., 2005). The expression of elav-like 1 in this structure in L. vulgaris is consistent with their proposed mechanosensory function other cephalopods.

Figure 4 with 2 supplements see all
Sensory cells.

(A) Confocal micrographs showing mRNA expression detected by hybridization chain reaction (HCR) for genes expressed in photoreceptors. Scale bars = 250μM. (B) Confocal micrographs showing mRNA expression detected by HCR for genes expressed in chemosensory organs and mechanosensory cells. Cartoons on the top-right corner of each picture indicate the position of the image on a frontal view (top) and the plane of acquisition on sagittal view (bottom). Scale bars = 250 μM. (C) Graphical representation summarizing the expression of identified genes in specific sensory organs and regions. Dorsal view (top) and transversal view (bottom). P/A, posterior/anterior; D/V, dorsal/ventral. (D) UMAP of the subclustered neuronal subset with the plotted expression level of the different identified sensory markers. (a3) arm pair 3, (de) dorsal epidermis, (L1), epidermal lines 1, (ol) optic lobe, (Oo) olfactory organ, (ret) retina.

We observed expression of cilia-associated protein (Cl27) at the surface of the arms (Figure 4B). The scRNAseq data revealed that the cells of this cluster frequently co-express the synaptic markers rab3 and synaptogyrin-1. Because of the expression of cilia-associated protein at the surface of the arms, co-expressed with neuronal markers and additional cilia/flagella markers such as kinectin-like and dynein light chain, we propose that these cells could be specific mechanoreceptors of the arms of L. vulgaris. The expression of marker genes used to identify sensory cell types is summarized in Figure 4C. The sensory markers that were identified were plotted onto the neuronal subclustered previously used in Figure 2. The results highlight the co-expression of neuronal markers and sensory markers, therefore supporting their sensory identity (Figure 4D).

Epidermis, connective tissue, and cephalopod-specific cell types

The epidermis of cephalopods is covered with pigmented cells organized in precise patterns. Additionally, the cephalopod-specific chromatophores, which are colored or iridescent, have been shown to be modulated by the nervous system (Dubas et al., 1986; Wardill et al., 2012; Gonzalez-Bellido et al., 2014). Cl08 and Cl09 could be identified as epidermal cells from the scRNAseq data because of the consistent expression of epidermal growth factor EGF-like. In Cl08 and Cl09, egf-like was co-expressed with an aminopeptidase (xaa-pro aminopeptidase) and feeding circuit activating peptide-like, a neuropeptide that was shown to be involved in the feeding circuit in the nervous system of Aplysia californica (Sweedler et al., 2002). Cl08 expressed the markers tyrosinase and tryptophan 2–3 dioxygenase, two enzymes involved in melanin synthesis pathway (Takeuchi et al., 2005; Crawford et al., 2020). HCR ISH experiments confirmed that egf-like is expressed in several cells across the epidermis of the animal (Figure 5B, Figure 5—figure supplement 2). Cl31 was characterized as epithelial due to the very specific expression of epithelial splicing regulatory protein (esrp). HCR ISH for this gene revealed that it is expressed specifically on the surface of the arms, on the epithelium covering the surface of the eyes and in cells of the epidermis, similarly to egf-like. Esrp is also expressed at a lower level across the surface of the epidermis (Figure 5B, Figure 5—figure supplement 1). Cl03, Cl04, and Cl12 showed specific enrichment for papilin. Studies on the role of papilin in D. melanogaster revealed that it is a secreted extracellular protein that acts as a metalloproteinase and is important for embryonic development (Kramerova et al., 2000; Kramerova et al., 2003). In L. vulgaris as well as other cephalopods, the presence of this protein is documented but its function is unknown. These clustered also showed elevated expression of gluthatione S-transferase, indicating a possible role of these cells in metabolizing reactive oxygen species. The location of papilin-expressing cells was assessed with HCR ISH, and we observed that this gene is expressed in a very particular set of cells of the skin that form characteristic rings around larger epidermal cells (Figure 5B). Although papilin is expressed everywhere on the epidermis of the head of L. vulgaris, it is absent from the surface of the olfactory organ (Figure 5—figure supplement 2). This suggests that the layering of the olfactory organ is composed of different cell types than the rest of the epidermis. Cl28 is characterized by the expression of a zinc metalloproteinase nas 14-like of unknown function. Due to the co-expression of collagen-related genes, myosins. and actins, these cells were assigned to connective tissue. HCR ISH experiments showed expression of nas 14-like throughout the body in large cells with extended processes (Figure 5C). The branched morphology of these cells is similar to fibroblasts. Nas-14 is expressed predominantly in the ventral epidermis but is also observed in the dorsal epidermis, the arms, and tentacles and is heavily expressed on the surface of the funnel (Figure 5—figure supplement 3). We analyzed the expression of reflectin-8, one of the several reflectins expressed in Cl20 and Cl22. However, in addition to having been shown to be also expressed in the olfactory organ (see Figure 4), HCR ISH for this reflectin showed broad expression at the surface of the epidermis but also inside the optic lobes and around the anterior basal lobes, superior buccal lobes, and peduncle lobes (Figure 5C, Figure 4—figure supplement 2). The expression pattern of reflectin-8 in the nervous system suggests that the involvement of this protein may go beyond its ability to modulate light reflection and could play a role inside the nervous system itself for the modulation of chromatophores and iridophores as is suggested in previous literature (Demski, 1992). This is also supported by the co-expression of reflectins and the neuronal markers scratch-2, elav-like 1, elav-like 2, and rhodopsin-1 in Cl20 and the expression of LIM homeobox in Cl22. A mesencephalic astrocyte-derived neurotrophic factor (madnf) is expressed in Cl23, Cl27, and Cl31, which were characterized as connective tissue, mechanosensory cells, and epidermis, respectively. HCR ISH experiments for this gene surprisingly show an expression pattern very similar to reflectin-8. The results showed that this gene is expressed throughout the nervous system, including in the optic lobes in the medial portion of medulla (Figure 5—figure supplement 4). While madnf-expressing cells of Cl27 co-express neuronal markers, the cells of Cl23 and Cl31 do not. Given the observed expression pattern of this gene, it is possible that these cells are a non-neuronal cell type of the nervous system and could be a type of glial cells. Suckers contain a specific cell type that makes up the hard structure of the sucker ring teeth, which are used for holding on to their prey. The sucker ring teeth consist of the protein suckerin, for which we found gene expression in Cl29, which is associated with epidermal tissues. The cells of the suckers expressed many uncharacterized genes, indicating that these cell types are highly specialized in cephalopods. Moreover, we noted co-expression of myeloperoxidase-like, a gene involved in antimicrobial activity. It has been shown that O. vulgaris suckers express peptides that display structural similarities with antimicrobial α-peptides (AMPs) and exhibited a significant role in antimicrobial activities (Maselli et al., 2020). It is likely that octopuses and other cephalopods exploit secreted AMPs as part of an innate defense mechanism, similarly to other aquatic animals (Thøgersen et al., 1992; Derby, 2014; Falanga et al., 2016; Gogineni and Hamann, 2018). HCR ISH experiments showed that the expression of suckerin-6 is specific to a cup-shaped group of cells in the inner sucker, likely becoming cells producing the sucker ring teeth (Figure 5C).

Figure 5 with 4 supplements see all
Epidermis, connective tissue, and cephalopod-specific cell types.

(A) Dot plot showing the average expression of marker genes of interest and proportion of cells in each cluster that express these genes. Clusters identified as epidermal, connective tissue, or reflectin+ are labeled in red. (B) Confocal micrographs showing mRNA expression detected by hybridization chain reaction (HCR) for genes involved in epidermis development or function. Cartoons on the top-right corner of each picture indicate the position of the image on a frontal view (top) and the plane of acquisition on sagittal view (bottom). Scale bars = 250μM. (C) Confocal micrographs showing mRNA expression detected by HCR for genes expressed in connective tissue or associated with iridophores and suckers. Cartoons on the top-right corner of each picture indicate the position of the image on a frontal view (top) and the plane of acquisition on sagittal view (bottom). Scale bars = 250 μM. (a3) arm 3, (dbl) dorsal basal lobe, (fl) frontal lobe, (L1) epidermal lines 1, (mbl) medial basal lobe, (ol) optic lobe, (Oo) olfactory organ, (s) suckers, (sbl) superior buccal lobe, (t) tentacles, (ve) ventral epidermis.

Discussion

We here provide molecular insights into the cell types present in the head of L. vulgaris by combining two state-of-the-art techniques HCR ISH and scRNAseq. The combination of these two techniques allowed us to provide insights into the different cell types that make up the L. vulgaris head, into the process of nervous system development and the identity of neuronal subtypes. The data produced for this article enabled the identification of 34 cell clusters that we used as a correlate for cell type diversity. Although the data can provide a broad overview of cell diversity, we note that few cell types that would be expected in the head of cephalopods such as components of the vascular and digestive systems could not be confidently identified in this study. The apparent lack of resolution in the clusters and their low number is likely to be in large part caused by the lack of completeness of the reference transcriptome, which resulted in a lower mapping rate of the single-cell reads. It has been shown in other cephalopods that having a high-quality genome largely contributes to improving cluster resolution (Gavriouchkina et al., 2022; Songco-Casey et al., 2022; Styfhals et al., 2022). Because of the lack of a reference genome for L. vulgaris, a non-redundant transcriptome was used as a mapping reference in this study, which has shown to be sufficient for identifying cell types and describing the expression some of their characteristic marker genes but is however not sufficient to affirm full coverage of cell type diversity of this organism. The use of a transcriptome instead of a genome as a reference has been previously shown to enable molecular description of cell types but can lack in resolution (Duruz et al., 2021). Additionally, it is likely that the transcriptome assembly could be improved by obtaining more RNA from L. vulgaris by including different embryonic stages and adult stages. We are convinced that re-running this analysis with a better reference would reinforce our analysis and provide better cell type resolution. In our analysis of the nervous system, we were able to recognize broad neuronal markers that are expressed throughout the nervous system while showing different levels of expression in predicted neuronal cell types. These newly identified neuronal markers will provide a new perspective to study the characteristics of the nervous systems of cephalopods. Particularly, the genes tetraspanin-8 and amyloid beta-binding-like show strikingly elevated expression pan-neuronally. Besides the identified sensory cells, the homogeneity of the neuronal clusters was rather surprising. Different neuronal subpopulations were expected to be identified by the expression of genes neurotransmitter synthesis pathway to distinguish neurons by their chemical modalities. Such fine distinctions between neuronal populations could be found in the nervous systems of E. berryi and O. vulgaris (Gavriouchkina et al., 2022; Songco-Casey et al., 2022; Styfhals et al., 2022). It is likely that this lack of resolution is due to missing genes in our reference transcriptome, which would classically contribute to distinguishing neuronal populations. The analysis of the expression of two highly expressed members of the ELAV family elav-like 1 and elav-like 2 revealed that although they can be found in most lobes of the nervous system, their expression levels differ in different clusters. Previous studies of ELAV proteins in S. officinalis revealed that the gene Sof-Elav1 was expressed in all the lobes of the nervous system but also during early embryogenesis while the other gene Sof-Elav2 was not specific to neurons (Buresi et al., 2013). In O. bimaculoides, a single neuronal ELAV protein was identified and was observed to be also expressed during neurogenesis and could be used as a neuronal marker to study nervous system development during embryogenesis, similarly to Sepia (Shigeno et al., 2015). It is interesting to note that contrarily to what was described in cuttlefish and octopus we show two cell type-specific elav genes, which are both expressed in neurons. This challenges the hypothesis of a duplication of the elav gene during lophotrochozoan evolution proposed after the observation of two elav genes in Capitella sp., one of which is expressed in neurons while the other is expressed mainly in the mesoderm (Meyer and Seaver, 2009). The identification of two neuronal elav genes in L. vulgaris raises the question of whether they are orthologous to the genes studies in Sepia and Octopus but their function is different. Alternatively, these two neuronal elav genes could the result of another duplication of the neuronal elav. Such questions would need to be further investigated by more thoroughly studying the sequence homology of these genes and reconstruct their phylogeny. We identified key transcription factors involved in stem cells and development. The expression of these genes likely reflects different states of neuronal differentiation and suggests a migration of neuronal progenitors towards the nervous system. We could show that proliferating cells that express of the transcription factor e2f3 are located in the lateral lips near the nervous system. We interestingly observed the expression of asc in the same region, together with several proliferation markers. Extensive study of neurogenesis in D. melanogaster has shown that the proneural complex asc can determine the neuronal fate of quiescent multipotent ectodermal cells (Campuzano and Modolell, 1992; Skeath and Doe, 1996). In vertebrates, on the other hand, asc is expressed in cells that have already undergone neuronal specification but are self-renewing (Bertrand et al., 2002; Soares et al., 2022). Studies of neurogenesis in the annelid Capitella teleta have shown that, similarly to vetrebrates, proneural genes such as asc are expressed in self-renewing neuronal progenitors (Sur et al., 2020). Our data suggests that asc is also expressed in proliferating cells, similiarly to vertrebrates and C. teleta but unlike Drosophila. This is consistent with documented co-expression of asc and the mitotic marker phosphor-histone H3 in O. vulgaris (Deryckere et al., 2021). The observation of soxB1 expression together with asc strengthened their characterization as neuronal progenitors due to the documented role of this gene in neuronal differentiation (Miyagi et al., 2009). However, strong soxB1 expression in differentiated neurons has been shown in the gastropod Patella vulgate and in the cephalopods S. officinalis and O. vulgaris (Le Gouar et al., 2004Focareta and Cole, 2016; Deryckere et al., 2021). In L. vulgaris, we observe similar results where soxB1 is expressed in the lateral lips together with asc but is also expressed in differentiated neurons. This was consistently shown by the HCR ISH experiments as well as the single-cell analysis. Therefore, the hypothesis proposed for the action of soxB1 in the early stages as well as the later stages of neuronal differentiation in O. vulgaris is consistent with our data (Deryckere et al., 2021). The unexpected expression of foxD1 in the nervous system of L. vulgaris, together with existing evidence of the role of this gene in neuronal development, motivated the grouping of the foxD1+ clusters with other clusters characterized as neuronal progenitors for further analysis (Herrera et al., 2004; Polevoy et al., 2017; Newman et al., 2018; Janssen et al., 2022). Our trajectory inference analysis using monocle3 and slingshot showed that foxD1 is expressed in cells that could represent an intermediate neuronal differentiation stage following the expression of asc but before the expression of amyloid beta-binding-like and the second peak of soxB1 expression. Such a role has been proposed in O. vulgaris in the cells expressing neuroD that are expressed in a gradient between the lateral lip and the central nervous system, which may indicate a pattern of migration of differentiating neurons (Deryckere et al., 2021). We could, however, not correlate these results with our hypothesis of an intermediate foxD1-expressing stage because neuroD expression could not be detected in specific clusters in our single-cell dataset. Although the trajectories we calculated point towards a role of FoxD1 in neuronal differentiation, the relatively poor cluster resolution and the lack of cells along the trajectory axes shown in Figure 3 suggest that this hypothesis would require further investigation that should include a more complete set of genes that have been involved with cephalopod neurogenesis. In our analysis of sensory cells, we provide valuable information regarding the different sensory systems of L. vulgaris and could identify photoreceptors and additionally describe newly identified markers of the olfactory organ and mechanosensory cells. The sensory cells of cephalopods have been extensively studied and are known to be very complex and can have multiple sensory modalities (van Giesen et al., 2020; Buresch et al., 2022). In O. bimaculoides and S. officinalis, several chemosensory receptors could be identified and characterized (van Giesen et al., 2020; Andouche et al., 2021). In our analysis of the L. vulgaris head, we were however not able to identify putative chemosensory genes, likely due to the absence of characteristic receptors in our transcriptome. The chemosensory function of cells analyzed here was only inferred because of their expression in the olfactory organ. A full analysis of chemosensory genes in L. vulgaris would be required to be able to truly describe the diversity of their sensory cell types. Finally, we further characterize the highly specialized epidermis of L. vulgaris and cephalopod-specific cell types such as suckers or iridophores. The expression of a few genes facilitated the characterization of some of these cell types as exemplified by the expression of egf-like in epidermal cells, suckerin-6 in suckers, and reflectin-8 in iridophores. However, many cluster-specific genes could not be assigned to a specific function. The abundance of papilin-expressing cells on the surface of the epidermis is striking and had not been previously documented. The gene madnf is known to play a role in the maturation and maintenance of dopaminergic neurons in D. melanogaster and C. elegans (Palgi et al., 2009; Richman et al., 2018). We found discrepancies between the observation of madnf in non-neuronal clusters and in a small cluster of mechanosensory cells and the expression patterns shown by HCR ISH, where expression is found in different areas of the brain, optic lobes, and olfactory organs. Similarly, the expression of reflectin-8 is not limited to the chromatophores but is also in the nervous system. The expression of these non-neuronal makers inside the nervous system results highlights the importance of considering non-neuronal cell types, even when aiming to answer neurobiological questions. Many cell types remain unidentified, and further characterization of these cells and the function of the genes that they differentially express could advance our understanding of cephalopod biology. The combination of high-resolution imaging rendered possible by HCR in conjunction with the data provided by single-cell transcriptomics enabled us to extend our understanding of the cell types of cephalopods. The experimental approach is efficient to identify and explore conserved sets of genes conferring cell type identity in these animals while comparing cell type markers across different clades. On the other hand, the accurate anatomical description gained with HCR ISH allowed us to find previously undescribed cell type-specific genes by identifying the anatomical features in which they are expressed, enabling the exploration not only of conserved cell types but also novel ones. By applying single-cell RNA sequencing in a member of this fascinating clade of mollusks and by correlating this data with high-resolution mRNA imaging techniques, we can further understand the main molecular features that drive cell type identity in L. vulgaris. The analysis of the expression of selected genes of interest can provide insight into their function, which can be a starting point for future investigations of cephalopod cell types. Further studies of the nervous system of squids and the different cell types that compose them could help fill important gaps in our understanding of the evolution of the nervous system and could provide insights into how brains complexified independently in very distant clades of animals. By comparing the molecular features of neuronal cell types of cephalopods with other species, we could better understand the evolutionary history of the unique brains of cephalopods and identify what confers them such unique capabilities. This is a significant step in the study of convergent evolution of brains and provides a valuable resource for future research on cephalopods.

Materials and methods

L. vulgaris embryo collection and maintenance

Request a detailed protocol

L. vulgaris eggs were collected during a sampling project of the European Marine Biological Resource Center (EMBRC) France at the Bay of Morlaix. The eggs were then incubated in a closed artificial seawater system (Aquarium Red Sea system reefer deluxe 170 Black). The salinity was between 34 and 36 particles per thousand (ppt), and pH was between 7.80 and 8.40. A 12 hr light/dark cycle was set with a crepuscular light and a dawn light before the day and the dark cycle; the used lamps were the Maxspect éclairage LED de la series RSX 100W. The eggs were suspended per bundle in a plastic structure made with PVC. Each bundle entailed between 15 and 30 arms. Each arm contained between 100 and 150 embryos. Great care was allocated to guarantee egg oxygenation; according to Le Gouar et al., 2004; Focareta and Cole, 2016; Deryckere et al., 2021, the oxygen consumption of one egg mass is twice the consumption of an adult animal of mean size. Therefore, two additional bubblers with air diffusers were added (Tetra APS 100 and Tetratec AS45). The stage of the animals was assessed based on the description of embryonic stages in Loligo paelii (Arnold, 1965). As the temperature was lowered to slow down development, the stage 28–30 animals (pre-hatchlings) were determined based on the presence of chromatophores but lack of visible ink sac and the size of the yolk sack approximately equivalent to the size of the head (Arnold, 1965).

RNA extraction

Request a detailed protocol

Ten embryos were collected into QIAzol Lysis Reagent (QIAGEN, Cat# 79306), homogenized, and centrifuged. RNA was isolated with chloroform and precipitated with isopropanol at –20°C. After washing with ethanol, the RNA was resolved in DEPC-treated water.

Transcriptome assembly

Request a detailed protocol

The quantity and quality of the extracted RNA were assessed using a Thermo Fisher Scientific fluorometer (Qubit 4.0) with the Qubit RNA BR Assay Kit (Thermo Fisher Scientific, Q10211) and an Advanced Analytical Fragment Analyzer System using a Fragment Analyzer RNA Kit (Agilent, DNF-471), respectively. The RNA was also tested by spectrophotometry using a Denovix DS-11 FX Spectrophotometer/Fluorometer to assess the purity of the RNA. Once all quality control tests confirmed high-quality RNA, the 'Procedure & Checklist – Iso-Seq Express Template Preparation for Sequel and Sequel II Systems' was followed (PN 101-763-800 version 02 [October 2019]). Two libraries were created from the same total RNA input (300 ng): one targeting 2 kb transcripts and one targeting >3 kb transcripts. Both libraries were SMRT sequenced using a Sequel binding plate 3.0, sequel sequencing plate 3.0 with a 20 hr movie time on a PacBio Sequel system on their own SMRT cell 1M v3 LR. The 2 kb library was loaded at 3 pM and generated 33.90 Gb and 409,187 ≥ Q20 Circular consensus reads (CCS). While the >3 kb library was loaded at 3 pM and produced 27.04 Gb and 248,304 ≥ Q20 CCS. All steps post RNA extraction were performed at the Next Generation Sequencing Platform, University of Bern. High-quality isoforms were generated from PacBio reads using the Iso-Seq analysis pipeline of single-molecule real-time (SMRT) Link (version 9), which generates circular consensus sequencing reads and then clusters and polishes the isoforms. Cd-hit-est version 4.6.8 (Fu et al., 2012) was used to remove redundancy due to different transcript isoforms. It was run using a sequence identity threshold of 95%. Afterward, the Coding GENome reconstruction Tool (cogent v 6.0.0, Huo et al., 2020) was used to further reduce redundancy and reconstruct the coding genome. In order to annotate the coding genome, the redundancy removed sequences were blasted against the nr database of NCBI using blastx version 2.9.0 (NCBI Resource Coordinators 2016). We aligned our reference non-redundant transcriptome to the benchmark universal single-copy orthologs database for Mollusca (BUSCOs) and obtained a score of 54.4% (49.8% single-copy, 4.6% duplicated, 1.9% fragmented). This indicated that the reference transcriptome has low redundancy (<5% duplicates) that is necessary for single-cell transcriptomic analysis but, on the other hand, is missing several genes.

Cell dissociation for single-cell RNA sequencing

Request a detailed protocol

Embryos were first dissected out of their eggs in stages 28–30 (pre-hatchling stage) and their heads were isolated and put in nuclease-free phosphate buffer (PBS) on ice. The heads were then centrifuged for 5 min at 4°C, and the supernatant was carefully removed and replaced with a papain enzyme solution at a concentration 1 mg/ml in nuclease-free PBS. The tissue was then incubated for 30 min at room temperature (20–25°C) under continuous agitation. During the incubation period, the solution was pipetted up and down every 10 min to facilitate dissociation. The cell suspension was then centrifuged at 2000 rpm for 5 min at 4°C. The supernatant was removed and replaced by 1 ml of PBS containing 0.04% bovine serum albumin (BSA) to stop the enzymatic reaction. This last step was repeated one more time to ensure removal of the enzyme. The suspension was then filtered through a 40 μm Flowmi Cell strainer (Bel-Art H13680-0040). The filtered suspension was centrifuged at 2000 rpm for 5 min at 4°C. The supernatant was discarded, and cells were suspended in 100 μl of PBS 0.04% BSA with added RNAse inhibitor (0.1 μL/ml) and further dissociation was ensured by gently pipetting the entire volume approximately 200 times. Cell concentration and viability were assessed using a DeNovix CellDrop Automated Cell Counter with an Acridine Orange (AO)/propidium iodide (PI) assay.

Single-cell capture and sequencing

Request a detailed protocol

Single-cell RNAseq libraries were prepared using a Chromium Single Cell 3′ Library & Gel Bead Kit v2 or v3 (10X Genomics), according to the manufacturer’s protocol (user guide). Two chips were loaded with the accurate volumes calculated based on the 'Cell Suspension Volume Calculator Table.' The initial single-cell suspension being estimated at >600 cells/μl, we targeted to recover a maximum 10,000 cells. Once GEMs were obtained, reverse transcription and cDNA amplification steps were performed. Sequencing was done on Illumina NovaSeq 6000 S2 flow cell generating paired-end reads. Different sequencing cycles were performed for the different reads, R1 and R2. R1 contained 10× barcodes and UMIs, in addition to an Illumina i7 index, and R2 contained the transcript-specific sequences. All steps were performed at the Next Generation Sequencing platform at the University of Bern.

Single-cell RNAseq analysis

Request a detailed protocol

10X Genomics Cell Ranger version 3.0.2 (Huo et al., 2020) was used to map the single-cell RNAseq reads from 10X Genomics Chromium to the above assembled coding genome and count the number of reads overlapping with each gene to produce a count matrix for each sample. Afterward, the R package scater v. 1.14.0 (McCarthy et al., 2017) was used to automatically filter low-quality cells. Additionally, the gene NP-062835.1 was removed due to its abnormally high expression. The cells were filtered to keep only those with a gene per cell count comprised between 200 and 2000. The data was normalized based on a deconvolution method integrated into the R package scran v. 1.14.5 (Lun et al., 2016). Integration analysis, dimensionality reduction, and clustering were done using the R package Seurat v. 3.1.3 (Stuart et al., 2019). The data was normalized using the normalization method 'LogNormalize' with a scale factor of 10,000. The function 'FindVariableFeatures' was run using the selection method 'vst'. The data was then scaled, and we run principal component analysis function 'RunPCA'. Clustering was performed using the functions 'FindNeighbors' and 'FindClusters' using 25 dimensions and a resolution of 2. Upregulated genes for each cluster were identified using the function 'FindAllMarkers'. The 100 most differentially expressed markers are provided in Supplementary file 1. For the analysis of the nervous system, neurons were subclustered by selecting the cells of the clusters obtained previously that contained neuronal markers (Cl5, Cl7, Cl11, Cl13, Cl15, Cl20, Cl22, Cl27, Cl32). The cells were then reclustered with the same parameters as described above but using six dimensions and a resolution of 1.2. The 100 most differentially expressed markers in these neuronal subclusters are provided in Supplementary file 2. Subclustering of the cells used for the developmental trajectories was selected by subsetting the clusters identified as either neuronal, asc+, foxD+, or stem cells (Cl5, Cl7, Cl10, Cl11, Cl13, Cl15, Cl16, Cl17, Cl18, Cl19, Cl21, Cl26, Cl27, Cl32). In Seurat, the cells were reclustered using the same parameters as above but using six dimensions and a resolution of 1. The 100 most differentially expressed markers obtained with Seurat in these subclusters are provided in Supplementary file 3.

Trajectory inference

Request a detailed protocol

Trajectory inference was used on the subset of cells that were reclustered from neurons, asc+ cells, foxD+ cells, or stem cells as described above. Monocle3 (version 1.2.9) was used to recluster the data and calculate trajectories. We used the R package 'Seurat_to_Monocle' to convert the Seurat object into a Monocle object (SeuratToMonocle3, version 0.1.0; https://github.com/Jackcava/Seurat_to_Monocle, Jackcava, 2022). The data was then preprocessed using the 'preprocess cds' function using PCA with 25 dimensions. Dimension reduction was performed using the 'reduce_dimension' function. Clustering was done using the function 'cluster_cells' with a resolution of 1e-3. Trajectories were calculated using the 'learn_graph' function and selecting Cl5 (see Figure 3—figure supplement 4) as the starting point. Trajectories were then plotted onto the UMAP using either cluster labels or pseudotime and the function "plot_cells" was used to show the expression of specific genes along the trajectory axis. The 100 most differentially expressed markers obtained using Monocle3 are provided in Supplementary file 4. Slingshot (Street et al., 2018, version 2.2.1) was used as an alternative to calculate trajectories. The 'slingshot' function was run with the embedded clustered Seurat object. Cl9 of the subclustered data (Figure 3—figure supplement 5) was selected as the starting cluster. The resulting trajectories were plotted onto the original clusters obtained with Seurat.

Hybridization chain reaction

Request a detailed protocol

Probes for HCR ISH were ordered and designed by Molecular Instruments (Choi et al., 2016). Whole pre-hatchlings L. vulgaris animals (stages 28–30, Arnold, 1965) were dissected out of their eggs and were then anesthetized in a mix of 1.7% MgCl2, 1% ethanol, and artificial seawater. After 20 min in the anesthetic solution, the animals were fixed in 4% formaldehyde for 48 hr at room temperature. Samples were then washed three times with phosphate buffer (PBS) and were dehydrated and stored in 100% ethanol at –20°C. Prior to HCR ISH experiments, samples were rehydrated in PBS in three 5 min steps. Our HCR ISH protocol used was adapted from the provided 'generic sample in solution' protocol (Choi et al., 2016) with the following modifications: a light proteinase K treatment was performed before pre-hybridization for 5 min at a concentration of 1 ul/ml in PBS at room temperature and subsequently washed three times with PBS. All volumes were reduced to 200 ul per well. Incubation times for hybridization and amplification were increased to 16–20 hr to increase signal (inappropriate for RNA quantification). DAPI was added together with one of the 5× SSCT washes. Double HCR ISH experiments were performed by combining B1-Alexa-488 hairpins with B2-Alexa-647 hairpins.

Confocal imaging

Request a detailed protocol

Images were taken on Leica Stellaris 8 falcon or a Leica TCS Sp5 laser scanning confocal microscopes. All mosaic images were taken by automatically stitching 3 × 3 or 3 × 4 confocal image stacks with a 10% overlap. Laser power and gain were adjusted depending on the experiment to ensure optimal visualization. Images were treated with ImageJ where the contrast was adjusted. All images are optimized for qualitative assessment and are therefore not appropriately acquired or treated for quantitative analysis or comparisons.

Data availability

Sequencing data have been deposited in GEO under accession code GSE200108.

The following data sets were generated

References

  1. Book
    1. Apte SS
    (2020) ADAMTS proteins: concepts, challenges, and prospects
    In: Apte SS, editors. ADAMTS Proteases: Methods and Protocols. New York, NY: Springer. pp. 1–12.
    https://doi.org/10.1007/978-1-4939-9698-8_1
    1. Arendt D
    (2003) Evolution of eyes and photoreceptor cell types
    The International Journal of Developmental Biology 47:563–571.
    https://doi.org/10.1387/ijdb.14756332
    1. Buresch KC
    2. Sklar K
    3. Chen JY
    4. Madden SR
    5. Mongil AS
    6. Wise GV
    7. Boal JG
    8. Hanlon RT
    (2022) Contact chemoreception in multi-modal sensing of prey by octopus
    Journal of Comparative Physiology. A, Neuroethology, Sensory, Neural, and Behavioral Physiology 208:435–442.
    https://doi.org/10.1007/s00359-022-01549-y
  2. Book
    1. Marini G
    2. De Sio F
    3. Ponte G
    4. Fiorito G
    (2017) 1.24 - behavioral analysis of learning and memory in cephalopods
    In: Byrne JH, editors. Learning and Memory: A Comprehensive Reference (Second Edition). Oxford: Academic Press. pp. 441–462.
    https://doi.org/10.1016/B978-0-12-809324-5.21024-9
    1. Young JZ
    (1976) The nervous system of Loligo. II. suboesophageal centres
    Philosophical Transactions of the Royal Society of London. Series B, Biological Sciences 274:101–167.
    https://doi.org/10.1098/rstb.1976.0041

Decision letter

  1. Sonia Sen
    Reviewing Editor; Tata Institute for Genetics and Society, India
  2. Didier YR Stainier
    Senior Editor; Max Planck Institute for Heart and Lung Research, Germany
  3. Sonia Sen
    Reviewer; Tata Institute for Genetics and Society, India
  4. Samuel Reiter
    Reviewer; Okinawa Institute of Science and Technology, Japan
  5. Astrid Deryckere
    Reviewer; Columbia University, United States

Our editorial process produces two outputs: (i) public reviews designed to be posted alongside the preprint for the benefit of readers; (ii) feedback on the manuscript for the authors, including requests for revisions, shown below. We also include an acceptance summary that explains what the editors found interesting or important about the work.

Decision letter after peer review:

Thank you for submitting your article "Molecular characterization of cell types in the squid Loligo vulgaris" for consideration by eLife. Your article has been reviewed by 3 peer reviewers, including Sonia Sen as the Reviewing Editor and Reviewer #1, and the evaluation has been overseen by a Reviewing Editor and Didier Stainier as the Senior Editor. The following individuals involved in the review of your submission have agreed to reveal their identity: Samuel Reiter (Reviewer #2); Astrid Deryckere (Reviewer #3).

The reviewers have discussed their reviews with one another, and the Reviewing Editor has drafted this to help you prepare a revised submission.

Essential revisions:

While appreciating the relevance of this there were a set of common concerns raised by the reviewers that need to be addressed. These are:

1. The reference transcriptome: A high-quality reference transcriptome will be critical for any downstream analysis. Could the authors please provide data for how they validated their reference transcriptome? For example, the authors could report the BUSCO score for assessing the completeness of their assembly, a cumulative density plot of AED scores for annotation quality, and any other metric they choose. Could the authors also report what percent of their single cell reads mapped to the reference?

2. The scRNAseq data: Could the authors please provide data that would allow readers to assess the quality of their single-cell data? At a minimum, could they please provide violin plots depicting the range of UMIs, genes, and mitochondrial genes per cell?

3. Cluster-specific markers: Could the authors please provide dot plots for the top differentially expressed genes for all clusters? In addition to this, the authors have described a set of genes that they did not identify as cell-specific markers (for example, neurotransmission-related ones). Could they please show dot plots for these as well?

4. Trajectory Analysis: Could the authors please explain how they have chosen the cells used in this analysis, the clustering resolution they have used (please provide a dimplot), and which starting point they have selected? In addition to this, could they please analyze the differentially expressed genes along the trajectory to make their claim stronger? Since this data challenges the current state of the field, it would be useful to demonstrate the robustness of their results by verifying it by:

4a. Using different trajectory analysis methods (Monocle or FateID).

4b. Analysing more marker genes (in addition to foxD1) along the differentiation axis.

5. The writing: Could the authors please pay attention to some aspects of the text to help the reader place their work in the broader context of the field? Specifically:

a. There are missing references (see detailed reviews below), including three other recently posted cephalopod scRNAseq preprints. Could the authors please incorporate these references?

b. On occasion, the chronology of the text and figures don't match. There are also some figures that are currently not mentioned in the text. Could the authors please fix this?

c. Many of the methods are not sufficiently well explained. Could the authors please pay attention to this, particularly for bioinformatics? Could the authors motivate this better?

d. The Results section contains many points that would belong better in the Discussion section.

e. Conversely, the discussion as it currently stands, is a summary of the results. Could the authors please revisit their discussion to place their work and findings in the context of the growing field of cephalopod neurobiology?

Review #1 (Recommendations for the authors):

1. I was concerned about the quality of the scRNAseq data for a few reasons:

- The number of genes per cells (200-2000) seems low to me.

- The clusters don't really resolve.

- Many genes that one would have expected to pick up, for example, neurotransmitter-related genes for the neuronal cluster, were not picked up.

So, could the authors please provide the violin plots depicting the range of UMIs, genes, and mitochondrial genes per cell?

2. It was not clear to me how the authors picked the marker genes. Could they provide the top 20 DEGs for each of the clusters?

3. The signal in each of the images could be enhanced. To orient the reader, it would also be useful to have an inset of each of the expression patterns in the whole head. Some of these are really nicely represented in the supplementary information.

4. The authors should refer to the three recent preprints that have used scRNAseq in other cephalopods and, in the discussion, place their work in this context.

5. In general, could the authors elaborate their method sections? Many fairly important aspects are too briefly described – for example, the reference transcriptome and its validation and the selection of marker genes, among others.

Reviewer #2 (Recommendations for the authors):

The core of this study is the analysis of scRNAseq data. Here the authors assembled a transcriptome de novo. I wonder about the quality of this assembly, as it will affect all downstream analysis. The other 3 studies (Styfhals et al. bioRxiv 2022.01.24.477459; doi: https://doi.org/10.1101/2022.01.24.477459; Songco-Casey et al. bioRxiv 2022.06.11.495763; doi: https://doi.org/10.1101/2022.06.11.495763; Gavriouchkina et al., bioRxiv 2022.05.26.490366; doi: https://doi.org/10.1101/2022.05.26.490366) all refer to the necessity of a high-quality reference genome for a high-quality transcriptome. If this is not the case here, it would be useful to quantify this.

I found the analysis of stem cells and progenitors quite interesting. The authors profile specific markers for cell types, map these cell types onto interesting anatomical regions, and analyze developmental trajectories defined transcriptomically.

In other parts of the paper, I was surprised that many cell types could not be defined based on certain marker genes/small combinations of genes. Figure 2b, for example, shows that clusters have a great deal of overlapping expression among 'marker' genes. This would deviate quite a bit from the other scRNAseq studies I am familiar with. I could see a few possible explanations: (1) This is a product of low-quality transcriptome assembly/single-cell data (2) This is a product of the immature developmental stage of the animal, and (3) This reflects a cephalopod-specific transcriptomic signature.

I doubt possibility 3, as the other 3 cephalopod scRNAseq papers describe cell types well differentiated (eg. By neurotransmitter expression, Figure 2 in Songco-Casey et al. 2022, Figure 3 in Gavriouchkina et al. 2022). The fact that specific neurotransmitters are often not well localized to specific cell types here calls for an explanation. Similarly, the lack of differentiation of neural types even when analyzed separately from other cell types (Figure 2c). Many phototransduction genes are expressed in Cluster 32, and these genes are expressed, specifically in the retina. They could potentially be used as a sanity test for the scRNAseq: Are these genes only found in cluster 32?

The lack of specific marker genes makes some of the FISH experiments in this manuscript less informative than others. Many plots show anatomical localization of a gene broadly expressed in many cell types, meaning we cannot conclude how transcriptomically identified clusters map spatially. In other cases, FISH is used to mark individual cell types. In some cases, it is not clear which situation we are looking at. It would be useful to have a quantification of how specifically the marker genes used in FISH experiments mark individual clusters.

The authors observe a similar anatomical expression of fmrf and reflectin-8 in the olfactory organ. Because these have been linked to chromatophore and iridiophore control, they argue that this suggests that chromatophores and iridiophores are modulated through the olfactory organs. I find this plausible but the logic questionable. FMRF has many functions unrelated to chromatophore modulation, so just observing its presence seems insufficient to make statements about the neural control of chromatophores. Reflectin-8 is shown later in the manuscript to be broadly expressed in the epidermis and brain (Figure 5c).

I was not clear why the expression of elav-like 1 in the epidermal lines 'informs on their potential mechanosensory function'. With elav-like 1 expressed very broadly across cell types and anatomical locations.

The expression of the cilia-associated protein, rab3 and synaptogyrin-1, and other cilia/flagella markers in cluster 27 suggested that this corresponds to a mechanosensory population. To test whether this is neural, the authors take the subset of cells in the dataset expressing amyloid β- binding-like (a neural marker) and see whether a subset of these also expresses the collection of sensory markers. This seems round about to me. Why not look directly at whether the cells in cluster 27 express amyloid β- binding-like? If the probability of co-expression is high, then concluding that these cells are neurons seems warranted. If the probability is lower, then the conclusion does not seem warranted. The author's existing analysis does not address the possibility that non-neurons (those not expressing amyloid β- binding-like) express the collection of sensory markers.

Line 63: 'Coleoid cephalopods, which have internalized or lost their shell over evolutionary time 64 were able to swim more easily and could exploit other niches.'

Not clear what 'other' refers to. The niches of ammonoids? Not clear whether it was the loss of shell that led to niche exploitation or the other way.

Line 93 "… as a first step to understanding cephalopod evolution". I'm not clear why a molecular comparison would be a first step to understanding cephalopod evolution. What then is the non-molecular evolutionary story given earlier in the introduction?

Figure 2 Are the Umap plots in 2c re-estimated on the subset of data in A? If so, I don't see a description of this in the methods section, and I am surprised at the lack of cell type differentiation among neurons.

Figure 4D is described as 'feature plots'. Is this UMAP? How are the plots on the right-hand side generated?

The developmental trajectory analysis using slingshot could be better explained in text/methods/figure legend and not easily interpretable from Figure 3e if one is not familiar with the technique

'Additionally, the cephalopod-specific chromatophores, which are colored or iridescent, have been suggested to be modulated by the nervous system'. Maybe I am misunderstanding the sentence, but as I read it this is an understatement of what we know. Over 50 years of work have profiled the neural control of chromatophores in detail.

Some of the supplementary figures are quite dark. Why not increase the brightness of the images? Others like Figure 4 supplement 1 are beautiful.

Reviewer #3 (Recommendations for the authors):

First of all, I would like to congratulate the authors on the extensive expression analysis study performed in this manuscript. However, I believe that this dataset offers so many more opportunities to better understand cell types in cephalopods, and in general, the data could be analysed and interpreted in more detail.

Since the field of scRNAseq is moving incredibly fast, the authors should be careful with claims of being first (e.g. Abstract line 34). Preprints of scRNAseq data in the nervous system of other cephalopods have been recently published (Styfhals 2022, Gavriouchkina 2022, Songco-Casey 2022).

The manuscript is missing key citations to support the results and to provide a comparative view of gene expression in other cephalopods. For example, the expression patterns of Elav, SoxB1, Asc have been reported during the development of octopus, squid, and cuttlefish, and could enhance the interpretation of the results observed in Loligo.

In general, the reader needs to be proficient in the brain areas of cephalopods in order to follow the results. Therefore, the manuscript would benefit from a short introduction of the main brain lobes in cephalopods to be accessible to a wider audience, either in the introduction or in the results in Figure 1C.

Please make sure to refer to the figures and figure panels chronologically. Multiple figure panels are not referred to throughout the text (e.g. the first figure described is Figure 1B and there is no notion to Figure 1A, similar to Figure 2, …).

The authors distinguish 33 cell clusters when using 25 pca's and a resolution of 2. The UMAP, however, shows most cells in one single, central blob. This is striking knowing how different connective tissue is from neurons. Did the authors try other methods for clustering and dimensionality reduction, or change the parameters in the current algorithm to try to resolve the data better? It would also be helpful to specify the filtering parameters for filtering out low-quality cells in the methods section.

In the first section of the results, the authors list a handful of enriched genes used to differentiate the major cell type categories. However, little explanation is given to why these genes were chosen and a comprehensive overview supporting the claims is missing. It would help the reader if they would add a plot (e.g. a DotPlot) to support these data. In addition, a UMAP plot with cells colored according to cell cycle could clearly show where the progenitor cells are located on the UMAP.

To avoid going back and forth between cell type categories in the first results paragraph, it would be helpful to describe Cl 27 (lines 202 onwards) and Cl 32 (lines 207 onwards) with the other nervous system cells, and Cl 29 (lines 213 onwards) with the epidermis section.

In order to interpret HCR expression patterns, the authors should add axis labels to the figures. In addition, I understand that the figure panels are small to add annotation of the brain regions, but since we are looking at different levels through the brain, either a schematic with the figure or annotation would really help.

The authors should explain the analysis performed to obtain the UMAP(?) plots in Figure 2C. Which cells were selected? Were they re-clustered? What do the clusters look like?

Did the authors perform image manipulation e.g. Figure 2C and 2D for the DAPI channel? Why is DAPI absent from the neurons in the tetraspanin-8 HCR, or the serotonin transporter staining?

In Figure 3B, it would be helpful to add some genes involved in proliferation, as exemplified in the text.

The authors claim that neural stem cells undergo neuronal differentiation and then migrate to the nervous system. However, evidence is missing as the authors do not show where progenitor cells stop dividing and start to differentiate.

The use of the marker foxD1 for differentiation neurons is not well supported. Additional marker analysis is required before this can be claimed.

The authors need to explain which clusters were used for trajectory analysis using Slingshot and how the data were reclustered. Were non-neuronal stem cells included? This might skew the analysis. Importantly, the method was not described in the methods section. It is critical for the readers to know which endpoints were selected in the slingshot algorithm. Off note, the slingshot is very dependent on clustering resolution and will try to connect all clusters. Have the authors tried a different pseudotime analysis method such as Monocle or FateID? Additionally, can the authors elaborate on which transcription factors are steering this differentiation? In the current state, the analysis is preliminary and contradicts present data on the origin and trajectory of neurons in the cephalopod brain (see comment below; Koenig et al. 2016, Deryckere et al. 2021).

With the current interpretation of the data, the authors claim that foxd1 cells are an intermediate developmental state of neurons. This would be very striking since the cortex of the optic lobe has been described as being the inner retina. Additionally, it is contradicting the expression of amyloid β-binding like in the optic lobe cortex, which the authors say is only expressed in differentiated neurons. Is amyloid β-binding-like expressed in clusters 15, 19, and 20? How about other markers for differentiated neurons? What would all those intermediate cells do there? Do the authors think they migrate into the medulla or even the central brain? This is contradicting existing literature on the origin of neurons in the cephalopod brain (Koenig 2016, Deryckere 2021) and requires in-depth discussion.

The authors prepared beautiful summary diagrams in Figures 3 and 4. However, they should be referred to in the main text.

After performing extensive HCR expression studies, can the authors go back to their scRNAseq data and better annotate the different clusters?

In order to better understand neuronal diversity and neural cell types (as the authors claim), the dataset has to be analyzed in more depth. 33 clusters/cell types seem low given the number of cells that were sequenced and the different tissues that were sampled. Can the authors identify the different brain regions when subclustering the neurons? Can they identify and annotate different neural cell types?

The authors present several striking findings. However, at the moment, the Discussion section explains the results obtained in the paper, but does not place the acquired data within the current state-of-the-art. Parts of the results should be moved to the discussion and the discussion should be elaborated on.

Specific comments:

Line 76: "by the use its specialized iridescent cells" misses an of

Line 97: it would be helpful for the readers to add the "pre-hatchling stage" to stage 28, as was done in the Results section for the audience to have a better idea of what stage was studied.

Throughout the text (starting line 219) please write "HCR ISH" instead of just "HCR". The chain reaction method is not specific to in situ hybridization and can also be used in immunohistochemistry and others.

Line 101: "in all animals" please be more specific; not all animals have neurons.

Line 110: "marking distinct stages of neural development", neural differentiation would be more appropriate.

Line 126: Add a reference to Figure 1A to show what the head region is.

Line 133 says 34 cell clusters were identified, but Figure 1B shows 33.

Line 145: "smooth muscle-like cells (Cl26)" should be Cl25.

Line 156-159: add a reference for this statement, maybe a specific example.

Line 167-171: "a gene whose function is homologous to …" and " exact function … is not known". This is contradictory. Do you mean that the gene's sequence is homologous?

Line 199-201: with the comments below on SoxB1 and Ascl1, this sentence needs to be revisited.

Line 239-240 and also further in the text: "in the medial portion of the optic lobes" do the authors mean the medulla (in contrast to the cortex (inner and outer granular layer))? Medial portion is not clearly defined and the cortex also stretches quite medially. (also line 278, 540).

Lines 263 and 271: please rearrange the supplementary figure so you can refer to panels A-D first. In general, please refer to figures chronologically.

Line 298 with Figure 2D: please switch panels of LIM and serotonin, so they follow the text.

Line 326: this part covers more than stem cells and progenitors (line 375 and onwards).

Line 327: please add an introductory sentence to this new paragraph; rephrase the first sentence (grammatically incorrect).

Line 338-340: this is an overstatement. What identity are the authors referring to? Genes involved in proliferation etc do not give identity to stem cells and Cl02 and Cl14 are also actively cycling. Do the authors suggest only Cl10 and Cl17 are stem cells? Can the authors show more evidence?

Lines 350-359: I would suggest moving this part to the discussion and incorporating the following comments:

Lines 347-352: The expression of Asc is not as well conserved as the authors suggest. Indeed, it has some proneuronal role, but the level of differentiation in the cell it is expressed is very different across the animal tree. In Drosophila as in other non-insect arthropods, Asc is expressed in quiescent ectodermal cells to drive them into the neural lineage (formation of a neural progenitor cell) (e.g. Cabrera et al., 1987; Skeath and Carroll, 1992). In contrast, in vertebrates and also in spiralia such as capitella teleta, Ascl1 is expressed in neural progenitor cells to drive them into differentiation (e.g. Bertrand 2002; Sur 2020). Ascl1 expression has been studied in octopus and seems to have a similar function as in capitella, not Drosophila (Deryckere 2021). This scRNAseq dataset has the power to test this hypothesis and figure out whether, in cephalopods, as is expressed in proliferating and not quiescent cells.

Lines 354-355: this conclusion does not reflect the data shown. The authors are showing proliferating stem cells, not the location of differentiation. In octopus, differentiation seems to happen in the transition zones (NeuroD cells, Deryckere 2021). It would be interesting to locate those cells in your dataset.

Lines 356-359: Again here, one needs to be careful to extrapolate the functionality of a gene in one species to other species. In invertebrates, SoxB1 is also expressed in differentiated neurons (e.g. Le Gouar et al., 2004; Semmler et al., 2010). Additionally, expression of SoxB1 has been described in Sepia (Focareta and Cole, 2016) and octopus (Deryckere 2021), indicating expression in neurons (with some interesting discrepancies). Furthermore, be careful with the use of the word "neuroblast" because it refers to different types of cells in different animals.

Line 361: "Interestingly" it seems a general theme for cephalopods, and maybe invertebrates in general (see comment above).

Line 367: Data are missing to show where the neural progenitor cells differentiate.

Line 370-374: in line with this, SoxB1 expression is also found in sensory epithelia of sepia and octopus (e.g. suckers and statocyst; Baratte, S. and Bonnaud 2009; Deryckere 2021).

Line 376-378: add a reference and elaborate. Regulators of development is a very general process. In addition "for this reason, these cells were presumed to be differentiating cells": this is an overstatement and needs additional gene expression patterns to support this presumption.

Line 382-385: optic lobe cortex = inner and outer granular layers, please adjust and rephrase. From Figure 3C upper panel, it seems that Foxg1 is also expressed in the medulla. Or is this all cortex? Annotation would be helpful here.

Line 386: this is an overstatement. Most neurons in the brain express elav1 and the function of foxD1 is very broad. Additional markers (like for example NeuroD) are required to support this claim.

Lines 394-400: Figure 3D is not mentioned in the text and seems substantial. The authors might also want to revise it with the comments above and in public review.

Lines 408-411: please add references to support these statements.

Lines 419-421: Why did the authors choose to map the expression of Pax9? Koenig et al. studied the eye in Doryteuthis and mapped the expression of several genes important for eye development. A reference to this work would be appropriate here. Pax9 is not an eye gene in vertebrates and invertebrates (in contrast to Pax2/6). In which clusters is it expressed (again here, a dotplot of the genes under study would be helpful)? The level of the image in Figure 4A for Pax9 is also different than the others, why is this? In which cells of the eye is Pax9 expressed?

Line 448: Figure 5B should be Figure 4B I think.

Line 457: can the authors verify that they selected cells with expression level >3? From the legend in Figure 4D, the expression level of 3 seems present in a lot of cells (it looks more like the selected level was 5). In addition, the authors should explain how they obtained the graphs on the right. It looks like they underwent an extra transformation? Or are they just distorted?

Line 459, and together with the comment above: It would be helpful to be able to check the cluster numbers. For example, are cilia-associated protein cells from Cl27? A dotplot or heatmap might be more convincing to support these statements.

In addition, please add a reference to Figure 4D here.

Line 468-470: egf-like does not seem to be expressed in Cl03 and 31 from Figure 5A as suggested in the text.

Lines 471-493: please add the discussed genes to the dotplot in Figure 5A or supplement.

Line 541: It seems that cells in Cl27 do express neuronal markers. Please revise.

Line 545 and onwards: please not only discuss the results obtained but place the discussion in the state-of-the-art literature. Several statements have already been proposed in other papers which deserve to be referenced (e.g. lines 559, 568). Multiple genes mapped with HCR in this manuscript have been described in other cephalopods and assigned a preliminary function (e.g. Buresi 2013, 2014; Shigeno 2015; Koenig 2016; Focareta 2016; Deryckere 2021 and others). Please acknowledge their presence in appropriate sections, and, if the authors which, they might add a paragraph describing consistency or differences in expression.

Line 562: in line with previous comments in the results, please revise. This is based on the functionality of asc in Drosophila, but not other animals. Additionally, in the opinion of the authors, what are the asc+ e2f3- cells?

Line 713: please add the parameters used to filter low-quality cells.

Line 722: please add additional information on how cells were subclustered (and re-clustered?). Please add a description of the slingshot analysis.

Figures:

In general, please rearrange so each panel comes chronologically in the Results section. Make sure to mention each panel in the text.

Add axes to the confocal images and consider annotation.

Figure 1:

A) Is this the ventral view instead of the dorsal view? The funnel is visible?

Please arrange the abbreviations alphabetically. Ibl and sbl are missing.

Figure 2:

Line 790: UMAP with clusters.

C) Add legend to the 1→6 expression? Bar. It is not clear how these featureplots were calculated. Which clusters were selected? Were they re-calculated? Please show the dimplot with cluster annotation. What is on the x- and y-axes? UMAP? t-SNE?

D) Please rephrase 'that were not represented in the scRNAseq data': do the authors mean not differentially expressed? Were the genes absent from the reference transcriptome?

Figure 3:

A) Why are Cl2 and Cl14 encircled?

C) Explain what the red encircled areas are.

Figure 4:

D) Double check that 3 is the threshold. The curly bracket is confusing (it reads like the right panels feed into the left). What are the x- and y-axes? Similar comment to Figure 2C.

https://doi.org/10.7554/eLife.80670.sa1

Author response

Essential revisions:

While appreciating the relevance of this there were a set of common concerns raised by the reviewers that need to be addressed. These are:

1. The reference transcriptome: A high-quality reference transcriptome will be critical for any downstream analysis. Could the authors please provide data for how they validated their reference transcriptome? For example, the authors could report the BUSCO score for assessing the completeness of their assembly, a cumulative density plot of AED scores for annotation quality, and any other metric they choose. Could the authors also report what percent of their single cell reads mapped to the reference?

We thank the reviewers for pointing this out. We have now included metrics about the reference transcriptome that was used for this study. We discuss the implications of using a transcriptome instead of a genome for the downstream analysis and emphasize this point as the main limitation of our analysis. We now include the suggested parameters and report the BUSCO score (s 876-882) and the mapping rate of the reads to the reference in figure 1 —figure supplement 1 and discuss that the low mapping rate is directly linked to the coverage of the reference transcriptome. We agree that a high-quality genome would be better for single cell RNAseq analysis, however transcriptomes provide an alternative if a genome is not available. We are indeed planning in the future to sequencing the Loligo vulgaris genome, given the genome size of many cephalopods this remains still likely a major effort.

2. The scRNAseq data: Could the authors please provide data that would allow readers to assess the quality of their single-cell data? At a minimum, could they please provide violin plots depicting the range of UMIs, genes, and mitochondrial genes per cell?

We now added a figure supplement (figure 1 —figure supplement 1) that shows violin plots for the gene-per-cell counts, number of UMIs and percentage of mitochondrial genes to better assess the quality and variability of the data.

3. Cluster-specific markers: Could the authors please provide dot plots for the top differentially expressed genes for all clusters? In addition to this, the authors have described a set of genes that they did not identify as cell-specific markers (for example, neurotransmission-related ones). Could they please show dot plots for these as well?

We added a supplementary figure with large dot plot that shows the top ten differentially expressed genes in each of the 34 clusters to help the reader identify the genes discussed in this study as well as potential other markers that are not discussed (figure 1 —figure supplement 2). We also provide dot plots for the neuronal markers that were not cluster-specific but for which we performed HCR ISH experiments. These plots are part of a new supplement (figure 2 —figure supplement 8).

4. Trajectory Analysis: Could the authors please explain how they have chosen the cells used in this analysis, the clustering resolution they have used (please provide a dimplot), and which starting point they have selected? In addition to this, could they please analyze the differentially expressed genes along the trajectory to make their claim stronger? Since this data challenges the current state of the field, it would be useful to demonstrate the robustness of their results by verifying it by:

4a. Using different trajectory analysis methods (Monocle or FateID).

This is indeed an important point. We agree with the reviewers that the explanation for the trajectories was not precise enough. We now provide more detailed information about how the trajectories were calculated and which parameters were used in a new section of the material and methods section titles "trajectory inference" (lines 961-981). We replicated the results obtained with slingshot by running the analysis using monocle3. We decided to include the monocle-based trajectories in the main figure (Figure 3) and the results obtained are described in the manuscript (lines 461-489).

We also identified additional genes that can be followed along the trajectory axes (Figure 3 —figure supplement 4). The analysis using slingshot is now part of the supplementary data and we could also follow the expression of several genes along the trajectory on a dot plot (Figure 3 —figure supplement 5, lines 490-501).

5. The writing: Could the authors please pay attention to some aspects of the text to help the reader place their work in the broader context of the field? Specifically:

a. There are missing references (see detailed reviews below), including three other recently posted cephalopod scRNAseq preprints. Could the authors please incorporate these references?

b. On occasion, the chronology of the text and figures don't match. There are also some figures that are currently not mentioned in the text. Could the authors please fix this?

c. Many of the methods are not sufficiently well explained. Could the authors please pay attention to this, particularly for bioinformatics? Could the authors motivate this better?

d. The Results section contains many points that would belong better in the Discussion section.

e. Conversely, the discussion as it currently stands, is a summary of the results. Could the authors please revisit their discussion to place their work and findings in the context of the growing field of cephalopod neurobiology?

We thank the reviewers for their suggestions for improving the writing. We added the missing references for the three recent cephalopod single-cell studies. We also added additional references to better support the rationale and the interpretation of the data. We checked the chronology of the figures and made sure that all figures and supplements are cited in the text. We developed the methods section to include more detail about the downstream analysis of the data and the bioinformatics. As stated above, we now also give more detailed information about the calculation of trajectories for figure 3 (Materials and methods, trajectory inference). We also included the information that was previously missing about the sub-clustering parameters used for figures 2 and 3 in the Materials and methods on lines 944-957. We developed the discussion to better place our findings within the existing literature and to further address possible issues with the analysis. Because the results are very descriptive by nature, we decided to keep some of the discussion-like elements in the Results section to provide the reader with a minimum of context before reaching the discussion. However, all the points mentioned in the results are then discussed more in depth in the discussion.

Review #1 (Recommendations for the authors):

1. I was concerned about the quality of the scRNAseq data for a few reasons:

- The number of genes per cells (200-2000) seems low to me.

- The clusters don't really resolve.

- Many genes that one would have expected to pick up, for example, neurotransmitter-related genes for the neuronal cluster, were not picked up.

So, could the authors please provide the violin plots depicting the range of UMIs, genes, and mitochondrial genes per cell?

Thank you for pointing this out. As discussed in point 1 of the joint reviews, we now provide the metrics for the quality of the RNAseq experiments as well as the reference transcriptome and the details concerning mapping rate of the single-cell reads (Figure 1 —figure supplement 1, lines 132-153). The quality of the metrics of the single cell RNAseq experiments are of high standards. Based on combined metrics of transcriptome and single-cell RNseq, we propose that the lack of some relevant marker genes and neuron-specific markers is caused by coverage of the reference transcriptome. This is likely due to low expression level and transcript abundance of some critical markers, indeed property that differs between transcritpomes and genomes. We now also discuss that the use of a transcriptome versus genome impacts the analysis.

2. It was not clear to me how the authors picked the marker genes. Could they provide the top 20 DEGs for each of the clusters?

Markers were chosen typically from the top differential expressed genes. We now provide a list of the top hundred DEGs for each cluster in supplementary table 1. Additionally, the expression of the top 10 differentially expressed genes for each cluster can now be seen on a dot plot included as figure 1 —figure supplement 2.

3. The signal in each of the images could be enhanced. To orient the reader, it would also be useful to have an inset of each of the expression patterns in the whole head. Some of these are really nicely represented in the supplementary information.

Thank you for raising this point. The contrast was enhanced on all the confocal micrographs. To help the reader understand the pictures shown on the main figures, we provide an inset showing the location of the micrograph in the animal as well as anatomical annotations, similarly to the figure supplements. We additionally include a description of the anatomy of the Loligo nervous system in the Results section on lines 253-266 with a supporting figure (figure 2 —figure supplement 1) to help orient the reader. When available, the overviews of the tiled confocal micrograph showing the whole head are available in the figure supplements.

4. The authors should refer to the three recent preprints that have used scRNAseq in other cephalopods and, in the discussion, place their work in this context.

We now cited these three papers and have placed out analysis in the context of these studies (Lines 98, 656 and 679).

5. In general, could the authors elaborate their method sections? Many fairly important aspects are too briefly described – for example, the reference transcriptome and its validation and the selection of marker genes, among others.

We indeed agree with this point. We have further and extended developed the results and Materials and methods section to include more details about the transcriptome quality assessment as well as the mapping rate of the reads to the transcriptome (Results: lines 132-153). We added a more detailed description of the clustering and sub-clustering process done in Seurat and we wrote a new section to describe the way the trajectories from figure 3 were obtained (Materials and methods: lines 961-981).

Reviewer #2 (Recommendations for the authors):

The core of this study is the analysis of scRNAseq data. Here the authors assembled a transcriptome de novo. I wonder about the quality of this assembly, as it will affect all downstream analysis. The other 3 studies (Styfhals et al. bioRxiv 2022.01.24.477459; doi: https://doi.org/10.1101/2022.01.24.477459; Songco-Casey et al. bioRxiv 2022.06.11.495763; doi: https://doi.org/10.1101/2022.06.11.495763; Gavriouchkina et al., bioRxiv 2022.05.26.490366; doi: https://doi.org/10.1101/2022.05.26.490366) all refer to the necessity of a high-quality reference genome for a high-quality transcriptome. If this is not the case here, it would be useful to quantify this.

I found the analysis of stem cells and progenitors quite interesting. The authors profile specific markers for cell types, map these cell types onto interesting anatomical regions, and analyze developmental trajectories defined transcriptomically.

In other parts of the paper, I was surprised that many cell types could not be defined based on certain marker genes/small combinations of genes. Figure 2b, for example, shows that clusters have a great deal of overlapping expression among 'marker' genes. This would deviate quite a bit from the other scRNAseq studies I am familiar with. I could see a few possible explanations: (1) This is a product of low-quality transcriptome assembly/single-cell data (2) This is a product of the immature developmental stage of the animal, and (3) This reflects a cephalopod-specific transcriptomic signature.

I doubt possibility 3, as the other 3 cephalopod scRNAseq papers describe cell types well differentiated (eg. By neurotransmitter expression, Figure 2 in Songco-Casey et al. 2022, Figure 3 in Gavriouchkina et al. 2022). The fact that specific neurotransmitters are often not well localized to specific cell types here calls for an explanation. Similarly, the lack of differentiation of neural types even when analyzed separately from other cell types (Figure 2c). Many phototransduction genes are expressed in Cluster 32, and these genes are expressed, specifically in the retina. They could potentially be used as a sanity test for the scRNAseq: Are these genes only found in cluster 32?

We thank the reviewer for these valuable comments. In light of the documented impact of reference quality on the downstream analysis, we now provide additional information regarding the quality of the single-cell reads and the reference transcriptome (Results: lines 132-153, figure 1 —figure supplement 1). As mentioned above (joint reviewer comments) we indeed fully agree that a reference genome would be beneficial for the downstream analysis of sc-RNAseq data. We aim to take this on the task of generating a Loligo vulgaris genome in the future. With the transcriptome as reference, we cannot exclude that some cell types do not resolve because of the developmental stage of the animal. The analysis of the completeness of the transcriptome as measured by its BUSCO score, it seems a likely explanation for the lack of cluster resolution. This is probably most relevant for the absence of some genes such as makers of neurotransmitter synthesis pathway or patterning transcription factors. However, we feel that the current analysis provides a substantial advance, particularly in an otherwise molecularly underexplored taxon. With the reference markers we can identify and characterizes most clusters. Most importantly we feel that the histological correlation and mapping of cell types in situ is powerful. When looking at the example of phototransduction genes, they are very specifically expressed in cluster 32 but they can also be found in other cell types at much lower expression levels. The cell-type specificity observed with HCR ISH in the retina may not be the full picture. As in most other species one might expect these genes to be also present in other cells and other tissues, as exemplified in recent years in model organisms such as Drosophila or mouse.

The lack of specific marker genes makes some of the FISH experiments in this manuscript less informative than others. Many plots show anatomical localization of a gene broadly expressed in many cell types, meaning we cannot conclude how transcriptomically identified clusters map spatially. In other cases, FISH is used to mark individual cell types. In some cases, it is not clear which situation we are looking at. It would be useful to have a quantification of how specifically the marker genes used in FISH experiments mark individual clusters.

Thank you for raising this important point. We now provide additional information regarding the cell type specificity of markers of interest in the supplementary table 1 and in figure 1 —figure supplement 2. Generally, the choice of marker genes was more based on the availability of information and literature regarding their function or cell type specificity rather than a purely mathematical selection. This explains why some of the genes picked in our analysis are not always the markers with the lowest p-value but are the ones that have been previously associated with cell types of interest, development and, in certain cases, that are candidates for marking novel cell types as exemplified by the sections about papilin and madnf expression. It is indeed an intriguing idea to start probing for instance cephalopod or Loligo specific genes in the future.

The authors observe a similar anatomical expression of fmrf and reflectin-8 in the olfactory organ. Because these have been linked to chromatophore and iridiophore control, they argue that this suggests that chromatophores and iridiophores are modulated through the olfactory organs. I find this plausible but the logic questionable. FMRF has many functions unrelated to chromatophore modulation, so just observing its presence seems insufficient to make statements about the neural control of chromatophores. Reflectin-8 is shown later in the manuscript to be broadly expressed in the epidermis and brain (Figure 5c).

This is indeed important and we fully agree with your point. We apologize for this overstatement. We have now reformulated this sentence. It now reads: " FMRFamide has been previously shown to be involved in the control of chromatophores in squids (Loi and Tublitz 2000; Sweedler et al. 2000; Satpute-Krishnan et al. 2006; Mobley et al. 2008)".

I was not clear why the expression of elav-like 1 in the epidermal lines 'informs on their potential mechanosensory function'. With elav-like 1 expressed very broadly across cell types and anatomical locations.

Several genes that we analyzed are highly expressed in some clusters but are also expressed at lower levels in other tissues. As discussed above with the example of phototransduction genes, HCR experiments tend to only reveal cell-types in which the gene is highly expressed. For this reason, we interpreted the high signal from the epidermal lines as an indication of the presence of neurons in this structure. The mechanosensory modality of these neurons is based on studies done in other cephalopods. We reformulated the sentence to make this clearer.

It now reads: " The expression of elav-like 1 in this structure in L. vulgaris is consistent with their proposed mechanosensory function other cephalopods."

The differential expression of the two elav genes is quite interesting and might reflect different functions as recent results in Drosophila of the three Elav-family members (Elav, Fne and Rbp9) suggest (Lee et al., 2021,Plos Genetics).

The expression of the cilia-associated protein, rab3 and synaptogyrin-1, and other cilia/flagella markers in cluster 27 suggested that this corresponds to a mechanosensory population. To test whether this is neural, the authors take the subset of cells in the dataset expressing amyloid β- binding-like (a neural marker) and see whether a subset of these also expresses the collection of sensory markers. This seems round about to me. Why not look directly at whether the cells in cluster 27 express amyloid β- binding-like? If the probability of co-expression is high, then concluding that these cells are neurons seems warranted. If the probability is lower, then the conclusion does not seem warranted. The author's existing analysis does not address the possibility that non-neurons (those not expressing amyloid β- binding-like) express the collection of sensory markers.

We thank the reviewer for pointing out this issue. We indeed realized that this is a circular argument in this case as cells were originally identified as neurons on the basis of the expression of markers such as amyloid-β binding. The intent was to visually represent the co-expression of these newly identified sensory markers with neuronal markers. We changed figure 4 to plot the expression of the sensory genes in the set of neuronal genes used in figure 2. We reformulated the sentence to highlight that this is not a proof of their neuronal identify but rather a demonstration that the hypothesis regarding sensory cell types is consistent with the hypotheses made in the previous sections. It now reads: " The results highlight the co-expression of neuronal markers and sensory markers, therefore supporting their sensory identity (figure 4D)." (lines 550-553).

Line 63: 'Coleoid cephalopods, which have internalized or lost their shell over evolutionary time 64 were able to swim more easily and could exploit other niches.'

Not clear what 'other' refers to. The niches of ammonoids? Not clear whether it was the loss of shell that led to niche exploitation or the other way.

We now specify that coleoid cephalopods were able to access resources in open waters, which was not the case for shelled cephalopods. It is not possible to confidently state that one event caused the other in this case. However, the exploitation of different niches, which was only possible because of the increased mobility of coleoid cephalopods, favored the differences in nervous system anatomy that can be observed today between Nautilus and coleoid cephalopods.

Line 93 "… as a first step to understanding cephalopod evolution". I'm not clear why a molecular comparison would be a first step to understanding cephalopod evolution. What then is the non-molecular evolutionary story given earlier in the introduction?

We agree that this statement is not accurate. We reformulated this sentence to state instead: "These comparisons constitute a major step for understanding the molecular basis of cephalopod evolution" (lines 92-93).

Figure 2 Are the Umap plots in 2c re-estimated on the subset of data in A? If so, I don't see a description of this in the methods section, and I am surprised at the lack of cell type differentiation among neurons.

The clusters of figure 2A were indeed re-clustered to obtain the plots of figure 2C. We have now included a description of this process and the parameters used for subclustering in the material and methods (lines 944-950). Additionally, the UMAP with the new clusters is available in figure 2 —figure supplement 8.

Figure 4D is described as 'feature plots'. Is this UMAP? How are the plots on the right-hand side generated?

We modified figure 4D and the plots shown there are UMAPs representing gene expression in the subset of cells as for figure 2 (figure 2 —figure supplement 8, Materials and methods: lines 944-948).

The developmental trajectory analysis using slingshot could be better explained in text/methods/figure legend and not easily interpretable from Figure 3e if one is not familiar with the technique

We have included a section in the material and methods section to explain how these results were obtained (Materials and methods: lines 961-981). In addition, results were replicated using a different trajectory inference method monocle (Results: lines 461489, Materials and methods: lines 961-976). The monocle data is now part of the main figure (figure 3) and the slingshot data obtained previously was moved to the supplements (figure 3 —figure supplement 5). More data related to the trajectories using monocle and slingshot are displayed on figure 3 —figure supplement 4 and figure 3 —figure supplement 5 respectively.

'Additionally, the cephalopod-specific chromatophores, which are colored or iridescent, have been suggested to be modulated by the nervous system'. Maybe I am misunderstanding the sentence, but as I read it this is an understatement of what we know. Over 50 years of work have profiled the neural control of chromatophores in detail.

We reformulated this sentence to acknowledge the existing knowledge on the subject. It now reads: " Additionally, the cephalopod-specific chromatophores , which are colored or iridescent, have been shown to be modulated by the nervous system (Dubas et al. 1986; Wardill et al. 2012; Gonzalez-Bellido et al. 2014)" (lines 559-562).

Some of the supplementary figures are quite dark. Why not increase the brightness of the images? Others like Figure 4 supplement 1 are beautiful.

We have now increased the contrast on many images to help the reader better visualize the expression patterns shown in the figures.

Reviewer #3 (Recommendations for the authors):

First of all, I would like to congratulate the authors on the extensive expression analysis study performed in this manuscript. However, I believe that this dataset offers so many more opportunities to better understand cell types in cephalopods, and in general, the data could be analysed and interpreted in more detail.

Since the field of scRNAseq is moving incredibly fast, the authors should be careful with claims of being first (e.g. Abstract line 34). Preprints of scRNAseq data in the nervous system of other cephalopods have been recently published (Styfhals 2022, Gavriouchkina 2022, Songco-Casey 2022).

Thank you for mentioning this. We have now cited the recent single-cell RNAseq preprints that were missing from the manuscript (Lines 98, 656 and 679) and have removed the claims of being the first analysis of this kind.

The manuscript is missing key citations to support the results and to provide a comparative view of gene expression in other cephalopods. For example, the expression patterns of Elav, SoxB1, Asc have been reported during the development of octopus, squid, and cuttlefish, and could enhance the interpretation of the results observed in Loligo.

We added important citations related to the expression and function of these genes in other cephalopods as well as in other clades. We have developed the discussion to better confront the results obtained in our study of Loligo vulgaris to the existing literature (discussion: lines 682 – 733).

In general, the reader needs to be proficient in the brain areas of cephalopods in order to follow the results. Therefore, the manuscript would benefit from a short introduction of the main brain lobes in cephalopods to be accessible to a wider audience, either in the introduction or in the results in Figure 1C.

This is indeed a good point. An introduction to the brain areas that are mentioned in the manuscript is now available at the paragraph "molecular mapping of the L. vulgaris nervous system" (lines 253 – 266). To better understand the spatial arrangements of the lobes, we also provide a rough schematic of the anatomy in a new figure supplement (Figure 2 —figure supplement 1).

Please make sure to refer to the figures and figure panels chronologically. Multiple figure panels are not referred to throughout the text (e.g. the first figure described is Figure 1B and there is no notion to Figure 1A, similar to Figure 2, …).

Thank you for pointing this out. We have now checked the chronology of the figures, so that they appear in the right order in the text.

The authors distinguish 33 cell clusters when using 25 pca's and a resolution of 2. The UMAP, however, shows most cells in one single, central blob. This is striking knowing how different connective tissue is from neurons. Did the authors try other methods for clustering and dimensionality reduction, or change the parameters in the current algorithm to try to resolve the data better? It would also be helpful to specify the filtering parameters for filtering out low-quality cells in the methods section.

We have indeed noticed that the resolution of the clusters using different tissues (such as in the head) is quite different from resolving similar cells in the brain, which we have mainly focused on in the past. We have initially tested several values for PCs and resolution and felt 25 PCs and resolution 2 gave the best result. This issue of having clusters too close to each other could not be resolved using different clustering parameters (PCA selection method, number of PCs, cluster resolution). Representing the data with TSNE instead of UMAP also showed similar issues. We propose in the discussion that such lack of resolution, although it is sufficient to analyze cell types, is due to missing genes in the reference transcriptome (Results: lines 367-372, discussion: lines 652-666, 680-682). It is actually an intriguing point we frequently discuss which way might be more adequate: splitting or grouping clusters; this may be partly a semantic issue or a matter of taste and mostly decide on a middle path.

In the first section of the results, the authors list a handful of enriched genes used to differentiate the major cell type categories. However, little explanation is given to why these genes were chosen and a comprehensive overview supporting the claims is missing. It would help the reader if they would add a plot (e.g. a DotPlot) to support these data. In addition, a UMAP plot with cells colored according to cell cycle could clearly show where the progenitor cells are located on the UMAP.

This is indeed a good point. As mentioned above, the genes were chosen among a list of differentially expressed genes in the clusters. The choice among the top markers is made based on previous knowledge on the association of the gene to a specific cell type. For example, genes of unknown function that have never been studied in other species were not selected in our analysis because they could not inform us on the identity of cell types, even though they may be statistically relevant. The choice of markers was generally subjective based on our own interests and pre-existing knowledge as the systematic analysis of all the top markers would have been practically impossible in our case. To provide the reader with a better idea of the relevance of the genes chosen in this study, we now provide a dot plot showing the top markers across all cluster (Figure 1 —figure supplement 2) in addition to the full list of differentially expressed genes (supplementary table 1).

To avoid going back and forth between cell type categories in the first results paragraph, it would be helpful to describe Cl 27 (lines 202 onwards) and Cl 32 (lines 207 onwards) with the other nervous system cells, and Cl 29 (lines 213 onwards) with the epidermis section.

We now grouped the description of Cl27 and Cl32 with the description of neuronal markers in this section and the description of Cl29 directly following the description of the epidermis markers (results: lines 198-205).

In order to interpret HCR expression patterns, the authors should add axis labels to the figures. In addition, I understand that the figure panels are small to add annotation of the brain regions, but since we are looking at different levels through the brain, either a schematic with the figure or annotation would really help.

We agree that the pictures did not provide enough information for the reader to understand the exact location of the picture. In addition to the insets, we now also added axes labels and annotation of the anatomical structures visible on the micrographs.

The authors should explain the analysis performed to obtain the UMAP(?) plots in Figure 2C. Which cells were selected? Were they re-clustered? What do the clusters look like?

The UMAP plots on Figure 2C are obtained from a re-clustering of the cells highlighted in Figure 2A. We know added information about this in results and details about the parameters used to obtain this was now added in the materials and methods (lines 944957). The clusters obtained with this method are also shown with a UMAP on figure 3 —figure supplement 8.

Did the authors perform image manipulation e.g. Figure 2C and 2D for the DAPI channel? Why is DAPI absent from the neurons in the tetraspanin-8 HCR, or the serotonin transporter staining?

The images were not manipulated but this is the result of an experimental issue. All the HCR ISH experiments were performed on whole mount samples and sometimes we encountered problems with the penetration of DAPI inside the denser tissues. This is visible on several pictures where the outside of the animals is strongly labeled with DAPI but the signal fades as we go deeper in the tissue. We realize that this may confuse the reader and added a word of caution regarding this problem in the results and the concerned figure legend (lines 306-309, lines 1075-1077).

In Figure 3B, it would be helpful to add some genes involved in proliferation, as exemplified in the text.

We added new plots showing the expression of additional markers in the scRNAseq along the proposed differentiation axis. The genes in question are pcna, elf2-like, histone H3, ets4, hes-4, chymotrypsin, synaptogyrin, elav-like1, elav-like 2 and reflectin-8. These plots are available in Figure 3 -supplements 4 and 5.

The authors claim that neural stem cells undergo neuronal differentiation and then migrate to the nervous system. However, evidence is missing as the authors do not show where progenitor cells stop dividing and start to differentiate.

We agree that currently there is no data supporting this hypothesis. Although we show that cells in different stages are present in different locations, the timeline of differentiation and migration cannot be inferred from our data but we propose that such a stage exists and is yet to be characterized. We reformulated these statements to be more representative of the evidence. It reads as follows: " Together, the overlapping expression patterns of e2f3, asc and soxB1 in domains associated with cephalopod neurogenesis strongly suggest that stem cells and neuronal progenitors are present the lateral lips of L. vulgaris suggesting that these cells are likely migrating to the nervous system at later differentiation stages " (results: lines 423-427 ).

The use of the marker foxD1 for differentiation neurons is not well supported. Additional marker analysis is required before this can be claimed.

This is indeed an important point, we apologize for not being clearer on this. In fact the, choice to include foxD1 in the analysis of neuronal differentiation was made after the HCR ISH experiments that showed its expression in the nervous system. In correlation with literature that link foxD1 genes with neuronal differentiation, we decided to check if we could computationally link the foxD+ clusters with neuronal clusters. We are aware that foxD1 genes are not a typical feature of differentiating neurons as it is also broadly expressed in other tissues in the development of other species (Herrera et al. 2004; Polevoy et al. 2017; Newman et al. 2018; Janssen et al. 2022). In the section about trajectory inference, we added additional markers of these cells (Figure 3 —figure supplements 4 and 5). We generally reformulated interpretations regarding foxD1 to avoid overstatements. While we are confident about this observation of where foxD1 is expressed in situ and in silico, we indeed agree that with the current data we may not define it as neuronal differentiation and have widely re-written this section and that it will in the future need to be thoroughly tested (discussion: lines 733-753).

The authors need to explain which clusters were used for trajectory analysis using Slingshot and how the data were reclustered. Were non-neuronal stem cells included? This might skew the analysis. Importantly, the method was not described in the methods section. It is critical for the readers to know which endpoints were selected in the slingshot algorithm. Off note, the slingshot is very dependent on clustering resolution and will try to connect all clusters. Have the authors tried a different pseudotime analysis method such as Monocle or FateID? Additionally, can the authors elaborate on which transcription factors are steering this differentiation? In the current state, the analysis is preliminary and contradicts present data on the origin and trajectory of neurons in the cephalopod brain (see comment below; Koenig et al. 2016, Deryckere et al. 2021).

Thank you for pointing out the lack of information regarding the trajectories and the need to replicate these results. We now include more detail about how the trajectories were calculated and which parameters were selected in the materials and methods in a new paragraph called "trajectory inference". We also include details about how the cells were chosen and re-clustered for figure 2 and 3 (Materials and methods: lines 944950). We now include in the main figure the analysis of the trajectories using monocle in addition to the previous analysis with slingshot that has been move to the supplements (figure 3 —figure supplement 5). The other cluster specific transcription factors that could be identified were also plotted onto the trajectory analysis on figure 3 —figure supplement 4 and 5. However, many genes that have been linked to neuronal differentiation (for example neuroD in Deryckere et al. 2021) could not be found in our data and are therefore not available to support our hypotheses. More generally, we agree that these results are preliminary and require further experimental validation (see discussion).

With the current interpretation of the data, the authors claim that foxd1 cells are an intermediate developmental state of neurons. This would be very striking since the cortex of the optic lobe has been described as being the inner retina. Additionally, it is contradicting the expression of amyloid β-binding like in the optic lobe cortex, which the authors say is only expressed in differentiated neurons. Is amyloid β-binding-like expressed in clusters 15, 19, and 20? How about other markers for differentiated neurons? What would all those intermediate cells do there? Do the authors think they migrate into the medulla or even the central brain? This is contradicting existing literature on the origin of neurons in the cephalopod brain (Koenig 2016, Deryckere 2021) and requires in-depth discussion.

We indeed agree on this point and as stated above we have included a more moderate interpretation. Our analysis of foxD1 originally showed two different things:

1. This gene is very specific to certain cell clusters and is highly expressed in these cells compared to other clusters (Figure 3A). These cells cannot be directly assigned a cell-type identity as they lack known cell-types markers, including markers of differentiated neurons.

2. This gene is highly expressed in different regions of the nervous system, including optic lobe cortex (Figure 3, figure 3 —figure supplement 2). Importantly, we do not think that the expression of this gene in these regions necessarily means that they are differentiated neurons.

Because these two observations do not clearly indicate a link between foxD1 and neurons, we decided to computationally assess whether there might be a developmental relationship between the cell types found in the nervous system and the progenitor-like cells present in the lateral lips. Contrary to our expectation, the trajectory analysis put foxD1 cells in between putative progenitors (expressing e2f3, asc and soxB1) and differentiated neurons (expressing elav-like 1, elav-like 2, tetraspanin-8, amyloid-β-binding) (We do find some cells expressing low levels of amyloid β-binding-like in FoxD expressing cells (clusters 16, 19, and 30)). Essentially, it means that there are more transcriptional similarities between foxD1 and neurons (respectively between foxD1 cells and asc+ cells) than there are between asc+ cells and neurons. Whether that means that foxD1 expression marks a specific stage neuronal differentiation is not clear and will needs to be further investigated. We moderated our conclusions regarding the analysis of foxD+ in the text and provide more discussion of the subject (lines 733-753). We agree that a theoretical intermediate cell stage could in any case not be defined by the expression of a single gene and that thorough analysis of more transcription factors would be necessary to support these claims. The absence of many neurogenesis-associated transcription factors from our data renders this analysis more difficult and a higher resolution dataset may be able to conciliate our findings with previous literature.

The authors prepared beautiful summary diagrams in Figures 3 and 4. However, they should be referred to in the main text.

Thank you for pointing this out, we have now made sure to refer to the diagrams in the text (lines 460, 553)

After performing extensive HCR expression studies, can the authors go back to their scRNAseq data and better annotate the different clusters?

In order to better understand neuronal diversity and neural cell types (as the authors claim), the dataset has to be analyzed in more depth. 33 clusters/cell types seem low given the number of cells that were sequenced and the different tissues that were sampled. Can the authors identify the different brain regions when subclustering the neurons? Can they identify and annotate different neural cell types?

This is an important point and we indeed hoped to find genes identifying neural subtypes or defining certain brain regions. We indeed identified some genes with regionalized expression, such as scratch and LIM, which to some degree can be associated to certain clusters (scratch 5, 15, 20; Lim high in 22). However, we did not uncover wide diversity of neuron subtypes. This may be partly due to the fact that we took entire heads for our analysis (as compared to only brain) and partly to the coverage of the transcriptome.

The HCR ISH studies provided a lot of additional information about the genes originally described in the scRNAseq data and the cell types that they represent. The interpretation of the results combining these two methods provides avenues for future studies but is usually not sufficient to confirm the identity of a previously undescribed cell type. The clustering parameters used initially resulting in a total of 34 clusters could have been changed to obtain more clusters but the markers genes of these new clusters could not have informed us on the possible identity of these cell types. The example of the lack of resolution for the nervous system well exemplifies this limitation in our study: several genes were shown to be expressed in specific lobes of the nervous system but this did not correlate with specific clusters in the scRNAseq data. We envision in the future to extend our investigations by including more developmental stages (and behavioral context) and to possibly assembling a Loligo genome to map the same reads to may resolve this issue (Results: lines 369-372, discussion: lines 652-666, 680-682).

The authors present several striking findings. However, at the moment, the Discussion section explains the results obtained in the paper, but does not place the acquired data within the current state-of-the-art. Parts of the results should be moved to the discussion and the discussion should be elaborated on.

We rewrote the discussion to better contextualize our findings within the existing knowledge about cephalopod cell types and the different functions of the genes discussed in this paper. We also cited many additional papers to better support our interpretation of the data.

Specific comments:

Line 76: "by the use its specialized iridescent cells" misses an of

Line 97: it would be helpful for the readers to add the "pre-hatchling stage" to stage 28, as was done in the Results section for the audience to have a better idea of what stage was studied.

We now made sure to specify pre-hatchling in every mention of the specific stages (lines 102, 133, 432, 835, 987).

Throughout the text (starting line 219) please write "HCR ISH" instead of just "HCR". The chain reaction method is not specific to in situ hybridization and can also be used in immunohistochemistry and others.

Line 101: "in all animals" please be more specific; not all animals have neurons.

We corrected this sentence to "in most bilaterians".

Following the recommendation of the reviewer, we now refer to the experiments as HCR ISH experiments.

Line 110: "marking distinct stages of neural development", neural differentiation would be more appropriate.

We corrected this and changed to "neural differentiation".

Line 126: Add a reference to Figure 1A to show what the head region is.

The dotted line on figure 1A represents the plane of amputation. Everything anterior to that was used for these experiments. We made that clearer in the legend of figure 1A: " The dotted line represents the plane of amputation. The region anterior to this line is what was used in this study and is referred to as the head region."

Line 133 says 34 cell clusters were identified, but Figure 1B shows 33.

We kept the default counting system used by Seurat where the count starts at 0. That is why there are 34 clusters but they are numbered to 33. This is specified in the text as follows: " This resulted in a total of 34 cell clusters numbered from 0 to 33" (lines148149).

Line 145: "smooth muscle-like cells (Cl26)" should be Cl25.

We thank the reviewer for noticing this mistake, the number was changed to 25 (now on line 159).

Line 156-159: add a reference for this statement, maybe a specific example.

We added the following citations to the statement referring to previous studies of connective tissue in cephalopods: Thompson and Kier 2001, Kier and Stella 2007.

Line 167-171: "a gene whose function is homologous to …" and " exact function … is not known". This is contradictory. Do you mean that the gene's sequence is homologous?

This sentence was reformulated to explain that the function of the Drosophila melanogaster papilin ortholog has been studied but the function in other species is unknown. It reads as follows: "… papilin, a gene which was suggested to play a role in the formation of extracellular matrix in Drosophila melanogaster (Fessler et al. 2004; Apte 2020). The function of this gene in other clades is however not known" (lines 185188).

Line 199-201: with the comments below on SoxB1 and Ascl1, this sentence needs to be revisited.

This sentence was revisited according to the discussed roles and expression patterns of soxB1 and asc: " The presence of a soxB1 in these cells indicates that they may proliferating neuronal progenitors, based on previous evidence of expression of soxB1 genes during neurogenesis in cephalopods (Miyagi et al. 2009; Focareta and Cole 2016)" (line 227-230).

Line 239-240 and also further in the text: "in the medial portion of the optic lobes" do the authors mean the medulla (in contrast to the cortex (inner and outer granular layer))? Medial portion is not clearly defined and the cortex also stretches quite medially. (also line 278, 540).

We did indeed mean the medulla when mentioning the medial portion of the optic lobes. We now specifically refer to the medulla and cortex in the text to use the correct terminology and avoid unnecessary confusion.

Lines 263 and 271: please rearrange the supplementary figure so you can refer to panels A-D first. In general, please refer to figures chronologically.

The figure was rearranged to match the chronology of the text.

Line 298 with Figure 2D: please switch panels of LIM and serotonin, so they follow the text.

The panels were switched in order to match the text (figure 4).

Line 326: this part covers more than stem cells and progenitors (line 375 and onwards).

This section was renamed stem cells and neuronal differentiation to include the proposed cell types of the neuronal lineage.

Line 327: please add an introductory sentence to this new paragraph; rephrase the first sentence (grammatically incorrect).

We reformulated this sentence to " Several marker genes expressed in specific clusters suggest an involvement in proliferation and possess some conserved markers of neuronal differentiation". We also added an introduction to the paragraph (line 376382).

Line 338-340: this is an overstatement. What identity are the authors referring to? Genes involved in proliferation etc do not give identity to stem cells and Cl02 and Cl14 are also actively cycling. Do the authors suggest only Cl10 and Cl17 are stem cells? Can the authors show more evidence?

This sentence was modified and references were added to support the claim. We indicate that cells of Cl02, Cl10, Cl14 and Cl17 are all proliferating cells. E2f3 is highly expressed in Cl10 and Cl17 and is also co-localized with asc in the lateral lips. The sentences were revisited as follows: " The role of e2f3 to promote cellular proliferation has been shown to be conserved between D. melanogaster, C. elegans and mammalians (Dynlacht et al. 1994; Attwooll et al. 2004; DeGregori and Johnson 2006). Because these clusters highly express proliferation markers but lack any distinctive features of differentiated cell types, we propose that they may be stem cells" (lines 392-397).

Lines 350-359: I would suggest moving this part to the discussion and incorporating the following comments:

Lines 347-352: The expression of Asc is not as well conserved as the authors suggest. Indeed, it has some proneuronal role, but the level of differentiation in the cell it is expressed is very different across the animal tree. In Drosophila as in other non-insect arthropods, Asc is expressed in quiescent ectodermal cells to drive them into the neural lineage (formation of a neural progenitor cell) (e.g. Cabrera et al., 1987; Skeath and Carroll, 1992). In contrast, in vertebrates and also in spiralia such as capitella teleta, Ascl1 is expressed in neural progenitor cells to drive them into differentiation (e.g. Bertrand 2002; Sur 2020). Ascl1 expression has been studied in octopus and seems to have a similar function as in capitella, not Drosophila (Deryckere 2021). This scRNAseq dataset has the power to test this hypothesis and figure out whether, in cephalopods, as is expressed in proliferating and not quiescent cells.

Lines 354-355: this conclusion does not reflect the data shown. The authors are showing proliferating stem cells, not the location of differentiation. In octopus, differentiation seems to happen in the transition zones (NeuroD cells, Deryckere 2021). It would be interesting to locate those cells in your dataset.

Lines 356-359: Again here, one needs to be careful to extrapolate the functionality of a gene in one species to other species. In invertebrates, SoxB1 is also expressed in differentiated neurons (e.g. Le Gouar et al., 2004; Semmler et al., 2010). Additionally, expression of SoxB1 has been described in Sepia (Focareta and Cole, 2016) and octopus (Deryckere 2021), indicating expression in neurons (with some interesting discrepancies). Furthermore, be careful with the use of the word "neuroblast" because it refers to different types of cells in different animals.

Line 361: "Interestingly" it seems a general theme for cephalopods, and maybe invertebrates in general (see comment above).

Line 367: Data are missing to show where the neural progenitor cells differentiate.

Line 370-374: in line with this, SoxB1 expression is also found in sensory epithelia of sepia and octopus (e.g. suckers and statocyst; Baratte, S. and Bonnaud 2009; Deryckere 2021).

We thank the reviewer for these very helpful comments. We know discuss the discrepancies in the role of asc across different phyla and can replace our findings in this context (lines 709-721):

"We interestingly observed the expression of asc in the same region, together with several proliferation markers. Extensive study of neurogenesis in Drosophila melanogaster has shown that the proneural complex asc can determine the neuronal fate of quiescent multipotent ectodermal cells (Campuzano and Modolell 1992; Skeath and Doe 1996). In vertebrates on the other hand, asc is expressed in cells that have already undergone neuronal specification but are self-renewing (Bertrand et al. 2002; Soares et al. 2022). Studies of neurogenesis in the annelid Capitella teleta has shown that, similarly to vetrebrates, proneural genes such as asc are expressed in self-renewing neuronal progenitors (Sur et al. 2020). Our data suggests that asc is also expressed in proliferating cells, similiarly to vertrebrates and Capitella teleta but unlike Drosophila.”

We also discuss the role of soxB1 and compare it to what has been observed in other cephalopods but also refer to the documented expression of this gene in differentiated neurons (lines 721-733):

“The observation of soxB1 expression together with asc strengthened their characterization as neuronal progenitors due to the documented role of this gene in neuronal differentiation (Miyagi et al. 2009). However, strong soxB1 expression in differentiated neurons has been shown in the gastropod Patella vulgate and in the cephalopods Sepia officinalis and Octopus vulgaris (Le Gouar et al. 2004; Focareta and Cole 2016; Deryckere et al. 2021). In L. vulgaris we observe similar results where soxB1 is expressed in the lateral lips together with asc but is also expressed in differentiated neurons. This was consistently shown by the HCR ISH experiments as well as the single-cell analysis. Therefore, the hypothesis proposed for the action of soxB1 in the early stages as well as the later stages of neuronal differentiation in Octopus vulgaris is consistent with our data (Deryckere et al. 2021).”

We also discuss the possible limitation of our analysis as some genes that have been shown to be involved in Octopus vulgaris neurogenesis could not be found in our data (lines 742-753):

“Such a role has been proposed in Octopus vulgaris in the cells expressing neuroD that are expressed in a gradient between the lateral lip and the central nervous system, which may indicate a pattern of migration of differentiating neurons (Deryckere et al. 2021). We could however not correlate these results with our hypothesis of an intermediate foxD1-expressing stage because neuroD expression could not be detected in specific clusters in our single-cell dataset. Although the trajectories we calculated point towards a role of FoxD1 in neuronal differentiation, the relatively poor cluster resolution and the lack of cells along the trajectory axes shown in figure 3 suggest that this hypothesis would require further investigation that should include a more complete set of genes that have been involved with cephalopod neurogenesis.”

A short interpretation of the data is left in the results to facilitate comprehension of the data but is then discussed more in depth in the discussion.

Line 376-378: add a reference and elaborate. Regulators of development is a very general process. In addition "for this reason, these cells were presumed to be differentiating cells": this is an overstatement and needs additional gene expression patterns to support this presumption.

We corrected this sentence to not overstate the identity of these cells. Although there is evidence of foxD1 being involved in the development of the nervous system, the hypothesis of these cells being a possible stage of neuronal differentiation is only prompted by the expression pattern observed by HCR ISH followed by the trajectory analysis described later in the manuscript. The sentence was modified as follows: " FoxD genes have been documented to be expressed in the nervous systems of Mus musculus, Xenopus laevi, C.elegans and D. melanogaster and to play a role in the development of neuronal cell types although not exclusively (Herrera et al. 2004; Polevoy et al. 2017; Newman et al. 2018; Janssen et al. 2022)" (lines 437-441).

Line 382-385: optic lobe cortex = inner and outer granular layers, please adjust and rephrase. From Figure 3C upper panel, it seems that Foxg1 is also expressed in the medulla. Or is this all cortex? Annotation would be helpful here.

This sentence was revisited to describe the expression pattern that is indeed observed strongly in the cortex but also in the medulla. The sentence reads al follows: " HCR ISH assessment of the expression of foxD1 showed that these cells are strongly expressed in the cortex of the optic lobes and in the medial basal and anterior basal lobes (figure 3C). The expression of foxD1 is very strong in the cortex of the optic lobes and in the medial and anterior basal lobes but is also sparsely located in the medulla (figure 3 —figure supplement 2)" (lines 441- 446).

Line 386: this is an overstatement. Most neurons in the brain express elav1 and the function of foxD1 is very broad. Additional markers (like for example NeuroD) are required to support this claim.

This sentence was revisited as follows: " The expression pattern of foxD1 indicated that these cell types are present in the nervous system and could therefore be types of neurons or neuronal progenitors" (lines 446-448). The claims about the role of foxD1 are stated later in the manuscript and are discussed in a critical manner (lines 733753).

Lines 394-400: Figure 3D is not mentioned in the text and seems substantial. The authors might also want to revise it with the comments above and in public review.

We modified figure 3D to better represent that soxB1 is expressed in both the lateral lips and in the nervous system together with markers of differentiated neurons. The HCR ISH experiments do not by themselves support the presence of an "in-between" stage as was suggested by the figure previously.

Lines 408-411: please add references to support these statements.

We have added references that support the conservation of the phototransduction cascade between cephalopods and other protostomes. The sentence reads as follows: " The phototransduction cascade that allows light stimuli to be transformed into an electrical signal in the neurons is well conserved among cephalopods and we were able to clearly identify photoreceptors in the scRNAseq data (Davies et al. 1996; Monk et al. 1996; Arendt 2003)" (lines 510-514).

Lines 419-421: Why did the authors choose to map the expression of Pax9? Koenig et al. studied the eye in Doryteuthis and mapped the expression of several genes important for eye development. A reference to this work would be appropriate here. Pax9 is not an eye gene in vertebrates and invertebrates (in contrast to Pax2/6). In which clusters is it expressed (again here, a dotplot of the genes under study would be helpful)? The level of the image in Figure 4A for Pax9 is also different than the others, why is this? In which cells of the eye is Pax9 expressed?

Pax9 was remaining for an initial screen of transcription factors that may be related to developmental processes. After careful re-analysis we decided to exclude this gene from the analysis as it is not a cluster-specific marker in the scRNAseq data and the expression pattern revealed by HCR ISH showed expression in what appeared to be non-neuronal tissue around the eye. Because there was no pre-existing evidence of this gene being involved in sensory cells, we realized that including of this gene in our analysis was not justifiable.

Line 448: Figure 5B should be Figure 4B I think.

This is indeed the case. We changed it to 4B.

Line 457: can the authors verify that they selected cells with expression level >3? From the legend in Figure 4D, the expression level of 3 seems present in a lot of cells (it looks more like the selected level was 5). In addition, the authors should explain how they obtained the graphs on the right. It looks like they underwent an extra transformation? Or are they just distorted?

Line 459, and together with the comment above: It would be helpful to be able to check the cluster numbers. For example, are cilia-associated protein cells from Cl27? A dotplot or heatmap might be more convincing to support these statements.

In addition, please add a reference to Figure 4D here.

It was pointed out to us in other comments out that the logic in validating the sensory cells this way was logically flawed. We have replaced the graph on figure 4D by plots of the expression of the sensory genes within the neuronal subset previously used in figure 2 (figure 2 —figure supplement 8, materials and methods – lines 944-948). This way, we show that the sensory cells are all included in the neuronal clusters although not all of them are cluster specific (for example fmrf, figure 2, figure 2 —figure supplement 8).

Line 468-470: egf-like does not seem to be expressed in Cl03 and 31 from Figure 5A as suggested in the text.

We accidentally mixed up some of the clusters and genes in this section. We thank the reviewer for spotting these mistakes. egf-like if strongly expressed in cl08 and cl09 whereas Cl03 and Cl04 are characterized by the expression of papilin. Cl31 expresses the epithelial splicing regulatory protein esrp.

Lines 471-493: please add the discussed genes to the dotplot in Figure 5A or supplement.

We now added the additional genes discussed in the text to the dotplot in figure 5A

Line 541: It seems that cells in Cl27 do express neuronal markers. Please revise.

We revised the sentence and stated the expression of neuronal markers in Cl27 as described in other sections. The sentence reads as follows: " While madnf-expressing cells of Cl27 co-express neuronal markers, the cells of Cl23 and Cl31 do not" (lines 619-620).

Line 545 and onwards: please not only discuss the results obtained but place the discussion in the state-of-the-art literature. Several statements have already been proposed in other papers which deserve to be referenced (e.g. lines 559, 568). Multiple genes mapped with HCR in this manuscript have been described in other cephalopods and assigned a preliminary function (e.g. Buresi 2013, 2014; Shigeno 2015; Koenig 2016; Focareta 2016; Deryckere 2021 and others). Please acknowledge their presence in appropriate sections, and, if the authors which, they might add a paragraph describing consistency or differences in expression.

As mentioned in response to public reviews, we have extensively developed the discussion to include more references, among which we can place our results.

Line 562: in line with previous comments in the results, please revise. This is based on the functionality of asc in Drosophila, but not other animals. Additionally, in the opinion of the authors, what are the asc+ e2f3- cells?

We now discuss these results, taking into account the differences in the function of asc across different animal phyla (lines 709-721). We propose that these cells are early self-replicating neuronal progenitors.

Line 713: please add the parameters used to filter low-quality cells.

The parameters used are now included in this section (lines 928-932). In addition, plots relevant to the assessment of the quality of the scRNAseq data are now provided in figure 1—figure supplement 1.

Line 722: please add additional information on how cells were subclustered (and re-clustered?). Please add a description of the slingshot analysis.

We now provide more detailed information about the clustering parameters for the full dataset as well as the subclustering used in figure 2,3 and 4 (lines 944-957).

Figures:

In general, please rearrange so each panel comes chronologically in the Results section. Make sure to mention each panel in the text.

Add axes to the confocal images and consider annotation.

We have now provided annotation and axes to facilitate the reader's understanding of the data. The panels have been rearranged to match the text.

Figure 1:

A) Is this the ventral view instead of the dorsal view? The funnel is visible?

This is indeed a ventral view, we now fixed this mistake.

Please arrange the abbreviations alphabetically. Ibl and sbl are missing.

Figure 2:

Line 790: UMAP with clusters.

C) Add legend to the 1→6 expression? Bar. It is not clear how these featureplots were calculated. Which clusters were selected? Were they re-calculated? Please show the dimplot with cluster annotation. What is on the x- and y-axes? UMAP? t-SNE?

D) Please rephrase 'that were not represented in the scRNAseq data': do the authors mean not differentially expressed? Were the genes absent from the reference transcriptome?

Figure 3:

A) Why are Cl2 and Cl14 encircled?

C) Explain what the red encircled areas are.

We have added the missing information in the figure legends.

Figure 4:

D) Double check that 3 is the threshold. The curly bracket is confusing (it reads like the right panels feed into the left). What are the x- and y-axes? Similar comment to Figure 2C.

As mentioned above figure 4D was modified to show the expression of the genes in the subset of neurons previously used in figure 2.

https://doi.org/10.7554/eLife.80670.sa2

Article and author information

Author details

  1. Jules Duruz

    Department of Biology, Institute of Zoology, University of Fribourg, Fribourg, Switzerland
    Contribution
    Formal analysis, Validation, Investigation, Visualization, Methodology, Writing - original draft
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0003-1860-9546
  2. Marta Sprecher

    Department of Biology, Institute of Zoology, University of Fribourg, Fribourg, Switzerland
    Contribution
    Formal analysis, Methodology
    Competing interests
    No competing interests declared
  3. Jenifer C Kaldun

    Department of Biology, Institute of Zoology, University of Fribourg, Fribourg, Switzerland
    Contribution
    Investigation, Methodology
    Competing interests
    No competing interests declared
  4. Al-Sayed Al-Soudy

    Department of Biology, Institute of Zoology, University of Fribourg, Fribourg, Switzerland
    Contribution
    Investigation, Methodology
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0002-7830-9660
  5. Heidi EL Lischer

    Interfaculty Bioinformatics Unit and Swiss Institute of Bioinformatics, University of Bern, Bern, Switzerland
    Contribution
    Investigation, Methodology
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0002-9616-2092
  6. Geert van Geest

    Interfaculty Bioinformatics Unit and Swiss Institute of Bioinformatics, University of Bern, Bern, Switzerland
    Contribution
    Investigation, Methodology
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0002-1561-078X
  7. Pamela Nicholson

    Institute of Genetics, University of Bern, Bern, Switzerland
    Contribution
    Investigation, Methodology, Project administration
    Competing interests
    No competing interests declared
  8. Rémy Bruggmann

    Interfaculty Bioinformatics Unit and Swiss Institute of Bioinformatics, University of Bern, Bern, Switzerland
    Contribution
    Investigation, Methodology, Project administration
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0003-4733-7922
  9. Simon G Sprecher

    Department of Biology, Institute of Zoology, University of Fribourg, Fribourg, Switzerland
    Contribution
    Conceptualization, Supervision, Funding acquisition, Writing - original draft, Project administration
    For correspondence
    simon.sprecher@unifr.ch
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0001-9060-3750

Funding

Schweizerischer Nationalfonds zur Förderung der Wissenschaftlichen Forschung (310030_188471)

  • Simon G Sprecher

Schweizerischer Nationalfonds zur Förderung der Wissenschaftlichen Forschung (IZCOZ0_182957)

  • Simon G Sprecher

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

Acknowledgements

We wish to acknowledge the staff of the station Biologique de Roscoff, France, for organizing the collection and shipment of L. vulgaris embryos. We would like to thank the staff of the Next Generation Sequencing facility of the university of Bern for helpful inputs and their help and with the single-cell sequencing experiments.

Senior Editor

  1. Didier YR Stainier, Max Planck Institute for Heart and Lung Research, Germany

Reviewing Editor

  1. Sonia Sen, Tata Institute for Genetics and Society, India

Reviewers

  1. Sonia Sen, Tata Institute for Genetics and Society, India
  2. Samuel Reiter, Okinawa Institute of Science and Technology, Japan
  3. Astrid Deryckere, Columbia University, United States

Publication history

  1. Preprint posted: March 29, 2022 (view preprint)
  2. Received: May 30, 2022
  3. Accepted: December 28, 2022
  4. Accepted Manuscript published: January 3, 2023 (version 1)
  5. Version of Record published: January 13, 2023 (version 2)

Copyright

© 2023, Duruz et al.

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

Metrics

  • 810
    Page views
  • 167
    Downloads
  • 1
    Citations

Article citation count generated by polling the highest count across the following sources: PubMed Central, Crossref, Scopus.

Download links

A two-part list of links to download the article, or parts of the article, in various formats.

Downloads (link to download the article as PDF)

Open citations (links to open the citations from this article in various online reference manager services)

Cite this article (links to download the citations from this article in formats compatible with various reference manager tools)

  1. Jules Duruz
  2. Marta Sprecher
  3. Jenifer C Kaldun
  4. Al-Sayed Al-Soudy
  5. Heidi EL Lischer
  6. Geert van Geest
  7. Pamela Nicholson
  8. Rémy Bruggmann
  9. Simon G Sprecher
(2023)
Molecular characterization of cell types in the squid Loligo vulgaris
eLife 12:e80670.
https://doi.org/10.7554/eLife.80670

Further reading

    1. Developmental Biology
    Marianne E Emmert, Parul Aggarwal ... Roger Cornwall
    Research Article Updated

    Neonatal brachial plexus injury (NBPI) causes disabling and incurable muscle contractures that result from impaired longitudinal growth of denervated muscles. This deficit in muscle growth is driven by increased proteasome-mediated protein degradation, suggesting a dysregulation of muscle proteostasis. The myostatin (MSTN) pathway, a prominent muscle-specific regulator of proteostasis, is a putative signaling mechanism by which neonatal denervation could impair longitudinal muscle growth, and thus a potential target to prevent NBPI-induced contractures. Through a mouse model of NBPI, our present study revealed that pharmacologic inhibition of MSTN signaling induces hypertrophy, restores longitudinal growth, and prevents contractures in denervated muscles of female but not male mice, despite inducing hypertrophy of normally innervated muscles in both sexes. Additionally, the MSTN-dependent impairment of longitudinal muscle growth after NBPI in female mice is associated with perturbation of 20S proteasome activity, but not through alterations in canonical MSTN signaling pathways. These findings reveal a sex dimorphism in the regulation of neonatal longitudinal muscle growth and contractures, thereby providing insights into contracture pathophysiology, identifying a potential muscle-specific therapeutic target for contracture prevention, and underscoring the importance of sex as a biological variable in the pathophysiology of neuromuscular disorders.

    1. Developmental Biology
    2. Genetics and Genomics
    Ankit Sabharwal, Mark D Wishman ... Stephen C Ekker
    Research Advance Updated

    The clinical and largely unpredictable heterogeneity of phenotypes in patients with mitochondrial disorders demonstrates the ongoing challenges in the understanding of this semi-autonomous organelle in biology and disease. Previously, we used the gene-breaking transposon to create 1200 transgenic zebrafish strains tagging protein-coding genes (Ichino et al., 2020), including the lrpprc locus. Here, we present and characterize a new genetic revertible animal model that recapitulates components of Leigh Syndrome French Canadian Type (LSFC), a mitochondrial disorder that includes diagnostic liver dysfunction. LSFC is caused by allelic variations in the LRPPRC gene, involved in mitochondrial mRNA polyadenylation and translation. lrpprc zebrafish homozygous mutants displayed biochemical and mitochondrial phenotypes similar to clinical manifestations observed in patients, including dysfunction in lipid homeostasis. We were able to rescue these phenotypes in the disease model using a liver-specific genetic model therapy, functionally demonstrating a previously under-recognized critical role for the liver in the pathophysiology of this disease.