L Douillard

V Yu Aristov

F Semond

P Soukiassian
email: psoukiassian@cea.fr

Surface Science

Keywords: Scanning tunneling microscopy, Silicon carbide, Superlattices, Surface relaxation and reconstruction, Surface structure, morphology, roughness and topography

We investigate the transition between b-SiC(100) surface reconstructions using scanning tunneling microscopy (filled and empty electronic states). We find the unexpected formation of Si-dimer lines that are self-assembled by pairs with a long-range order resulting in an even 8×2 surface reconstruction. This very uncommon surface ordering results in the combination of alternating periodic 3×2 and 5×2 unit cells. In contrast to the b-SiC(100)3×2 surface reconstruction, the Si-dimers are symmetric. Such atomic line pairs provide a novel type of one dimensional sub-nanostructure potentially having interesting characteristics.

Recent progress in silicon carbide (SiC ) surface (a), the most studied but presenting numerous different polytypes, and a cubic one ( b), more and interface experimental investigations has strengthened the already important fundamental promising but difficult to grow [2]. To date, few detailed investigations have been performed on the as well as technological interests devoted to this promising refractory material. In fact, SiC is b-SiC(100) surfaces due in part to important difficulties encountered in high-quality surface expected to bring rapidly an appropriate solution in high-power, high-frequency and high-temper-preparation [START_REF] Kaplan | Properties deposition. This could lead to the discovery of of Silicon Carbide, EMIS Datareview Series, INSPEC, nanostructures having novel and interesting[END_REF][START_REF] Bermudez | [END_REF][START_REF] Soukiassian | Physique IV (Paris) 7, C6 In conclusion, using scanning tunneling micro[END_REF]. Most of the knowledge is based on well-established techniques such as low electron ature electronic devices and sensors [1,2]. Representing so far the only stable and long-range energy diffraction (LEED), Auger electron spectroscopy (AES) and X-ray photoelectron ordered form of the IV-IV compound family, SiC has two crystallographic phases, a hexagonal one spectroscopy ( XPS) [START_REF] Kaplan | Properties deposition. This could lead to the discovery of of Silicon Carbide, EMIS Datareview Series, INSPEC, nanostructures having novel and interesting[END_REF][START_REF] Bermudez | [END_REF][START_REF] Soukiassian | Physique IV (Paris) 7, C6 In conclusion, using scanning tunneling micro[END_REF]. It is only recently that advanced probes such as scanning tunneling microscopy (STM ), medium energy ion scattering structions [START_REF] Kaplan | Properties deposition. This could lead to the discovery of of Silicon Carbide, EMIS Datareview Series, INSPEC, nanostructures having novel and interesting[END_REF][START_REF] Bermudez | [END_REF][START_REF] Soukiassian | Physique IV (Paris) 7, C6 In conclusion, using scanning tunneling micro[END_REF][START_REF] Hara | Ohdoreconstruction based on the formation of unexmari[END_REF][7][START_REF] Long | the 3×2 and 5×2 recon[END_REF][START_REF] Bermudez | of symmetric 475. Si-dimers. Such atomic lines that self-organize by[END_REF][10][11][12]. These include: (1) Si-rich 3×2 [START_REF] Kaplan | Properties deposition. This could lead to the discovery of of Silicon Carbide, EMIS Datareview Series, INSPEC, nanostructures having novel and interesting[END_REF][START_REF] Bermudez | [END_REF][START_REF] Soukiassian | Physique IV (Paris) 7, C6 In conclusion, using scanning tunneling micro[END_REF][START_REF] Hara | Ohdoreconstruction based on the formation of unexmari[END_REF][7][START_REF] Long | the 3×2 and 5×2 recon[END_REF][START_REF] Bermudez | of symmetric 475. Si-dimers. Such atomic lines that self-organize by[END_REF][10][11][12]. The surface quality was checked by LEED to exhibit sharp diffraction patterns. Other experi-and derived (2n+1)×2 series, (2) stoichiometric

S U R F A C E S C I E N C E L E T T E R S
Si-terminated c(4×2)/2×1, (3) stoichiometric mental details can be found elsewhere [11,12]. We first show in Fig. 1 the effect of thermal C-terminated c(2×2) and (4) C-rich 1×1 [START_REF] Kaplan | Properties deposition. This could lead to the discovery of of Silicon Carbide, EMIS Datareview Series, INSPEC, nanostructures having novel and interesting[END_REF][START_REF] Bermudez | [END_REF][START_REF] Soukiassian | Physique IV (Paris) 7, C6 In conclusion, using scanning tunneling micro[END_REF][START_REF] Hara | Ohdoreconstruction based on the formation of unexmari[END_REF][7][START_REF] Long | the 3×2 and 5×2 recon[END_REF][START_REF] Bermudez | of symmetric 475. Si-dimers. Such atomic lines that self-organize by[END_REF][10][11][12]. The (2n+1)×2 (n≥1) reconstructions, i.e. 3×2, annealing sequences on the b-SiC(100)3×2 surface structure using STM (filled electronic states). 5×2,... were found to result from the long range ordering of Si dimer rows lying on a full Si plane Fig. 1a shows a 150 A ˚×150 A ˚topograph of the Si-rich 3×2 reconstruction exhibiting dimer rows [10]. These Si dimer rows, representing a total of 1/(2n+1) Si monolayer (ML), consist of Si dimers with the characteristic oval shape spots associated with silicon dimers [11]. Fig. 1b andc show repre-stacked perpendicularly to the row axis. Arranged with a lateral spacing of n times the surface lattice sentative STM topographs of the 3×2 surface modifications (Fig. 1a) after short and longer parameter (a=3.08 A ˚), these Si lines define the whole (2n+1)×2 (n≥1) reconstruction sequence, annealing times at 1150°C. On the latter STM images, one can clearly observe that annealing ultimately ending in the stoichiometric Siterminated c(4×2) reconstruction (n +2).

results in Si atom removal leading to the formation of straight parallel lines arranged in a more or less Since only odd Si-rich b-SiC(100) reconstructions have been observed so far, a 3×2 to c(4×2) dense pattern. In Fig. 1b, two different types of line ordering are observed. The first one includes transition model having no possible even reconstruction was proposed [10].

a dense array (marked A in the margin) with a lateral distance between lines of about 9.2 A ˚, i.e. In this Letter, we provide evidence of a novel Si dimer row ordering pattern. This arrangement 3a (a being the surface lattice parameter) and corresponds to a remaining 3×2 reconstruction consists of long-range ordered Si-dimer lines that are self-assembled by pairs upon thermal annea-area. The second one is made of a less dense array where lines are very interestingly set by pairs that ling. This unexpected pair arrangement which, as far as we know, has never been observed before are regularly spaced on the surface. The effect of further annealings can be traced in Fig. 1c, which on any surface, occurs between the 3×2 and 5×2 surface orderings and results in an even 8×2 exhibits two different types of line patterns with, as in Fig. 1a andb, some dimer vacancies. Indeed, reconstruction of the b-SiC(100) surface. Contrary to the case of b-SiC(100)3×2, which only exhibits we observe the previously depicted line pair arrangement (marked B) together with another asymmetric dimers all tilted in the same direction [11], the Si dimers are symmetric.

and less dense single line arrangement. In the latter case, the measured lateral distance between lines The experiments were performed using an Omicron scanning tunneling microscope (STM ) at close to 15.4 A ˚(5a) indicates clearly a 5×2 reconstruction area [10]. The two latter STM topo-working pressures of 4×10-11 Torr. The STM topographs (filled and empty electronic states) graphs suggest the existence of a new b-SiC(100) surface ordering coexisting with the 3×2 ( Fig. 1b) were recorded in the constant current mode, and were reproducible and not tip-dependent. Single-and 5×2 (Fig. 1c) reconstructions. As a matter of fact, upon annealing, part of the top silicon atom domain b-SiC thin films were prepared at LETI (CEA-Technologies Avance ´es) by chemical vapor desorption leads to the formation of less dense Si dimer row reconstructions having odd (2n+1)×2,

deposition (C 3 H 8 , SiH 4
) growth on vicinal (4°) Si(100) surfaces. Surface preparation was achieved n≥1, LEED patterns that ultimately results in the formation of single self-organized Si atomic lines by annealing at 1150°C to remove native oxide as well as possible contaminants. Silicon surface stoi-as recently observed [13].

In order to explore the long-range order of this chiometry was restored by room-temperature Si deposition followed by brief annealing at 1000°C. interesting and novel surface structure, we performed LEED measurements. As can be seen in Increasing annealing times allowed the complete sequence of b-SiC(100) reconstructions to be Fig. 2a, an even 8×2 LEED pattern is precisely observed after annealing the 3×2 surface and obtained, ranging from Si-rich 3×2 to C-rich 1×1 before the 5×2 reconstruction appearance. This demonstrates the long-range ordering by pairs of these Si dimer lines in a novel 8×2 surface reconstruction. The 150 A ˚×150 A ˚STM topograph (filled electronic states) displayed in Fig. 2b shows tively. In Fig. 3a (filled states), Si dimer row pairs are regularly set in a periodic arrangement with the characteristic oval shape spot associated with independent Si dimers. Within a line pair, the line-

S U R F A C E S C I E N C E L E T T E R S
to-line distance labeled as intra-pair distance d 1 (Fig. 3a) is 9±1 A ˚. The inter-pair distance d 2 (Fig. 3a) measuring the distance between the centers of two adjacent pairs is 25±1 A ˚. The dimerto-dimer distance within a line is 5.8±0.9 A ˚, as for the 3×2 surface reconstruction [11]. The intrapair distance d 1 is close to three times the surface lattice parameter a (3.08 A ˚), whereas the inter-pair distance d 2 is close to 8a, and the dimer to dimer distance approaches 2a. The latter distance is equivalent to that observed for the odd (2n+1)×2 reconstructions and therefore supports the Si dimer stacking rows as being the common features of all Si-rich b-SiC(100) (2n+1)×2 reconstruction types. This means that Si dimer lines belonging to the same pair have the same position as for the 3×2 reconstruction. Considering two neighboring line pairs, the closest lateral line-to-line distance scales at 8-3=5, i.e. five times the surface unit cell parameter a as for the 5×2 reconstruction. Thus, the periodicity of this line pair reconstruction is 3a+5a=8a and 2a perpendicular and parallel to the dimer rows respectively, in excellent agreement with the even 8×2 LEED pattern (Fig. 2a).

We now turn to the empty electronic states STM topograph (Fig. 3b), which exhibits a slightly more complicated picture. One can observe large straight bands subdivided into three parts: a large central line surrounded by two narrower ones. Two adjacent central lines are separated by the same distance as that between two adjacent row pairs in the filled state image, i.e. one inter-pair distance d 2 =25.3±0.6 A ˚. The intra-pair distance d 1 corresponds to the distance between the two adjacent Si atoms forming a dimer [11]. For the 3×2 5×2 surface reconstructions and results in an unexpected even 8×2 reconstruction. Fig. 4 reconstruction, the Si dimers were found to be asymmetric and all tilted in the same direction, displays a schematic view of this novel b-SiC(100)8×2 surface reconstruction together with which seems to facilitate identification of the two Si atoms' contributions [11]. The height profile the 3×2 and 5×2 reconstructions for the sake of comparison. The 8×2 reconstruction consists of displayed in Fig. 3c perpendicular to the Si dimer rows exhibits a regular shape, which indicates that dimer line pairs made from stacked dimers having a direction perpendicular to the line axis. Within the two Si atoms are at the same height, in contrast to the b-SiC(100)3×2 reconstruction where one line pair, the dimer rows are separated by 3a, as for the 3×2 reconstruction. The closest distance dimers are asymmetric [11]. This dimer symmetric character for the line pairs is likely to result, in between two lines belonging to two different pairs (5a) corresponds to a 5×2 periodicity. From a the middle part of the dimer row pair, in a significant overlap between dangling bond orbitals structural point of view, the line pair reconstruction can be seen as a mixture of 3×2 and 5×2 that would in turn give an unresolved central line, as observed in Fig. 3c. In addition, this empty-unit cells. Such a structural array gives a Si coverage, H, of 1/4 ML, intermediate between 1/3 ML state topograph stresses that the dimer direction is perpendicular to the row axis as shown for the for the 3×2 [11] and 1/5 ML associated with the 5×2. Therefore, this new line pair reconstruction, 3×2 reconstruction [11]. This picture, which, once again [11], is based on direct identification of Si which is derived from the 3×2 surface reconstruction by selective Si atom organization upon ther-atom position, rules out previous suggestions [10] of dimer direction parrallel to the dimer row. mal annealing, is not only able to coexist with the Si-rich 3×2 and 5×2 surface reconstructions, but Our present STM results indicate that, in addition to the known single dimer row reconstruction also appears as a novel intermediate structure between the two latters. It demonstrates that the mode, the b-SiC(100) surface undergoes a novel reconstruction based on an interesting Si dimer only constraint about dimer line lateral position is that the closest distance between two neighboring line pair array. This unique Si-dimer pair arrangement takes place between the b-SiC(100) 3×2 and dimer rows has to be an odd multiple of the a surface parameter, which, contrary to previous Technologies Avance ´es, Grenoble) for providing b-SiC(100) samples. understanding [10], allows even (2n)×2 recon-

S U R F A C E S C I E N C E L E T T E R S
structions since the sum of two odd numbers always gives an even number.

Finally, one can anticipate that, due to their unique structure, these Si atomic line pairs would become an ideal system for metal or molecule References (inorganic and organic) adsorption studies. In this context, such a low-dimensional structure surface

  Fig. 2. (a) 8×2 LEED photograph (E p =89 eV ) of the b-SiC(100) atomic line pair reconstruction. (b) 150 A ˚×150 A flled states STM topograph of the atomic line pair b-SiC(100) reconstruction. Tip bias is V t =+3.55 V at a current I t =0.33 nA.

Fig. 1 .

 1 Fig. 1. (a) b-SiC(100)3×2 surface: 150 A ˚×150 A ˚STM topo-

Fig. 3 .Fig. 4 .

 34 Fig. 3. (a) 100 A ˚×100 A ˚filled electronic state STM topograph

Acknowledgements

Dujardin, L. Pizzagalli, C. Joachim, Phys. Rev. Lett. 78 (1997) 907. [13] P. Soukiassian, F. Semond, A. Mayne, G. Dujardin, Phys. The authors are grateful to L. di Cioccio, C. Rev. Lett. 79 (1997) 2498. Jaussaud and C. Pudda at LETI (CEA-


Keywords: Molluscan genome, nuclear receptors, GPCRs, gene expansion, simultaneous hermaphroditism

The great pond snail Lymnaea stagnalis has served as a model organism for over a century in diverse disciplines such as neurophysiology, evolution, ecotoxicology and developmental biology. To support both established uses and newly emerging research interests we have performed whole genome sequencing (~ 176 x depth), assembly and annotation of a single individual derived from an inbred line.

These efforts resulted in a nal assembly of 943 Mb (L50 = 257; N50 = 957,215) with a total of 22,499 predicted gene models. As a rst step towards understanding the hermaphroditic reproductive biology of L. stagnalis, we identi ed molecular receptors, speci cally nuclear receptors (including newly discovered 2xDNA binding domain-NRs), G protein-coupled receptors, and receptor tyrosine kinases, that may be involved in the cellular speci cation and maintenance of simultaneously active male and female reproductive systems. A phylogenetic analysis of one particular family of GPCRs (Rhodopsin neuropeptide FMRFamide-receptor-like genes) shows a remarkable expansion that coincides with the occurrence of simultaneous hermaphroditism in the Euthyneura gastropods. As some GPCRs and NRs also showed qualitative differences in expression in female (albumen gland) and male (prostate gland) organs, it is possible that separate regulation of male and female reproductive processes may in part have been enabled by an increased abundance of receptors in the transition from a separate-sexed state to a hermaphroditic condition. These ndings will support efforts to pair receptors with their activating ligands, and more generally stimulate deeper insight into the mechanisms that underlie the modes of action of compounds involved in neuroendocrine regulation of reproduction, induced toxicity, and development in L. stagnalis, and molluscs in general.

Introduction

The ability of animals to perceive and react to their environments is key to their survival and evolutionary tness. The study of such sensory abilities across the metazoan tree of life has led to deep insights across many elds and has deepened our understanding of how the abilities to taste, hear, touch, see and smell evolved [1][2][3][4] . Similarly, highly sensitive and consistent cellular responses to a local environment are prerequisites for the operation of robust developmental and reproductive programs. With the growing library of metazoan genome sequences, it is clear that a rich diversity of molecular receptors underlies many aspects of these sensory abilities [5][6][7] . For example, the family of receptors known as G proteincoupled receptors (GPCRs) is one of the largest protein superfamilies in the human genome 8 . However, our general understanding of the ligands that bind these receptors in most cases is lacking, and concerted efforts to pair these receptors with their activating ligands (increasingly with new technologies such as arti cial intelligence), are underway [9][10][11][12] . Given a history across a variety of disciplines (from behavioural, to neural and molecular -reviewed in 13 ), the great pond snail L. stagnalis has the potential to make signi cant contributions to these efforts. Many of these previous investigations were focused on basic biological processes such as respiration, heart rate, egg laying, mating, feeding and more, and have also resulted in Lymnaea becoming a standard species for ecotoxicology studies (OECD TG243 [START_REF] Ducrot | Development and validation of an OECD reproductive toxicity test guideline with the pond snail Lymnaea stagnalis (Mollusca, Gastropoda)[END_REF] ).

More recently, knowledge of fundamental neuro-biological processes has provided important insight into learning and memory in L. stagnalis, and how these processes are disturbed by human and/or environmental factors (e.g. [START_REF] Hussein | Effect of photoperiod and light intensity on learning ability and memory formation of the pond snail Lymnaea stagnalis[END_REF]16,[START_REF] Rivi | What can we teach Lymnaea and what can Lymnaea teach us[END_REF] ). Likewise, ecotoxicological compounds and immune elicitors in the environment have been shown to affect this species in terms of behavioural, developmental, endocrine, immunological and oxidative stress responses (e.g., 18,[START_REF] Brix | Investigations into the mechanism of lead toxicity to the freshwater pulmonate snail, Lymnaea stagnalis[END_REF][START_REF] Bouétard | Impact of the redox-cycling herbicide diquat on transcript expression and antioxidant enzymatic activities of the freshwater snail Lymnaea stagnalis[END_REF][START_REF] Côte | Genetic variation of Lymnaea stagnalis tolerance to copper: A test of selection hypotheses and its relevance for ecological risk assessment[END_REF][START_REF] Karakaş | Accumulation and histopathological effects of cadmium on the great pond snail Lymnaea stagnalis Linnaeus, 1758 (Gastropoda: Pulmonata)[END_REF][START_REF] Seppälä | Examining adaptive evolution of immune activity: opportunities provided by gastropods in the age of 'omics[END_REF][START_REF] Seppälä | Transcriptome pro ling of Lymnaea stagnalis (Gastropoda) for ecoimmunological research[END_REF] ). Importantly, some of these responses have been linked to molecular pro les in the central nervous system (e.g. [START_REF] Feng | Transcriptome analysis of the central nervous system of the mollusc Lymnaea stagnalis[END_REF][26][START_REF] Van Nierop | Identi cation and Functional Expression of a Family of Nicotinic Acetylcholine Receptor Subunits in the Central Nervous System of the Mollusc Lymnaea stagnalis[END_REF][28] ). For example, recent work has highlighted the role of the neurotransmitter serotonin during development in the egg (e.g. [START_REF] Ivashkin | Serotonin Mediates Maternal Effects and Directs Developmental and Behavioral Changes in the Progeny of Snails[END_REF][START_REF] Ivashkin | Transglutaminase Activity Determines Nuclear Localization of Serotonin Immunoreactivity in the Early Embryos of Invertebrates and Vertebrates[END_REF] ), for which the extensively documented ontology development of this species has proven valuable (e.g. 31, 32- 35 ). Lymnaea has also proven to be a valuable model for investigations focused on the cellular and molecular mechanisms of memory [36][START_REF] Martens | Stressful stimuli modulate memory formation in Lymnaea stagnalis[END_REF][38] , biomineralisation [START_REF] Cerveau | A survey of miRNAs involved in biomineralization and shell repair in the freshwater gastropod Lymnaea stagnalis[END_REF][START_REF] Herlitze | Molecular modularity and asymmetry of the molluscan mantle revealed by a gene expression atlas[END_REF][START_REF] Hohagen | An ancient process in a modern mollusc: early development of the shell in Lymnaea stagnalis[END_REF] , and has been instrumental in elucidating the key genes involved in body-plan chirality, thanks to a naturally occurring sinistral mutant line that could be back-crossed, and experimentally and transcriptionally investigated alongside the normal, dextral type [START_REF] Abe | The development of CRISPR for a mollusc establishes the formin Lsdia1 as the long-sought gene for snail dextral/sinistral coiling[END_REF][START_REF] Kuroda | Diaphanous gene mutation affects spiral cleavage and chirality in snails[END_REF][START_REF] Shibazaki | Body handedness is directed by genetically determined cytoskeletal dynamics in the early embryo[END_REF][START_REF] Davison | Formin Is Associated with Left-Right Asymmetry in the Pond Snail and the Frog[END_REF] . A similar approach can be used to identify the gene that is responsible for the recently discovered shell colour variant, referred to as Ginger 46 .

A fascinating feature of L. stagnalis that has received increasing attention over the past decades is its reproduction. As a simultaneous hermaphrodite, these snails can mate in the male and female role, but within a single mating interaction one individual must choose which role to perform; roles can subsequently be immediately swapped with the same or a different partner [START_REF] Koene | Coolidge effect in pond snails: male motivation in a simultaneous hermaphrodite[END_REF] . The particulars of this behaviour allow for experimentally disentangling male and female processes and underlying motivations [48][START_REF] Moussaoui | Female behaviour prior to additional sperm receipt in the hermaphroditic pond snail Lymnaea stagnalis[END_REF][START_REF] Daupagne | Disentangling female postmating responses induced by semen transfer components in a simultaneous hermaphrodite[END_REF][START_REF] Álvarez | Cognition and Its Shaping Effect on Sexual Con ict: Integrating Biology and Psychology[END_REF] , something that is di cult to do in many other hermaphroditic species [START_REF] Charnov | Simultaneous hermaphroditism and sexual selection[END_REF][START_REF] Anthes | Bateman Gradients in Hermaphrodites: An Extended Approach to Quantify Sexual Selection[END_REF] . As a result, the great pond snail has played a central role in driving evolutionary research on hermaphroditic reproductive processes including sel ng, mate choice, sexual con ict, and sperm competition [START_REF] Nakadera | Receipt of Seminal Fluid Proteins Causes Reduction of Male Investment in a Simultaneous Hermaphrodite[END_REF][START_REF] Koene | Male accessory gland protein reduces egg laying in a simultaneous hermaphrodite[END_REF][56][START_REF] Puurtinen | Predominance of outcrossing in Lymnaea stagnalis despite low apparent tness costs of self-fertilization[END_REF][58][START_REF] Coutellec | Heterosis and inbreeding depression in bottlenecked populations: a test in the hermaphroditic freshwater snail Lymnaea stagnalis. Genetic load in a freshwater snail[END_REF] .

To advance these elds, the Lymnaea research community requires an annotated genome sequence to facilitate the use of contemporary sequencing methods (RNA-Seq, single cell sequencing, ATACseq), proteomics, comparative genomics and functional genomics (e.g., RNAi and CRISPR-Cas9). Moreover, such a resource will be the foundation for the identi cation and characterisation of novel molecular receptors. The process of ligand/receptor de-orphanisation will enable the study of their role in basic biological processes and inform the ecological risk assessment of chemicals such as suspected endocrine disrupters (see the OECD framework for endocrine disruption; 60,61 ). In addition, such an inventory will provide a framework to investigate whether this hermaphrodite (and others) deviates from separate-sexed species in its use and diversity of receptors, and how such life history strategies may have evolved.

To promote further progress, we sequenced, assembled and annotated the genome of L. stagnalis, and searched for gene family expansions using orthology analyses. We focused on identifying and analysing nuclear receptors (NRs), G-protein coupled receptors (GPCRs), and receptor tyrosine kinases (RTKs). These three highly conserved superfamilies were selected based on their pivotal roles in endocrine systems (NRs: 62 ), signalling pathways and reproduction (GPCRs: 63 ); (RTKs: 64,65 ), and possible crosstalk among the signalling pathways that are activated by these different types of receptors 66 . For a subset of these receptors, we performed in situ hybridisation (ISH) on early developmental stages and histological sections of sex-speci c organs of adult L. stagnalis to localise their expression. Our annotation effort complements a previous effort devoted to ion-channels and ionotropic receptors in the central nervous system 67 . Here we focus on the largest family of receptors, GPCRs, to assess the difference between simultaneously hermaphroditic gastropod molluscs (such as L. stagnalis) and gastropods that manifest one sex at a time. Given that GPCRs are involved in numerous behaviours across the Metazoa including copulation and egg-laying 68 , we anticipate that our approach will contribute to evolutionary hypotheses underlying the evolution and diversi cation of receptors and their roles in maintaining the simultaneous hermaphroditic state.

Results

Genome sequencing and assembly

The Lymnaea stagnalis genome was sequenced using a combination of both Illumina paired-end and mate-pair reads. The resulting draft assembly spans 943 Mb and comprises 6,640 scaffolds with lengths ranging from 2,000 to 7,106,496 bp, as outlined in Table 1. This assembly is 250 Mb lower than the previously measured genome size of 1.193 Gb (C-value = 1.22 pg 69 ). This discrepancy is consistent with the general trend observed in eukaryotes (see [START_REF] Elliott | What's in a genome? The C-value enigma and the evolution of eukaryotic genome content[END_REF] ).

Mitochondrial genome

BLAST searches returned mitogenomic sequence similarities for the full length of Lsta_scaffold2639. As indicated by a 127 bp sequence overlap at the termini, this scaffold represents a circular sequence. Nanopore reads con rmed that the assembled sequence is contiguous, and gene annotations revealed the complete standard gene complement for a metazoan mitogenome. The 13,834 bp mitogenome of L. stagnalis (Genbank MW221941) is 28.21% GC-rich, different from the 37.55% AT content of the nuclear genome. Several protein coding genes are delimited by downstream tRNA genes, causing the ORF to terminate with a truncated stop codon (T-or TA-) that is completed by polyadenylation of the mRNA transcripts. This mitogenome sequence has 94% identity with another, independently characterised L. stagnalis mitogenome (13,807 bp, MT874495). The different lengths of these two mitogenomes are mostly due to indels in regions with tRNA-and rRNA genes, a two bp indel occurs at the 5' start of ND4L.

These mitogenome sequences align closely but contain abundant single nucleotide variants (including several non-synonymous substitutions) within protein coding genes. The organisation of the 22 tRNA genes, two rRNA genes plus 13 protein coding genes (Supplementary Figure SF1) of the L. stagnalis mitogenome is generally similar to that of other lymnaeid (and hygrophilid) snail species, with a few minor differences in the location of tRNA genes.

Receptor identi cation

Receptor Tyrosine Kinases A total of 32 gene predictions (including manually curated versions) were identi ed as Receptor Tyrosine Kinases (RTKs), distributed across 18 of the 20 families or types (see Supplementary Table ST1). Among these, 12 gene predictions did not contain expected (canonical) domains. We found no gene prediction representing type III or XI RTKs, while two gene predictions could not be assigned to a speci c RTK type.

Nuclear receptors

In total, 47 nuclear receptors (NRs) were identi ed from automatic and manual gene predictions, 39 of which containing the complete sequence of both conserved domains, i.e., the DNA-binding domain (DBD) and ligand binding domain (LBD) (Supplementary Tables ST2 andST3). All 6 traditional NR families were represented as follows: NR1 (28 predictions), NR2 (12), NR3 (3), NR4 (1), NR5 (2), NR6 (1). However, in terms of NR nomenclature, these corresponded to only 29 different NRs, as several gene predictions were found to have the same annotation (e.g., NR1F1_RORA, represented by four different gene predictions, that share low levels of domain sequence similarity (from 28.9-41.5%).

Signi cantly, we found two NRs with 2DBDs and one LBD. The occurrence of such an atypical NR structure was rst reported in 2006 in the parasitic platyhelminth Schistomosa mansoni [START_REF] Wu | Schistosoma mansoni (Platyhelminthes, Trematoda) nuclear receptors: sixteen new members and a novel subfamily[END_REF] and further con rmed in other animal species, including protostomes and deuterostomes (both non-chordate and chordates; see [START_REF] Wu | Updated knowledge and a proposed nomenclature for nuclear receptors with two DNA binding domains (2DBD-NRs)[END_REF] ), suggesting an origin that predates the protostome-deuterostome split. Using the metazoan-scale orthology analysis performed in this study, the two sequences identi ed in L. stagnalis belong to the same orthogroup (OG00004255), which is composed of 33 sequences from molluscan (Gastropoda and Bivalvia), annelid, and platyhelminth taxa, plus the hemichordate Saccoglossus kowalevskii. However, in this OG, no protein other than the two L. stagnalis' possess 2 DBDs, which may re ect con icts in orthology search due to this very particular NR structure. Although a stable phylogeny of 2DBD-NRs is still to be established, Wu and LoVerde [START_REF] Wu | Updated knowledge and a proposed nomenclature for nuclear receptors with two DNA binding domains (2DBD-NRs)[END_REF] proposed to group them into a new family (NR7). Accordingly, we found higher sequence similarity among the four DBDs of these two 2DBD-NRs than with any other L. stagnalis NR or with their corresponding best hits (see green clade in Fig. 1 and Supplementary Figure SF2).

To take the analysis of the unusual 2DBD-NRs further, we concatenated the DBD and LBD domains from 60 other 2DBD-NRs representative of the metazoan groups studied by Wu and Loverde [START_REF] Wu | Updated knowledge and a proposed nomenclature for nuclear receptors with two DNA binding domains (2DBD-NRs)[END_REF] (see Supplementary Table ST4). A phylogenetic analysis of these domains suggests that the two L. stagnalis receptors are members of two different groups (A vs B; Fig. 2 and Supplementary Figure SF3), albeit with low statistical support (see branch bootstrap values). Therefore, following the recommendation of Wu et al. 2023, L. stagnalis has two NR7 genes, GSLYST00003556001 (NR7A) and GSLYST00001311001 (NR7B), corresponding respectively to C. gigas 2DBDγ and 2DBDδ [START_REF] Vogeler | The nuclear receptor gene family in the Paci oyster, Crassostrea gigas, contains a novel subfamily group[END_REF] .

Among the incomplete NRs, four genes were found to lack either the entire LBD (1 gene) or the DBD (3 genes) and should thus be classi ed as members of the NR0 family (although this grouping is not monophyletic). We tentatively annotated them following their BLAST best-hit (H. sapiens and/or nr). On this basis, the single DBD-only gene was annotated as NR1A2_THRB, despite the lack of a true thyroid system in Mollusca but see [START_REF] Morthorst | Thyroid-like hormone signaling in invertebrates and its potential role in initial screening of thyroid hormone system disrupting chemicals[END_REF] , and the three LBD-only genes as NR1D2_REVERBB, NR1H3_LXRA, and NR2F1_COUP-TFI, respectively. Note that NR1D2 and NR1H3 annotations were shared by at least one other gene prediction which contained both the DBD and LBD. Another speci c case was GSLYST00019008001 (and its related transcript). Although we could not nd an open reading frame (ORF) that would generate both the DBD and LBD, two alternative ORFs would generate sequences with either of the domains. Ambiguously, each of these sequences has a different BLAST best-hit in H. sapiens: NR1C1_PPARA for the DBD-only ORF, and NR1H3_LXRA for the LBD-only ORF. In addition to the observation that NR1C1 and NR1H3 are already represented by two and three other gene predictions respectively, this suggests that GSLYST00019008001 and its corresponding transcript do not constitute a functional NR.

Despite extensive searches for alternative ORFs, 5 NRs were found to have an incomplete domain sequence (lacking part of either the DBD or the LBD). These might represent non-functional gene copies, as suggested by the observed density of nucleotide repeats in the problematic areas. While this should be con rmed through speci c re-sequencing, the fact that some of these and other genes with complete domains share the same functional annotation (see NR1C3_PPARG and NR2E1_TLL, Supplementary Table ST2) tends to support the hypothesis of non-functionality. One last speci c situation is that of NR2E3 (PNR), which is represented by two protein predictions, neither of which have both domains in full (an incomplete DBD, and the LBD is either complete or missing).

G protein-coupled receptors

InterProScan identi ed a signi cant GPCR hit for a total of 892 gene predictions, four of which were discarded due to incompleteness or other con icting annotations (see Supplementary Table ST5). The remaining 888 gene predictions are distributed among the ve main GPCR families (GRAFS system [START_REF] Schiöth | The GRAFS classi cation system of G-protein coupled receptors in comparative perspective[END_REF] )

and the ancestral cAMP 76 : class A -Rhodopsin-like (784); class B -adhesion/secretin (65); class C -Glutamate (31); class E -cAMP (3); class F -Frizzled (5). Of these, 418 were members of ve orthogroups (OGs) found to be expanded within the L. stagnalis genome (orthology analysed at metazoan level, Z-score > 2, Table 3; see also Supplementary Information for further details), while the remainder (470) was distributed across 135 non-expanded OGs. Most gene predictions belonging to the ve expanded GPCR OGs were identi ed as Rhodopsin-like (class A) 7 trans-membrane (7TM; IPR17452). These ve expanded OGs encompassed a high proportion (89 to 100%) of molluscan proteins from the six taxa included in the orthology analysis, and four of these OGs were predominantly exclusive to gastropods (85 to 100%; Table 3). 7TM-class A Peptide GPCRs: a predominance of FMRFamide R-like Three of the ve expanded OGs with a GPCR hit (OG000008, 0G0001159, OG0002276) consisted exclusively of Rhodopsin-like receptors from the group of (7-transmembrane-domain) 7tmA_FMRFamide_R-like receptors (peptide), representing in total 424 L. stagnalis proteins. However, OG0002279, consisting of 49 exclusively Gastropoda proteins, was heterogeneous with only 2 GPCR hits for L. stagnalis proteins. The 22 other L. stagnalis proteins in this OG included BLASTp annotations against bacterial proteins, casting doubt on the statistical reliability of this OG. By contrast, OG0002601, although not expanded in L. stagnalis, was clearly an FMRFamide receptor OG (45 proteins), apparently absent in Deuterostomia and Arthropoda. Molluscs were the most represented taxa in this group, with 38 proteins, 22 of which belong to Gastropoda, and 4 to L. stagnalis. OG0010705 is also a small, nonexpanded, FMRFamide receptor OG (8 proteins), with 7 molluscan representatives, all gastropods from the Heterobranchia order. However, the fact that one protein from this OG belongs to the arthropod F. candida questions the validity of this small group as a true OG. Finally, an extra-set of 44 L. stagnalis gene predictions from various and heterogeneous OGs were annotated as FMRFamide_R-like, although only 18 gene predictions had a complete 7TM motif (as determined with NCBI Conserved Domains database search tool). In summary, about 450 GPCRs could be con dently designated as FMRFamide_Rlike. For completion, an OG of biogenic amine GPCRs (G0004949), which was also expanded in L. stagnalis, contained 13 receptors to muscarinic acetylcholine and one receptor to adrenaline.

We hypothesised that FMRFamide receptor proliferation may have supported the evolution of simultaneous hermaphroditism in gastropods. To test this hypothesis, we re ned the analysis of orthology by focusing on gastropods, and used available genomic resources to balance hermaphroditic species (A. californica, B. glabrata, Bulinus truncatus, Candidula unifasciata, Elysia chlorotica, L. stagnalis, Radix auricularia, Lottia gigantea) against gonochoristic species (separate genders, Batillaria attramentaria, Gigantopelta aegis, Haliotis rubra, Littorina saxatilis, Pomacea canaliculata, Potamopyrgus antipodarum) (see Supplementary Information, Orthology analyses). Note that L. gigantea is a protandrous sequential hermaphrodite while all other species are simultaneously hermaphroditic. This analysis generated nearly 29,500 OGs (Fig. 3; Supplementary Table ST6). Because all simultaneous hermaphrodites in our analysis are in the Euthyneura infraclass of molluscs, we also identi ed various

OGs speci c to this clade, including FMRFamide receptors, c-type lectins, cytochrome P450 cyclodipeptide synthase-associated, (see Supplementary Table ST7 for domain annotation of the main euthyneura-speci c OGs).

The search for FMRFamide-R among L. stagnalis orthologs retrieved 483 genes distributed in 89 different OGs. FMRFamide-R groups previously found expanded in L. stagnalis appeared no more enriched, which is not surprising as all members of the analysis were gastropods. We excluded Radix auricularia from this analysis because that dataset was abnormally small (less than 18,000 protein predictions in total), and that species was recurrently absent from OGs otherwise well represented by other Hygrophila species (casting further doubt on the completeness of its predicted proteome). Excluding R. auricularia, the number of FMRF-R amide genes was signi cantly higher (ANOVA F-value [1,11] Peptide GPCRs: a high level of ligand diversity OG0000001 is comprised of peptide GPCRs and is the second largest OG determined by OrthoFinder with 3,444 genes representative of all taxa included. Although not expanded in L. stagnalis, this OG contains 139 L. stagnalis gene predictions, all assigned as peptide GPCRs, and predicted to bind a great diversity of molecules (more than 40 categories; see supplementary Table ST8). OG0000121 grouped gene predictions which shared four annotations with OG0000001 (in bold in Supplementary Table ST8) or were identi ed as other peptide receptors (neurotensin, somatostatin, FF neuropeptide, NPY peptide, opioid).

Among the rhodopsin-like gene predictions, a structural pattern was observed in the genome, with a spatial clustering of similar genes. Whether all members of a given cluster are functional genes or pseudogenes cannot be determined based on our data, yet it is to be noted that they do not all possess a complete 7TM motif and these clusters were also often interspersed with repetitive regions. In line with its high level of occurrence, FMRFamide_R-like was either the main or exclusive annotation of 30 of the 46 largest groups (size: 4 to 26 tandem GPCR genes), and also appeared in other groups. Among these, the proportion of transmembrane domains with at least 6 helices varied from 50 to 100%. Furthermore, 10 of the 30 groups of FMRFamide_R-like genes contained genes which also exhibited an additional viral chemokine GPCR hit (PHA02638 or PHA03087). Details on cluster annotations are given in Supplementary Table ST9. ST10). Several other predicted proteins with signi cant matches to GPCRs from either of these classes were also identi ed but discarded due to incomplete or total lack of a transmembrane domain.

Spatial expression of FMRFamide receptors

While the association between FMRFamide receptor expansion and simultaneous hermaphroditism is strong (as indicated by OG composition), it is problematic that the Euthyneura are monophyletic and are also the only simultaneous hermaphrodites used in these analyses. Therefore, to clarify whether some of these receptors are associated with the hermaphroditic condition in L. stagnalis (as opposed to other euthyneuric traits such as detorsion), we searched for differences in the expression of receptors in our transcriptomic datasets, and accordingly a selection was used for in situ expression analysis. The differential expression of some FMRFamide receptors between male (prostate) and female (albumen) reproductive organs drove our selection of 4 of these receptors (GSLYST00019397001, GSLYST00001331001, GSLYST00016909001, GSLYST00019383001). In addition, we investigated the expression patterns of three differentially expressed "traditional" nuclear receptors (GSLYST00008490001, GSLYST00009627001 and GSLYST00020139001) and the two 2DBD-NRs (GSLYST00001311001 and GSLYST00003556001). We also included ovipostatin (GSLYST00000615001) as a reference because it is known to be differentially upregulated in the prostate gland relative to the albumen gland [START_REF] Swart | Temporal expression pro le of an accessory-gland protein that is transferred via the seminal uid of the simultaneous hermaphrodite Lymnaea stagnalis[END_REF] . The spatial expression pro les for these genes were striking in their similarity in early veligers; in general, these genes were robustly expressed in ectodermal cells of the foot, head and mantle edge (Fig. 4). The albumen and prostate glands are fully differentiated organs with separate reproductive functions, and we were able to detect qualitatively distinct spatial expression patterns for these genes in one or both reproductive organs. As expected, ovipostatin generated a strong signal in sections of the prostate gland, while effectively no signal was detected in the albumen gland (Fig. 4). In addition, the expression patterns of the two 2DBD-NRs in early embryos (early and late blastulae, early and late gastrulae) indicates that these gene products play likely critical roles in early developmental processes (Fig. 5).

Discussion

Our characterisation of the genome of L. stagnalis expands the available range of organisms to study molluscan biology in the genomic era and provides a fundamental resource for research involving L. stagnalis in its role as a longstanding, invertebrate model organism. Speci cally, our analysis revealed that the nuclear genome of L. stagnalis encodes an abundance of receptors, including DNA-binding-, 2DBD-, RTK-receptors and GPCRs. We focused our attention on the identi cation of receptors involved in biological processes because such information is lacking in many elds that employ Lymnaea as a model. For example, in ecotoxicology mechanistic evidence for toxicity pathways is critically lacking in L. stagnalis, even though it is used as a model for its ecological value (functional representativity within food-webs and in other ecological interactions). The resources provided in this work will be of particular interest for the study of endocrine disruptors (nuclear receptors, domain annotation and phylogeny) and more generally for xenobiotic responsive pathways involving GPCRs (e.g., 78 ) or RTKs 65 .

We have focused our analysis of the L. stagnalis genome on the identi cation of molecular receptors, and the expression, evolution and diversi cation of a selection of these. We have also characterised the mitogenome (see Supplementary Information), two unique double DNA-binding domain (2DBD) containing genes [START_REF] Wu | Updated knowledge and a proposed nomenclature for nuclear receptors with two DNA binding domains (2DBD-NRs)[END_REF] , and investigated the in situ expression pro les of several receptors during embryogenesis and in sex-speci c adult organs of this simultaneous hermaphrodite. These ndings complement a previous study that concentrated on ion-channels and ionotropic receptors in the central nervous system (CNS, 67 ). Perhaps most notably, we observe an association between the evolution of simultaneous hermaphroditism and the diversi cation of GPCRs.

Identi cation of potential male/female receptors

As L. stagnalis is a simultaneous hermaphrodite, meaning each individual can execute male and female reproductive processes (see below), we searched for receptors that are potentially involved in regulating this simultaneous state. We used RNAseq data from albumen gland (AG; female) and prostate glands (PG; male) to look for sex-speci c expression patterns. Indeed, some GPCRs and NRs were exclusively expressed in one or the other of these organs (Supplementary Table ST11). These potential sex-speci c receptors are particularly interesting because these are potential targets of seminal uid proteins (LyACPs, Lymnaea Accessory Gland Proteins) which have been shown to speci cally in uence male or female reproductive processes in the sperm recipient (e.g. [START_REF] Koene | Male accessory gland protein reduces egg laying in a simultaneous hermaphrodite[END_REF][START_REF] Nakadera | Receipt of Seminal Fluid Proteins Causes Reduction of Male Investment in a Simultaneous Hermaphrodite[END_REF] ). In the future, it will be fruitful to test the regulatory status of these receptor genes, how transient their expression is, and whether this is in uenced by mating status (e.g., recently inseminated). We expect that such efforts will not only result in deorphanising these receptors but will also present new opportunities for understanding the proximate mechanisms via which LyACPs induce their effects. This also argues for the generation of more wellannotated chromosome level gastropod genomes, and for the functional analysis of speci c gene products.

Receptor expression

The broadly similar spatial expression patterns of the 2 NRs, 3 GPCRs and ovipostatin in early veliger larvae is intriguing (Fig. 4), and hints at potential functional commonalities for these distinct genes at this stage of development. These similarities are made more remarkable by the fact that ovipostatin is not a receptor, but rather a secreted accessory gland protein produced by the prostate gland that is understood to repress oviposition behaviours in Lymnaea [START_REF] Swart | Temporal expression pro le of an accessory-gland protein that is transferred via the seminal uid of the simultaneous hermaphrodite Lymnaea stagnalis[END_REF] . These genes expand their expression from an apparently exclusive ectodermal pattern in early veligers, to the endodermally-derived adult reproductive organs we investigated (Fig. 4). The presence of mRNA transcripts for the selected receptors in one or both of these terminally differentiated organs implies that they are functionally required in these contexts, and that they therefore contribute to the maintenance of simultaneously functioning male and female reproductive states. As mentioned above, gene-speci c (and likely reproduction-condition-speci c) assays are required to elucidate the functional roles of these genes in these contexts. However, these results provide a foundation from which future investigations can be started.

The expression of the 2DBD-NR retinoic acid receptor-"beta" (as referred to in [START_REF] Vogeler | The nuclear receptor gene family in the Paci oyster, Crassostrea gigas, contains a novel subfamily group[END_REF] ) gene (GSLYST00001311001) in distinct cells of gastrula embryos (Fig. 5) corresponds with earlier reports on the effect of lithium (which is known to interact synergistically with retinoic acid [START_REF] Holtz | Lithium chloride inactivates the 20S proteasome from WEHI-3B D+ leukemia cells[END_REF] ) on gastrulation and development in L. stagnalis 80 . The expression pattern of the 2DBD-NR subfamily 1 group D gene (GSLYST00003556001) in numerous nuclei of gastrula embryos also suggests an important role in early development, however, pending functional data (such as knock-downs or knock-outs) the speci c ontogenetic roles of these genes remain unknown.

GPCR evolution/expansion

In general, G protein-coupled receptors (GPCRs) are an important class of membrane receptors that are involved in many biological processes, including external (e.g., pheromones, ions, odours, toxins, taste, and light) and internal (e.g., neuropeptides, proteins, lipids, neurotransmitters, and nucleotides) sensory signal detection (e.g. 81 ). As a superfamily, GPCRs have evolved by tandem local duplication during metazoan evolution (prior to two rounds of whole-genome duplication in vertebrates 82 ). In this context it is signi cant that we found around 450 FMRFamide R-like (Rhodopsin-like) GPCRs in the Lymnaea stagnalis genome. The presence of this expanded set of receptors can be explained by the fact that they are the receptors for the many FMRFamide-like (neuro)peptides that are involved in regulating many biological processes, including the hermaphroditic condition (e.g. 83 ). These peptides belong to a class that contain the C-terminal Arg-Phe-NH2 (RFamide; RFamide-related peptides or RFRPs). The evolutionary origins of these (neuro)peptides dates back to a common bilaterian ancestor 84 . The rst of these peptides to be discovered was the cardio-excitatory tetrapeptide FMRFamide, identi ed in the neural ganglia of the bivalve mollusc, Macrocallista nimbosa 85 . Since then, FMRFamide has been reported to be present in many species, including the ecdysozoan protostome Caenorhabditis elegans 86 , and various lophotrochozoan molluscs, including Cephalopods 87, 88 such as Sepia o cinalis 89 and Octopus vulgaris 90, 91 , and Gastropods including Cornu aspersum, Biomphalaria glabrata and Biomphalaria alexandrina 92, 93 , to name but a few. In addition, many RFRPs have since been discovered and are involved, besides cardiovascular regulation, in osmoregulation, digestion, locomotion, feeding and reproduction (e.g. [START_REF] Zatylny-Gaudin | Characterization of a novel LFRFamide neuropeptide in the cephalopod Sepia o cinalis[END_REF] ).

The mechanism that supported the apparent expansion of these receptors in the Euthyneura remains unknown. Whole genome duplication (WGD) events have been identi ed within several molluscan clades (at the base of the Cephalopoda, and within the Stylommatophora, the Neogastropoda, Tonnoidea, Cypraeidae, Capulidae and Heterobranchia) but apparently not at the origin of Euthyneura [START_REF] Hallinan | Comparative analysis of chromosome counts infers three paleopolyploidies in the mollusca[END_REF][96][START_REF] Liu | Giant African snail genomes provide insights into molluscan whole-genome duplication and aquatic-terrestrial transition[END_REF] . With respect to GPCRs, due to the high level of gene clustering observed in L. stagnalis, it seems unlikely that a WGD event at the origin of the Euthyrneura would be responsible for the FMRFamide receptor proliferation we observe in euthyneuran species. However, this hypothesis should be tested further, using appropriate tools and chromosome-level datasets in the Mollusca and elsewhere in the animal kingdom.

Interestingly Rondón et al. 98 recently reported an expansion of GPCR genes in two lineages of caenogastopods, one of which belongs to a lineage that underwent a whole genome duplication. As in the Lymnaea genome they also observe clustering of expanded GPCRs, which argues against a genome duplication event being a direct explanation to the observed expansion in GPCRs.

Based on our analyses of molluscan RFRPs and GPCRs, we postulate that the patterns of their abundance may in part have supported the transition from a separate-sexed state to a hermaphroditic condition. Our phylogenetic analysis supports an expansion of GPCRs in parallel with the occurrence of hermaphroditism in the molluscan species with well annotated genomes that were available at the time of writing for comparison. Supporting this observation is the fact that many molluscan reproductive processes are regulated by RFRPs and other (neuro)peptides. As in vertebrates, separate sexed molluscs can use the same substances for regulating male and female reproductive processes, because each body only contains one sex (e.g., Octopus vulgaris [START_REF] Di Cristo | Presence of two neuropeptides in the fusiform ganglion and reproductive ducts of Octopus vulgaris: FMRFamide and gonadotropinreleasing hormone (GnRH)[END_REF] ; Sepia o cinalis: 89 ; Sepiella japonica: [START_REF] Li | Identi cation, Characterization, and Expression Analysis of a FMRFamide-Like Peptide Gene in the Common Chinese Cuttle sh (Sepiella japonica)[END_REF] ; Haliotis discus hannai: 100 ). However, the evolution of simultaneous hermaphroditism means that male and female processes must be regulated concurrently within one body (e.g., Lymnaea stagnalis, Cornu aspersum, Aplysia californica 87 ). This requires largely nonoverlapping neurobiological wiring, regulatory peptides and accompanying receptors to avoid the simultaneous expression of con icting reproductive processes 101 . Hence, excitatory and inhibitory crosstalk between the male and female regulatory systems needs to be accompanied by an expansion of the receptors involved, so that they can be mutually exclusively used for executing one sexual role (while potentially inhibiting the other role).

It is worth illustrating this with the multitude of roles that FMRFamide and its co-expressed peptides ful l 102-105 next to being cardio-excitatory 106 . These tetra-and heptapeptides have been shown to be important for the ne coordination of the execution of male mating behaviour and insemination 105,107 .

Importantly, FMRFamide has also been shown to inhibit the activity of the caudo-dorsal cells (CDC), which are normally responsible for the release of egg laying hormone during their synchronous bursting (CDCH 108,109,110 ). Interestingly, the single, previously identi ed receptor for an RFamide is a ligand for the decapeptide TPHWRPQGRF-NH 2 (LyCEP/luqin/tachykinin 111 ). This peptide is cardioexcitatory, but also inhibits CDC activity, and by extension egg laying. Thus, when the CDCs are activated and egg laying is triggered, LyCEP release is likely suppressed, which would agree with our observation of a reduction in heartrate during egg laying (Smelik et al. unpubl.). This illustrates how the different RFRPs are interlinked and can collectively stimulate one physiological function while simultaneously inhibiting con icting functions.

Conclusion

Here we report a full genome assembly and annotation for the great pond snail Lymnaea stagnalis. This snail has long served as a model organism to a wide variety of scienti c disciplines, and we envisage this resource will assist traditionally distinct communities to more e ciently collaborate, provide deeper insight into long-standing questions, and to stimulate novel research directions. The annotation of this genome highlights a remarkable expansion of receptors critical for a wide variety of biological processes, and hints at a causal link with the evolution of simultaneous hermaphroditism in gastropods. The identi cation of such receptors is also crucial for understanding the chemosensory lifestyle that many gastropods lead and will likely provide insight into the mechanisms that underlie phenomena such as pheromone release, egg-laying, predator detection, mate location, mate-partner recognition, avoidance behaviours, memory, parasite infection and co-habitation behaviours. We envisage that our ongoing efforts to identify novel receptors, in combination with parallel advances in structural in silico modelling of protein structure, protein receptor localisation techniques and functional studies of receptors will ultimately de-orphanise the many receptors present in the genome of L. stagnalis and will provide insight not only into their mechanisms of action, but into their evolutionary histories.

Materials and Methods

Snail line generation

The individual snail used for genome sequencing and assembly was derived from an inbred line previously used to map and identify the snail chirality locus [START_REF] Davison | Formin Is Associated with Left-Right Asymmetry in the Pond Snail and the Frog[END_REF] . This line was derived by alternating rounds of self-fertilisation and full-sib mating, over at least ten generations. The inbred snail line is freely available from the authors upon request.

Extraction of high molecular weight DNA One whole adult snail was ash-frozen and stored at -80°C. Upon thawing, the soft tissue was removed from the shell, cleaned in a physiological water bath and ground with a mortar and pestle under liquid nitrogen. The ground material was resuspended in 15 mL SEB1X buffer in the presence of 0.2% βmercaptoethanol and gently stirred for 15 min at 4°C. This solution was then passed through a strainer in order to remove as much mucus as possible. 0.05 vol of a SEB1X-Triton mixture was then added to the 14 mL of ltrate, incubated for 10 min at 4° C, and centrifuged for 20 min at 600g at 4° C. The supernatant was aspirated, and the resulting pellet of nuclei was washed 2 times with SEB1X buffer, then resuspended in 300 µL of S buffer in the presence of 30 µL of SDS. Buffer composition is available as Supplementary Information.

The protocol was immediately continued by embedding the nuclei in plugs. Low melting point agarose (2%) was prepared and added to the nuclei suspension, volume by volume. After gentle homogenisation with a vortex, avoiding bubbles, the mixture was poured into plugs at a rate of 100 µL per plug and left to set on ice for 15 min. The plugs were then transferred into a 50 mL tube in the presence of 10 mL lysis buffer per plug and incubated for 1h at 37°C. After removal of the lysis buffer, plugs were incubated in 10 mL of digestion buffer and proteinase-K at 40 µg/mL nal concentration at 50°C overnight. The following day, the proteinase-K solution was removed and inactivated by washing with a mixture of 10 mL TE 10:1 pH 8 and 23 µL PMSF (phenylmethylsulphonyl uoride) 0.1 M (40 µg/mL nal). Two washes of 30 min at room temperature were followed by 3 washes for 10 min each at room temperature with 20 mL of TE 10:1 pH 8. The plugs were then stored in 30 mL EDTA 0.5 M pH 8 at 4°C until extraction.

High-molecular-weight DNA was extracted from the plugs by β-agarase treatment. After 5 washes in 20 mL TE 1X buffer followed by 2 washes in 20 mL MES buffer, each plug was brie y dried on paper, placed in a 1.5 mL tube and incubated for 7 min at 70°C in a dry water bath before being incubated at 42°C overnight in the presence of 16 µL β-agarase. The following day, the entire volume (~ 500 µl) was transferred onto a dialysis lter using a broad tip and left for 30 min at room temperature. DNA quality was checked using Opgen QCards (Opgen, Rockville, MD, USA), and concentration determined via Qubit Fluorometric Quantitation (Thermo Fisher Scienti c, Waltham, MA, USA).

A set of cDNA libraries were constructed using various tissues and organs from individuals of the same strain used for genome sequencing (see Supplementary Table ST12). For the stomach, digestive gland, reproductive system, and foot, stranded mRNAs were generated using Illumina® TruSeq stranded mRNA Library Prep, whereas for nervous ganglia, head and mantle tissue, a kit adapted to low input RNA was used (SMARTer® Universal Low Input RNA Kit).

Transcriptome assembly

Illumina reads were quality ltered and trimmed using the procedure described in Alberti et al. 112 . In addition, for RNA-Seq data, ribosomal RNA-like reads were detected using SortMeRNA 113 and Oases 0.2.08 114 , using a k-mer size of 63 bp. Reads were mapped back to the contigs with BWA-mem and coverage was computed based on solely consistent paired-end reads. Uncovered regions were extracted and used to identify hypothetic chimeric contigs. In addition, open reading frames (ORF) and domains were searched using respectively TransDecoder (Haas, BJ.

https://github.com/TransDecoder/TransDecoder) and CDsearch 115 . Uncovered regions that do not contain an ORF or a domain were tagged as mis-joins and contigs were split. Read strand information was used to correctly orient RNA-Seq contigs.

Error-correction of cDNA nanopore reads Raw reads were corrected using Illumina RNA-Seq reads and the NaS software 116 with the following parameters: --ilmn_size 150 --mode fast --nb_proc 4 --t 2 --k 32 --untgl_seq_size 300. As the computation time is correlated with the number of Illumina reads, we decided to reduce the complexity by adding a normalisation step of the Illumina short-reads using the normalise-by-median.py script of the khmer package 117 with the following parameters: -C 20 -k 20 -N 4 -x 10e9.

DNA sequencing genome assembly

Library preparation and whole-genome sequencing Four Illumina PE were prepared starting from 100 ng to 250 ng DNA. Four independent DNA fragmentations were performed using the E210 Covaris instrument (Covaris, Inc., USA) in order to generate fragments mostly around 400 and 500 bp. Libraries were constructed using the NEBNext DNA Sample Prep Master Mix Set (New England Biolabs, MA, USA). DNA fragments were PCR-ampli ed using Platinum Pfx DNA polymerase (Invitrogen) and P5 and P7 primers. Ampli ed library fragments were size selected on 3% agarose gel around 400 bp or on 2% agarose gel around 600bp and 700 bp. Library traces were validated on Agilent 2100 Bioanalyzer (Agilent Technologies, USA) and quanti ed by qPCR using the KAPA Library Quanti cation Kit (KapaBiosystems) on a MxPro instrument (Agilent Technologies, USA).

The paired end (PE) libraries were sequenced with 250bp paired-end reads on an Illumina HiSeq 2500 instrument (Illumina, San Diego, CA, USA).

The Mate Pair (MP) libraries were prepared using the Nextera Mate Pair Sample Preparation Kit (Illumina, San Diego, CA). Brie y, genomic DNA (4 µg) was simultaneously enzymatically fragmented and tagged with a biotinylated adaptor. Tagmented fragments were size selected on a Sage Science Electrophoretic Lateral Fractionator or SageELF (Sage Science, MA, USA). This system allows production of narrowly sized MP libraries, isolating 12 different discrete size fractions from a single sample loading. We selected 8 fractions (from 3Kb to 15 Kb) to continue the MP library preparation. Size-selected fractions were circularised overnight with a ligase. Linear, non-circularized fragments were digested, and circularised DNA was fragmented to 300-1000-bp size range using Covaris E210. Biotinylated DNA was immobilised on streptavidin beads, end-repaired, 3`-adenylated, and Illumina adapters were added. DNA fragments were PCR-ampli ed using Illumina adapter-speci c primers and puri ed. Finally, libraries were quanti ed by qPCR and libraries pro les were evaluated using an Agilent 2100 bioanalyzer (Agilent Technologies, USA). Each library was sequenced using 100 base-length read chemistry on a paired-end ow cell on the Illumina HiSeq2000 (Illumina, USA).

Genomic DNA library coverage is summarized in Supplementary Table ST13.

Genome pro le and assembly K-mers were counted with jelly sh1 118 to draw the genomic pro le. Illumina paired-end and mate-pair reads were assembled in contigs using spades v. 3.5.02 119 with kmer = 127. Contigs larger than 500 bp were assembled in scaffolds with sspace2.03 120 and gaps were lled with gap ller4 121 Scaffolds longer than 2 kb were kept for genome annotation.

Genome repeats and transposable elements

A strategy was developed to predict and then mask repetitive DNA before gene annotation. Due to the characteristics of the obtained assembly (~ 943 Mbps; 6640 scaffolds), Transposable Elements (TEs) were de novo annotated with a method speci cally designed to overcome computational challenges of TE discovery in large genomes and using an automated process, as implemented in the package REPET V2.5 122,123 . The process started with focusing on a subset of the genome su cient to identify "easy-tond" TEs. In brief, we eliminated stretches of ns > 11 from the scaffolds and kept contigs larger than 16 kbps. The resulting subset contained ~ 655 Mbps (representing 69.44% of the initial genome) and 16475 contigs, further used as effective DNA sequences to build TE consensus without N stretches. The de novo detection was rst self-aligned onto the genome subset and blasted against itself using the TEdenovo pipeline similarity approach implemented in REPET and based on Blaster default parameters. Matching sequences were clustered using programs speci c to interspersed repeats (Grouper, Recon, Piler; see REPET package). A consensus sequence was then built for each cluster, based on the multiple alignment of all members of that cluster. All consensus sequences were then classi ed using the classi cation system of Wicker et al. 124 , as implemented in PASTEClassi er tool in REPET package.

Redundancy was removed and TEs were grouped by family, resulting in a ltered library of 5085 nonredundant consensus sequences representative of "easy-to-nd" TEs. TE annotation was performed using an iterative process with REPET (Blaster, Censor, RepeatMasker, see REPET package). First, the complete "easy-to-nd" TE library was aligned by similarity searching on the initial genome (943 Mbps), obtaining a TE coverage of 42.6%. This led to the retrieval of 23211 TE full length copies (FLC, i.e., TE annotations aligned on more than 95% of their cognate consensus). From these, we used the corresponding consensus sequences and built a FLC TE sub-library made of 2643 sequences. Next, the FLC TE library was aligned against the initial genome (943 Mbps). Despite lower coverage (37.89% of the assembly), this alignment identi ed a larger number of FLCs (24190) than did the previous one. Repeats detected ab initio and annotated TEs, representing 326 Mbps in total, were masked on the genome for further annotation.

The genome was rst self-aligned and blasted against itself (similarity approach, using Blaster default parameters). Matches were clustered using programs speci c to interspersed repeats (Grouper, Recon, Piler). A consensus sequence was built for each cluster, based on the multiple alignment of all sequences thereof. All consensus sequences were classi ed using the classi cation system of Wicker et al. 124 , as implemented in PASTEClassi er. Redundancy was removed and TEs were grouped by family, resulting in a ltered library of 5085 non-redundant consensus sequences (easy-to-nd TEs). TE annotation was performed using an iterative process with TEannot pipeline (Blaster, Censor, RepeatMasker, see REPET package). First, the complete library of easy-to-nd TEs was aligned by similarity searching on the initial genome (943 Mb), with a coverage of 42.6%. This led to the retrieval of 23,211 full length copies (FLC, i.e., TE annotations aligned on more than 95% of their cognate consensus). Next, the sub-library of FLC consensus sequences (n = 2,643) was aligned against the initial genome (943 Mb). Despite lower coverage (37.89%), this alignment identi ed a larger number of FLCs (24,190) than the previous one. Summary statistics on TE annotation are given in supplementary Tables ST14 andST15.

Repeats detected ab initio and annotated TEs, representing 326 Mb in total (34.6% of the assembly), were masked on the genome for further annotation. Gene annotation was based on the search for exon/intron structure on the genome assembly and using both ab initio (using SNAP) and evidence driven (RNAseq) gene prediction.

Genome annotation

The general annotation work ow is summarised on Supplementary Figure SF4. Following genome masking (described above), gene prediction was performed using both ab initio and evidence-driven (RNA-Seq and conserved proteins) methods. Transcript contigs (obtained as described above) were rst aligned against the repeat-masked genome with BLAT. For each contig, we selected a single alignment with the highest BLAT score having at least an identity percent of 90%. Additionally, any alignment containing an intron larger than 100Kb was divided into two sub-alignments. The BLAT alignments de ned genomic regions where we performed spliced alignments with the corresponding RNA-Seq contig to accurately identify splice sites. These alignments were performed using est2genome 125 and those with at least 90% identity and a length ratio of aligned contig ≥ 85% were selected.

We used error-corrected cDNA nanopore reads generated as described above to calibrate the SNAP ab initio gene prediction tool 126 . First, ORFs were detected in corrected nanopore reads using TransDecoder (default parameters, build r20131117), and these were mapped onto the genome assembly with the same work ow used to align RNA-Seq contigs. For each genomic locus, we selected the longest nanopore read as the representative isoform, and all the representative isoforms were used to train SNAP and to generate a L. stagnalis con guration le. Finally, gene structures were predicted on the masked genome assembly using SNAP and the newly created parameter le.

Gastropod protein sequences extracted from Uniprot were mapped on the masked genome using a strategy similar to the one applied to RNA-Seq contigs. Protein sequences were rst masked for low complexity using segmasker 127 and aligned using BLAT. The best match and matches with a score ≥ 90% of the best match score, were retained and alignments containing introns ≥ 100kb were split to create distinct alignments. Alignments were then re ned using Genewise (default parameters), which is more precise for intron/exon boundary detection 128 . Alignments were selected if more than 80% of the length of the protein was aligned to the genome. In order to catch poorly conserved genes, we also relaunched the same work ow in which the rst alignment step was performed using BLASTx instead of BLAT and only the 5 best matches per genomic locus were retained.

Gene model prediction was performed using Gmove 129 . This annotation tool combines various evidence from RNAseq and conserved protein alignments, as well as ab initio predictions. It produces a consensus on gene structure without any prior calibration. Gmove creates an oriented graph (with exons as nodes and introns as vertex) from which all possible paths are extracted, representing potential gene models.

Coding frames consistent with the available protein evidence were searched and in cases where no protein alignment exists, the longest coding frame is selected. Single-exon genes and genes located in the intron of another gene were rejected if their coding sequence was shorter than 300 nt. Then, genes that overlap TEs on 50% of their lengths were also removed. Among the remaining genes, we retained all genes with a coding sequence greater than 300 nt as well as shorter genes with at least one match to UniProt proteins (e-value < 10-10).

We analysed the completeness of the predicted gene catalogue using BUSCO software (v5.2.2). The odb10 eukaryotic BUSCO dataset contains single-copy orthologous genes found in > 90% of known eukaryotic species' genomes 130 . Annotation of the predicted genes were done via BLAST (E-value threshold of 1E-5) search against Nr, SwissProt and KOG databases and search against InterPro using RunIprScan (http://michaelrthon.com/runiprscan/), a command line utility for interfacing with the InterProScan server at the European Bioinformatics Institute. Protein domain annotation was completed with the Conserved Domain Database (NCBI CD search tool).

Orthology and gene family analysis

An orthology analysis of the predicted proteome of L. stagnalis was performed using OrthoFinder at two phylogenetic levels. First, as a broad level of investigation, a set of 19 speci c proteomes was chosen from Genbank to represent a large diversity of metazoan lineages: Placozoa (Trichoplax adhaerens), Porifera (Amphimedon queenslandica), Platyhelminthes (Schistosoma mansoni), Mollusca (Biomphalaria glabrata, Aplysia californica, Lottia gigantea, Crassostrea gigas, Octopus bimaculoides), Annelida (Helobdella robusta), Nematoda (Caenorhabditis elegans), Arthropoda (Daphnia pulex, Hyalella azteca, Folsomia candida, Drosophila melanogaster), Echinodermata (Apostichopus japonicus), Hemichordata (Saccoglossus kowalevskii), and Chordata (Branchiostoma oridae, Ciona intestinalis, Homo sapiens). Orthogroups were identi ed at this scale of diversity. Gene family expansion was statistically assessed in L. stagnalis using a Z-score based test 131,132 .

Next, a re ned analysis was performed at the level of Gastropoda, with a speci c focus on GPCRs, to test for the possible implication of FMRFamide receptor expansion (as inferred from the rst analysis) in the evolution of simultaneous hermaphroditism in this class of molluscs. The analysis was based on 14 gastropod species, including 6 species with separate sexes (Haliotis rubra, Gigantopelta aegis, Batillaria attramentaria, Pomacea canaliculata, Potamopyrgus antipodarum, and Littorina saxatilis), one sequential hermaphrodite (Lottia gigantea), and 6 simultaneous hermaphrodites (Aplysia californica, Candidula unifasciata, Elysia chlorotica, Bulinus truncatus, Biomphalaria glabrata, Radix auricularia, and Lymnaea stagnalis). From this analysis, the gene contribution of each species to all 89 orthogroups with L. stagnalis members annotated as "FMRFamide_receptor-like" was summed up and used to compare the mean number of genes per species across reproductive modes (simultaneous hermaphroditism vs other; hermaphroditism vs separate sexes) using an analysis of variance.

(see supplementary information, orthology analyses).

Mitochondrial

Sequence similarity searches of the genome assembly identi ed a scaffold containing the complete mitochondrial genome of L. stagnalis. Constituent NaS Nanopore reads were checked to evaluate contiguity of the scaffold. The initial automated annotation of the mitogenome (MITOS; http://mitos.bioinf.uni-leipzig.de/index.py; 133 was manually checked using Snapgene software (www.snapgene.com) and applying criteria for annotation of molluscan mitogenomes as recommended by Ghiselli et al. 134 . The gene order of the mitogenome was compared among L. stagnalis and representative hygrophylid gastropods available in GenBank.

Whole mount in-situ hybridisation A total of 10 genes were selected for whole mount in situ hybridisation (WMISH) against male (prostate gland) and female (albumen) reproductive organs. We also surveyed late developmental stages for all genes, and also gastrulation stages for the two 2DBD-NRs. These included two DNA-binding-domain containing nuclear receptors (GSLYST00001311001 and GSLYST00003556001), 4 GPCR genes (GSLYST00019397001, GSLYST00001331001, GSLYST00016909001, GSLYST00019383001) and 3 nuclear receptors (GSLYST00008490001, GSLYST00009627001 and GSLYST00020139001), and as a positive control the seminal uid protein ovipostatin (GSLYST00000615001; [START_REF] Swart | Temporal expression pro le of an accessory-gland protein that is transferred via the seminal uid of the simultaneous hermaphrodite Lymnaea stagnalis[END_REF] . Fragments of these genes were ampli ed from cDNA using previously described standard procedures 135 and the primers described in Supplementary Table ST16. PCR fragments were cloned into a TA vector (pGEM-T Easy) and veri ed via Sanger sequencing. DIG labelled anti-sense riboprobes were synthesised as previously described 135 . A para n Tissue Microarray was 3D printed and constructed using STL les as described in Pazaitis and Kaiser 136 . Snails for these experiments were maintained as previously described [START_REF] Herlitze | Molecular modularity and asymmetry of the molluscan mantle revealed by a gene expression atlas[END_REF] . The albumen and prostate glands of mature L. stagnalis were dissected out as previously described 137 and xed for 2 hours in 0.1% glutaraldehyde in PBS with 0.1% Triton at room temperature. These tissues were then washed twice with PBS + 0.1% Tween-20, and then dehydrated through a series of ethanol with an overnight incubation in 100% ethanol. The following day the tissue was exchanged into 100% xylene and incubated for a further 24 hours before being perfused with molten para n overnight. Tissue samples were placed into the recipient para n block and allowed to set for several hours to overnight (the pairing of albumen and prostate glands from individual snails was always maintained). Microtome sections (10 µm) were collected onto coated microscope slides and allowed to adhere at 37°C overnight. Following depara nisation in xylene and a rehydration series, in situ hybridisation was performed as previously described 135 .

Snails and experimental conditions for mating observations reproductive organ RNAseq

Ten adult snails of three months old and with a shell length of 30 mm were collected from an agesynchronised mass rearing culture at Vrije Universiteit (VU) in Amsterdam. These snails were isolated for 8 days by placing each individual in a transparent, plastic container with a volume of 625 mL lled up to 460 mL. These containers were perforated with slits on opposing sides to let water ow through and were placed in a laminar ow tank with running, low-copper water. The water temperature was 20 ± 1°C and the light conditions were 12h:12h, as per standard culturing conditions at VU Amsterdam. These snails were each fed one disk of broad-leaf lettuce, measuring approximately 19 cm 2 , daily. On day seven the containers were replaced with clean ones to induce egg laying, so that this behaviour would not interfere during the mating observations on the next day. On day eight, snails were randomly placed in pairs in the same type of containers, but without slits, that were lled up to 460 mL with water of the same temperature. These containers were places on a lab bench for mating observations 6h, e.g. 138 . Once unilateral copulation was completed, the sperm donor (male-acting individual) and the sperm recipient (female-acting) were euthanised by injecting 2-3 mL of 50 mM MgCl 2 through the foot (needle gauge: 0.3 mm x 13 mm). The (male) prostate gland and the (female) albumen gland of donors and recipients were dissected out and snap frozen for RNAseq and differential gene expression analyses. Reads (see earlier RNAseq description) were quality controlled using FastQC and MultiQC. Genome-guided alignments of this data was performed using HISAT2, followed by a quality control of the alignment using QualiMap. Reads that mapped to annotated genes were counted using HTseq, after which speci c receptors (NR and GPCR) were selected for further investigation.

Declarations Figures

L. stagnalis NR relationships. Phylogenetic relationships among L. stagnalis NR families: ML tree (Jones-Taylor-Thornton substitution model, Gamma distributed substitution rate with 4 categories, statistical support from 1,000 bootstraps) based on ClustalW-aligned concatenated sequences of the DBD and most conserved part of the LBD of L. stagnalis NRs and their BLAST best-hits (alignment available upon request). Family NR0 was excluded from this analysis as it could not be aligned with DBD-only and LBDonly sequences. The triangles represent collapsed clades, with the number of branches/sequences indicated in brackets (see details on the full tree provided in Supplementary Figure SF2).

Figure 2 Phylogenetic

 2 Figure 2

Figure 4 Representative

 4 Figure 4

Table 1

 1 Summary of the L. stagnalis genome assembly.Repetitive elements spanned 37.89% of the 943 Mb-genome. Using a semi-automated process (described in the materials and methods section), we annotated the L. stagnalis genome with 2,643 transposable element (TE) consensus sequences (provisionally available via INRAE data repository:(https://entrepot.recherche.data.gouv.fr/privateurl.xhtml?token=1ad8e7c4-fa43-4002-b99f-c236eae7ba01), giving 883,147 TE copies(24,190 of which are full-length copies (FLCs). In the consensus library, 1,517 TEs were classi ed as retrotransposons (Class I) and 774 as DNA transposons (Class II), 243 were not classi ed and 105 as Potential Host Genes (PHGs). Most of these consensus sequences (2,604, i.e., 99%) have one or more full-length copy, indicating that a consensus sequence is likely a reasonable estimate of the common ancestor of any one TE. The resulting annotation was used to mask the TE landscape in the genome.

	Metric	Final assembly	Final assembly
		Contigs	Scaffolds
	Assembly size 873,879,680	942,996,421
	Number	64,200	6,640
	N50	35,351	957,215
	L50	6,921	297
	N90	7,115	151,064
	L90	27,651	1,178
	Largest	315,608	7,106,496
	GC%		37.55%
	%N's		7.33%
	Repeat content

Annotation

Following repeat masking of the genome, we employed ab initio, transcriptome-and homology-based methods to predict a set of 22,499 protein-coding genes (see Table

2

for detailed gene characteristics). This number is somewhat lower than in other annotated gastropod genomes (e.g., 23,851 genes in L. gigantea, 28,683 in A. californica, 36,943 in B. glabrata), and lower than the average number of genes found in a panel of 20 representative Metazoa (see OrthoFinder analysis).

Table 2

 2 Characteristics of the L. stagnalis genome annotation.

	Parameter	value
	Assembly size	943 Mbp
	# predicted genes	22,499
	# predicted genes without intron	3,492 (15.5%)
	# exon per gene avg:med	7.59:5
	Gene size (nt) avg:med	15,299:8,502
	CDS size (nt) avg:med	1,375:993
	Intron size (nt) avg:med	1,611:604
	Complete BUSCOs	98.5% (vs, odb10 Eukaryota)
	Complete and single-copy BUSCOs	97.3%
	Complete and duplicated BUSCOs	1.2%
	Fragmented BUSCOs	1.2%
	Missing BUSCOs	0.3%
	Total BUSCO groups searched	255
	# predicted genes with domain*	19,387 (86%)
	# predicted genes with InterPro* domain* 16,622 (74%)
	# predicted genes with Pfam* domain*	15,088 (67%)

*at least one protein domain

Table 3

 3 Summary statistics for GPCR Orthology Groups, with a focus on molluscan representatives. Percentage values indicate for each orthogroup the proportion of proteins belonging either to Mollusca, Gastropoda, or to L. stagnalis. Abbreviations: OG, Orthogroup. A Z-score value (number of standard deviations from the mean) > 2.0 indicates signi cant expansion.

	OrthoGroup	InterProScan	OG	Mollusca	Gastropoda	L.
			size			stagnalis
							Z-
							score
	OG0000008	Rhodopsin_R	1,215	1085	1035	369	2.537
		(peptide,
		FMRFamide)		(89%)	(85%)	(30%)
	OG0001159	Rhodopsin_R	72*	71	71	24	2.336
		(peptide,
		FMRFamide)		(99%)	(99%)	(33%)
	OG0002276	Rhodopsin _R	49*	48	48	17	2.452
		(peptide,
		FMRFamide)		(98%)	(98%)	(35%)
	OG0002279	Rhodopsin_R	49	49	49	22	3.106
		(2 GPCRs)		(100%)	(100%)	(45%)
	OG0002601	Rhodopsin_R	45	38	22	4	< 2
		(peptide,
	Absent in	FMRFamide)
	Deuterostoma and
	Arthropoda
	OG0010705	Rhodopsin (peptide,	8	7	7	4	< 2
		FMRFamide)
	OG0000001	Rhodopsin_R	3,444	911	561	139	< 2
		(peptide)
	OG0004949	Rhodopsin_R	30	30	27	14	3.347
		(biogenic amine)
				(100%)	(90%)	(47%)
	OG0000121	Rhodopsin_R	286	268	115	29	< 2
		(peptide and
		biogenic amine)
	OG0000012	Adhesion_R	1,072	310	186	42	< 2

*in OG0001159, the single non-molluscan member was B. oridae, in OG0002276 C. elegans.

Acknowledgements

Funding Information

This work was supported by the Genoscope, the Commissariat à l'Énergie Atomique et aux Énergies Alternatives (CEA) and France Génomique (ANR-10-INBS-09-08). JMK and YN were funded by Open Competition Grant OCENW.KLEIN.062. DJJ was funded by Deutsche Forschungsgemeinschaft (DFG project #528314512).

Data availability

The inbred line of snails is available from the authors. The Illumina and nanopore sequencing data, as well as the genome assembly and gene predictions, are available in the European Nucleotide Archive (ENA) under the project PRJEB67819. Moreover, all data are freely available from the Genoscope website: www.genoscope.cns.fr/lymnaea.

Alignment available upon request. See also Supplementary Figure SF3 for detailed clade composition). Figure 3 Gastropod proteome orthology. Left: phylogenetic relationships among 14 gastropod species (topology taken from 139 ) used for gene orthology assessment (colour code: blue, simultaneous hermaphrodite; purple, sequential hermaphrodite: pink, separate sexed). Right: pairwise correlation of species contribution in terms of number of proteins to 29,428 orthogroups produced by OrthoFinder. Correlations are particularly strong within the set of euthyneuran species used in the analysis. The heat map gradient on the right indicates the level of (positive) correlation. asterisk), the mantle (white arrowhead) and the cephalic lobes (yellow arrowheads). In adult tissues ovipostatin is highly expressed in the prostate and effectively undetectable in the albumen gland as expected. Qualitative differences are also evident for many of the receptors surveyed, and hint at possible roles in maintaining these simultaneously active reproductive organs. Representative scale bars are 100 mm.

Supplementary Files

This is a list of supplementary les associated with this preprint. Click to download. SupITFSCIREP.pdf