Nuclear genome of Bulinus truncatus, an intermediate host of the carcinogenic human blood fluke Schistosoma haematobium

Some snails act as intermediate hosts (vectors) for parasitic flatworms (flukes) that cause neglected tropical diseases, such as schistosomiases. Schistosoma haematobium is a blood fluke that causes urogenital schistosomiasis and induces bladder cancer and increased risk of HIV infection. Understanding the molecular biology of the snail and its relationship with the parasite could guide development of an intervention approach that interrupts transmission. Here, we define the genome for a key intermediate host of S. haematobium—called Bulinus truncatus—and explore protein groups inferred to play an integral role in the snail’s biology and its relationship with the schistosome parasite. Bu. truncatus shared many orthologous protein groups with Biomphalaria glabrata—the key snail vector for S. mansoni which causes hepatointestinal schistosomiasis in people. Conspicuous were expansions in signalling and membrane trafficking proteins, peptidases and their inhibitors as well as gene families linked to immune response regulation, such as a large repertoire of lectin-like molecules. This work provides a sound basis for further studies of snail-parasite interactions in the search for targets to block schistosomiasis transmission. The snail Bulinus truncatus is an intermediate host of the carcinogenic human blood fluke Schistosoma haematobium. Here the authors report the genome of Bu. truncatus, explore protein groups inferred to play a role in its interaction with the schistosome parasite, and identify expansions in gene families linked to immune response regulation.

T he phylum Mollusca (molluscs) is represented by at least 70,000 species (55 families), which are essential invertebrates of terrestrial or marine ecosystems 1 . Many act as intermediate hosts (i.e., vectors) of parasites of vertebrates including humans 2 . Key representatives of the latter group (intermediate hosts) are aquatic snails that transmit parasitic trematodes (flukes) which cause some of the most chronic and destructive neglected tropical diseases (NTDs) of humans, including clonorchiasis, opisthorchiasis and schistosomiasis 3,4 . These NTDs (called trematodiases) affect~280 million people worldwide 3,4 and often cause serious complications, particularly malignant cancers, in chronically affected people 5,6 .
Although the control of such trematodiases presently relies heavily on the chemotherapeutic treatment of people with a single drug, called praziquantel 7 , the key to effective and sustainable control is to prevent or block disease transmission. One possible way of achieving this is to break the transmission cycle at the level of entry or exit of the parasite from the snail intermediate host.
Thus, understanding the fundamental biology of the intermediate host and its relationship with the parasite, particularly at the molecular level, could assist significantly in finding new methods to interrupt transmission and prevent disease.
While there have been major advances in our understanding of the biology and molecular biology of the parasitic trematodes that cause clonorchiasis, opisthorchiasis and schistosomiasis [8][9][10][11] , this is not the case for the snail intermediate hosts that transmit these diseases. Indeed, the only studies conducted to date 12,13 have reported draft genomes for Biomphalaria glabrataa key intermediate host of the blood fluke Schistosoma mansoni, which causes hepato-intestinal schistosomiasis in humans, to underpin studies of snail-schistosome interactions 12,[14][15][16][17][18] . However, curated genomic resources are lacking for the vast majority of trematode-transmitting snails 19 , which represents a major blind spot for health-related research.
As some of our work over the years has focused on Schistosoma haematobium-the causative agent of urogenital schistosomiasis, affecting 100 million people worldwide 20 -we have become particularly keen to create a multi-omic tool box for studying the biology of the key snail intermediate host-Bulinus truncatus-of this carcinogenic blood fluke and the molecular interplay between this snail host and the larval stages of the fluke that asexually replicate within it 21 .
Recently, we defined a mitochondrial genome that represents an established laboratory line (designated BRI) 22 and a transcriptome of the adult stage of this line of Bu. truncatus 23 . In the present study, we logically extend this work, and build on recent success in the assembly of relatively large genomes (400-550 Mb) of invertebrates using third-generation (i.e. long-read or longrange) technology 11,24 , to define a reference genome for Bu. truncatus using this technology, combined with secondgeneration (short-read) sequencing and advanced bioinformatics. This work will underpin fundamental studies of the biology, ecology and population genetics of Bu. truncatus and the interactions between this snail and S. haematobium, with broader implications for comparative genomic and molecular investigations of other snail/pathogen systems.

Results
Genome assembly. From a total of 128.5 Gb (~100-fold coverage) of Illumina short-read, 11 Gb of Oxford Nanopore (~10-fold coverage) and 48 Gb of Hi-C sequence data (Supplementary Table 1), a draft genome (designated Btru.v1) was assembled for Bu. truncatus; it contains 11,176 contigs and 523 scaffolds and has a total length of 1.2 Gb (N50 = 5 Mb; L50 = 68; longest contig = 36.5 Mb) and a GC content of 36.3% (Table 1). Scaffolding with Hi-C data gave~8 .3 million contacts (Supplementary Fig. 1 and Supplementary Table 2),~7.3 and 1 million of which were within and between chromosomes, respectively; the~7 million intra-chromosomal contacts were <20 kb in length. In total, 914 of 954 (95.8%) metazoan BUSCO genes were identified in this genome, indicating that the assembly represents a substantial proportion of the complete genome (Table 1).
Subsequently, ploidy was assessed using short read data representing an individual adult specimen of Bu. truncatus ( Supplementary Fig. 2). Using 21-mers, coverage-peaks were inferred using a diploid ( Supplementary Fig. 2a) or tetraploid ( Supplementary Fig. 2b) model. Although these models were not unequivocally supported by k-mer frequency coverage or readmapping to the Btru.v1 assembly ( Supplementary Fig. 2c), the "smudgeplot" did infer diploidy (A/B) with minor triploidy (AAB), with a high probability (91%) that the data matched a diploid model and a genome size estimate of~1.2 Gb (Supplementary Fig. 2d).
Genome annotation. A large portion (51.0%) of the genome (Btru.v1) of Bu. truncatus was repetitive (Supplementary Table 3) and included transposable elements (~23.7% DNA transposons and 6.1% retrotransposons). Although most DNA transposons could not be classified, hobo-activator hAT transposon-like elements were highly represented (11.4% of the assembly). Retrotransposons were mainly long terminal repeats (LTR; 4.8%) and long interspersed nuclear elements (LINEs; 1.4%). The remaining repeat content included unclassified (18.0%) or simple (2.9%) repeat elements (Supplementary Table 3). The distribution of total repetitive elements across the Btru.v1 genome was relatively even (Supplementary Table 3 and Supplementary Fig. 3). In each 500 kb genomic region of the Bu. truncatus genome, LTRs, LINEs, and DNA transposons were encoded on average 47.8, 16.5, 375.8 times, respectively, with no clear association with gene model density (Supplementary Table 3 and Supplementary  Fig. 3). The highest frequency of encoded LTRs was observed in scaffolds: HiC_scaffold_34 (n = 174 elements in 500 kb); HiC_scaffold_309 (n = 163); and HiC_scaffold_10 (n = 153). The highest frequency of encoded LINEs was observed in scaffolds:  25 , were used for the evidence-based prediction of protein-coding genes from the Btru.v1 genome. A total of 75,434 genes were predicted in the masked assembly, reflecting marked complexity. To enable the characterisation of this gene set, we defined gene clusters ( Supplementary Fig. 4a) and selected those in clusters 2 (n = 12,209 genes), 4 (n = 5238) and 5 (n = 9674) as each having salient, common features (i.e. more than one exon, GC content and complexity within predicted protein) ( Supplementary Fig. 4b). Thus, a total of 26,292 protein-encoding genes were characterised for Bu. truncatus (Table 2) and annotated ( Table 3).
The statistics for the annotated genes in the Btru.v1 genome were similar to those of Bi. glabrata (BB02 strain) 12 : mean gene lengths (11,860 Fig. 6; Table 2; Supplementary Data 1), with 14,801 genes having transcriptional support using both data sets and displaying a direct association between RNA short-read and long-read TPM values (adjusted Rsquared: 0.6589; p-value: < 0.01; Supplementary Fig. 6). Long-read data had a mean TPM value of 35.4 and a high coverage (>80%) for 16,216 genes, compared with a TPM value of 28.1 and a high coverage (>80%) for 2507 genes for short-read data.
Within the predicted gene set, we identified 905 (95%) of 954 complete, conserved metazoan genes by BUSCO analysis, suggesting that this set represents most of the genome. This result is similar to that (89%) inferred for complete BUSCOs in Bi. glabrata (BB02 strain); however, the Btru.v1 genome was predicted to have more BUSCO gene duplicates (17.6%), but fewer fragmented (2.7%) or missing (2.4%) genes ( Table 2). For Bu. truncatus, we inferred 23,248 (88.4%) or 18,139 (69.0%) genes encoding proteins with sequence homology to proteins present in the non-redundant UniProt TrEMBL or SwissProt databases, respectively (Table 3).
Excretory/secretory (ES) proteins are reported to play central roles in snail-schistosome interactions 18 . Within the Bu. truncatus gene set, 1096 (4.2%) genes were predicted to encode extracellular ES proteins, based on the presence of a signal peptide domain (for 3121 genes; 11.9%) and absence of one or more transmembrane domains (for 18,507 genes) ( Table 3; Supplementary Data 1). The full complement of ES proteins (i.e. the secretome) was inferred to represent 434 peptides with KEGG annotations that could be assigned to the following protein groups: peptidases and their inhibitors (n = 102), exosome (67), membrane trafficking (48), glycosaminoglycan-binding proteins (47) and lectins (26). Of these 434 proteins, 330 were orphans, 118 had homology to individual proteins in Bi. glabrata (BB02 strain), and 251 had transcriptional support in the adult stage of Bu. truncatus.
Marked synteny in mollusc genomes. The annotated proteincoding gene set was first used to assess conserved gene order (synteny) in Btru.v1 genomic scaffolds compared to genes encoded in the reference genomes for Bi. glabrata, Achatina immaculata (gastropod; family Achatinidae) and Pecten maximus (bivalve; family Pectinidae); the relationship among these snails is summarised in phylogenetic analyses of concatenated amino acid sequence data inferred from 2315 SCOs using maximum likelihood (ML) and Bayesian inference (BI) methods (Supplementary Fig. 7). Upon pairwise comparison, many blocks of nucleotide sequence were aligned across conserved gene regions among Bu. truncatus, Bi. glabrata (S1316-R1 strain), Ac. immaculata and P. maximus (Fig. 1), with most conservation seen between Btru.v1 and Bi. glabrata ( Fig. 2a; Table 4). Almost all genome scaffolds contained regions that aligned in 594 syntenic blocks of ≥3 single-copy orthologs spanning 557. Protein ortho-groups and expansions. Next, the annotated Bu. truncatus protein-coding gene set was compared to genome annotations available for Bi. glabrata (BB02 strain), Aplysia californica, Elysia chlorotica and P. maximus (Table 5 and Fig. 2); the relationship among these snails is summarised in phylogenetic analyses of concatenated amino acid sequence data inferred from 2315 SCOs using maximum likelihood (ML) and Bayesian inference (BI) methods ( Supplementary Fig. 7) Upon comparison, 17,866 ortho-groups containing one or more proteins were identified. Most ortho-groups (n = 13,302) were shared by Bu. truncatus (n = 20,349 proteins) and Bi. glabrata (n = 19,853 proteins);  Fig. 1 Synteny between genomes. Synteny and contiguity of the genome (Btru.v1) of Bulinus truncatus with the draft genome of Biomphalaria glabrata and the chromosomal-level reference genomes of each Achatina immaculata and Pecten maximus, respectively. Scaffolds are arranged in circular (circos) plots with reference scaffolds for Ac. immaculata or P. maximus linked to inferred syntenic blocks of the Bu. truncatus genome using distinctly-coloured bars. Linked syntenic blocks between scaffolds of Bu. truncatus and Bi. glabrata are shown (light blue). a Synteny between the reference genomes of Bu. truncatus and Bi. glabrata established based on the positions of 595 syntenic blocks each containing three or more single-copy orthologs (i.e. 4047 of a total of 11,051). b Synteny between the reference genomes of Bu. truncatus and Ac. immaculata, established based on the positions of mapped Bu. truncatus proteins in 243 syntenic blocks each containing six or more single-copy orthologs (i.e. 2234 of a total of 8517). c Synteny between the reference genome of Bu. truncatus and P. maximus established based on the position of 291 syntenic blocks containing three or more single-copy orthologs (i.e. 1139 of a total of 4685).
Protein groups inferred to be involved in the snail-schistosome relationship. Given the role of Bu. truncatus as a key intermediate host for S. haematobium and related schistosome taxa 21 , we explored protein groups inferred to be involved in the immune/defence system of the snail and/or interactions with schistosomes (  Table 5).
Proteins of the Guadeloupe resistance complex (GRC) 26 and the polymorphic transmembrane cluster 2 (PTC2) 13 , proposed to be loci associated with parasite recognition and/or reduced susceptibility to schistosome infection 13,26 , were inferred for Bu. truncatus. The 11 GRC homologs predicted ( Fig. 3; Table 6) were linked to cellular processes (lysosome and adherens junction) and metabolic processes required for glycosaminoglycan degradation, and five PTC2 homologs ( Table 6, Fig. 3 and Supplementary  Table 5) linked to metabolic processes required for glycosphingolipid biosynthesis.
The baculovirus inhibitor of apoptosis (IAP) repeat (BIR) protein homologs (n = 117) identified in Bu. truncatus were inferred to be involved in apoptosis, ubiquitin-mediated proteolysis, necroptosis and NF-kappa beta-signalling pathways ( Table 6). Most BIR proteins (n = 25) of Bu. truncatus were represented within ortho-group OG0000004 containing 43, 1, 9  Fig. 3). These findings are consistent with earlier reports showing an expansion of BIR proteins in molluscs 12 .
Likely to be involved in reduced snail susceptibility to schistosome infection 18 are cathepsin homologs (n = 21) of Bu. truncatus, commonly involved in transport and catabolism (lysosome, phagosome and autophagy), cell growth and death (apoptosis) and antigen-processing/presentation pathways ( Table 6). Ten of these cathepsins were within ortho-group OG0000133, together with 9, 5, 1 and 4 respective orthologues of Bi. glabrata, Ap. californica, E. chlorotica and P. maximus (Fig. 3). Of the chitinases (n = 103) predicted to be involved in amino sugar and nucleotide sugar metabolism in Bu. truncatus (Table 6) and reduced snail susceptibility to schistosome infection, 36 were within ortho-group OG0000085, but had no ortholog in any other mollusc species studied (Fig. 3). In addition, of the calmodulins (n = 42) that likely associate with pathogen interactions 27 (Tables 6), 6 were within ortho-group OG0000309, also with 6, 4 and 4 respective orthologues of Bi. glabrata, E. chlorotica and P. maximus (Fig. 3).

Metabolism Genetic Information Processing
Signalling and Cellular Processes * * * ** Fig. 2 Orthologous protein groups. Orthologous groups of one or more protein(s) in gastropod taxa with available gene annotation. Pecten maximus is an outgroup (bivalve). a Pairwise comparison of the orthologous groups common to Bulinus truncatus and Biomphalaria glabrata. b Genes in orthologous groups common to Bu. truncatus and Bi. glabrata. c UpSet plot of the intersections of unique or shared orthologous groups inferred from protein data sets for Bu. truncatus, Bi. glabrata, Aplysia californica, Elysia chlorotica and P. maximus. d Phylogenetic tree inferred from single copy orthologs aligned among selected gastropod taxa, including all available annotated genes. As the topology of the maximum likelihood (ML) and Bayesian inference (BI) trees was the same, the ML tree is displayed and shows nodal support values for both BI (pp) and ML (bootstrap). e Genes (annotated by KEGG) predicted to be unique to Bu. truncatus (coloured), or shared between or among the five molluscan species included here.
represented fibrinogen C-terminal domain proteins within orthogroup OG0000016 (38 proteins), also containing 8, 14 and 2 orthologs of Bi. glabrata, Ap. californica and P. maximus, respectively (Fig. 3).  (Fig. 4b, Supplementary Fig. 8). Unlike sFReDs, some Class A proteins contained N-terminal α-helices (Class A α) or β-sheets and α-helices (Class A βα). The presence of α-helices in Class A α proteins in close proximity to the N-terminus was suggestive of a signal peptide, as they often preferentially adopt an α-helical form 28 . Class B and C proteins matched Bi. glabrata FREP proteins, with a single FReD domain linked to one (Class B) or two (Class C) N-terminal immunoglobulin superfamily (IgSF) domains by α-helices (Supplementary Table 6 and Supplementary  Fig. 8). Class D proteins were Bi. glabrata FREM-like proteins (Supplementary Table 6). Class E and F were Bu. truncatus FReD-like proteins with complex structures in the N-terminus but had no structural homology to known, predicted Bi. glabrata proteins (Fig. 4b). Some (Class A-like) proteins contained a partial FBG domain (Supplementary Table 6), and, thus, had been excluded from the original phylogenetic analysis (Fig. 4a).
A detailed exploration of fibrinogen-related proteins in these phylogenetic groups based on their predicted structures provided enhanced insights. Group 1 contained 36 Bi. glabrata proteins, 21 of which had a predicted signal peptide based on primary amino acid sequence (Fig. 4a). Group 1 proteins were classified as Class B FREPs, (n = 9), Class C FREPs (n = 20) or sFReDs with (Class A α n = 4; Class A βα n = 1) or without (Class A; n = 1) additional N-terminal tertiary structures. Group 2 included 27 Bu. truncatus proteins, 16 of which had signal peptides. Most Group 2 Bu. truncatus proteins were classified as sFReDs without (Class A; n = 14) or with additional N-terminal structures (Class A α n = 4; Class A βα n = 7). Many proteins in Group 2 were assigned to ortho-group OG0000016, including one Bu. truncatus protein (Btru_033719) with a distinct N-terminal domain encoding several, ordered β-sheets (Class E). The three Bi. glabrata proteins within Group 2 all encoded a signal peptide and were classified as sFReDs (n = 2) or Class Aα (n = 1) (Fig. 4a). Group 3 contained only two Bi. glabrata proteins and were classified as Class Aα. Group 4 included 58 Bu. truncatus and 5 Bi. glabrata proteins, 33 and 3 of which had signal peptides, respectively. Proteins in Group 4 were most diverse, including FReDs assigned to ortho-groups OG0000103, OG0000581, OG0010992 and OG0015404. Bi. glabrata proteins in Group 4 were all predicted to be sFReDs and represented a distinct subgroup with three Bu. truncatus proteins (Fig. 4b; pp = 0.85). Most Bu. truncatus in Group 4 were sFReDs with (Class Aα n = 13; Class A βα n = 2) or without additional N-terminal structures (Class A; n = 37) or were predicted to encode only a partial FReD domain (n = 5). One FReD-like protein (Btru_048110) in this group encoded a novel N-terminal domain (Class F; Fig. 4b), but clustered with Bu. truncatus sFReDs (Class A). Group 5 proteins were all classified as Class Aα (4 Bu. truncatus and 2 Bi. glabrata), two of which were assigned to ortho-group OG006406. Only one protein (Bu. truncatus Btru_048819) in this group had a predicted signal peptide. Group 6 proteins all had signal peptides; two were Bi. glabrata proteins (one sFReD and one in Class C). Interestingly, the tertiary structure of a Class C-like Bu. truncatus protein (Btru_007887) in this group encoded two IgSF-like domains (cf. Figure 4b), which had not been detected previously based on amino acid sequence homology (Supplementary Data 1); this Bu. truncatus protein had a shorter FReD domain and a longer chain of α-helices than seen for Bi. glabrata Class C FREPs ( Supplementary Fig. 8). Group 7 contained three Bu. truncatus and three Bi. glabrata protein, all of which were classified as Class A βα, mostly within ortho-group OG0004042, most of which were not predicted to encode a signal peptide.

Discussion
Here, we present the first draft nuclear genome (Btru.v1) of a key representative of the snail genus Bulinusa complex of at least 37 species presently divided into four main groups 29 . This draft genome of the BRI-laboratory line of Bu. truncatus, originally from Egypt, was assembled using a combination of second-(short-read) and third-generation (long-read) sequence data. The draft assembly shared a high degree of contiguity with the available reference genomes (chromosomes) of the giant land snail, Ac. immaculata (Gastropoda: Styllommatophora) and a marine scallop, P. maximus (Bivalvia: Pectinida). Most gene models were strongly supported by short-and long-read RNA sequence data, and high-quality models enabled an exploration of the molecular biology of Bu. truncatus (BRI strain) and comparative investigations of the genomes of selected molluscs, including Bi. glabrata-which is an intermediate host of S. mansoni.
Bu. truncatus is proposed to have four sets of chromosomes (tetraploid) 29,30 , thought to have arisen via alloploidy by ancestral hybridisation of closely related diploid species 29 . This proposal is supported by the karyotype and zymograms of an Egyptian isolate of Bu. truncatus 31 . However, alternate hypotheses are that this tetraploidy might have resulted from evolutionary saltation(s) upon nuclear fusion of genomes following hybridisation of two distinct diploid species (i.e. allo-tetraploidy) or whole-genome duplication in a diploid ancestor (i.e. auto-tetraploidy). The polymorphism seen in the Btru.v1 genome of Bu. truncatus,  which seemed to reflect a diploid organism ( Supplementary  Fig. 2), was consistent with that observed in Ac. immaculata which, indeed, underwent whole-genome duplication in a diploid ancestor 32 . This evidence indicates that Bu. truncatus is an autotetraploid snail, with limited chromosomal divergence. Limited genetic divergence among the four sets of chromosomes and/or the use of a single restriction enzyme (DpnII) for Hi-C library construction, are probable reasons for the Bu. truncatus scaffolds being shorter than the expected chromosome lengths. Highcoverage long-range and long-read genome sequencing and comparative karyotypic studies of key members of the Bulinus complex, using Btru.v1 as a reference, should establish their chromosomal evolution. A genome-wide analysis revealed more genomic synteny between Bu. truncatus and Bi. glabrata 1316-R1 strain 13 (snail hosts of schistosomes) than between Bu. truncatus and Ac. immaculata (land snail) and P. maximus (marine scallop). This finding is consistent with the evolutionary relationships of planorbid/bulinid gastropods 33 . The synteny of large regions (>46% of aligned genomes) of the Bu. truncatus and Bi. glabrata genomes has loci likely central to the susceptibility of these snails to respective schistosome species 14 Fig. 3 Gene expansions in key protein groups. Expansion of protein ortho-groups in Bulinus truncatus predicted to relate to snail-schistosome interactions, based on published information (see Table 6). Cluster dendograms showing orthogroups for Bu. truncatus with the largest gene expansions inferred from protein data sets of Bu. truncatus, Biomphalaria glabrata, Aplysia californica, Elysia chlorotica and Pecten maximus as defined using OrthoFinder. The locations of genes linked to the Guadeloupe resistance complex (GRC) or polymorphic transmembrane cluster 2 (PTC2) are indicated (boxed). NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-022-28634-9 ARTICLE NATURE COMMUNICATIONS | (2022) 13:977 | https://doi.org/10.1038/s41467-022-28634-9 | www.nature.com/naturecommunications homologous to loci encoding members of the GRC and PTC2 complexes associated with reduced susceptibility of Bi. glabrata to S. mansoni infection were identified 13,34 . Detailed comparative studies are now needed to fill the void in our knowledge of the molecular basis of susceptibility in Bu. truncatus to S. haematobium.
More broadly, a large proportion of proteins inferred for Bu. truncatus were orthologous to those predicted for Bi. glabrata (BB02 strain) and formed ortho-groups to the exclusion of proteins of the other species of gastropod (Ap. californica and E. chlorotica) or bivalve (P. maximus) studied. These distinct planorbid/bulinid protein groups may be linked to one or more evolutionary events which led to the speciation within the superfamily Lymnaeoidea. Expansions of some protein-gene families/groups in Bu. truncatus were relatively consistent with those in Bi. glabrata 12,14,[35][36][37] . For instance, shared expansions were seen for some genes encoding GPCRs, proteins involved in membrane trafficking, and peptidases and their inhibitors. The large repertoire of genes encoding GPCRs and transport proteins likely relates to the diversity of receptors in olfactory organs of gastropods with key chemosensory roles 12,38 .
Key expansions of select gene families could have arisen during the evolution of snail defences against pathogens, including . Proteins with predicted signal peptide (SP) domains are indicated with cyan boxes. A distinct colour represents each FReD class, with gene accession numbers boxed. # Bu. truncatus ortho-groups inferred from proteins of Bu. truncatus, Bi. glabrata, Aplysia californica, Elysia chlorotica and Pecten maximus (cf. Figure 3). b Tertiary structure models for FReDs of Bu. truncatus employed for enhanced classification (A-F). C-terminal domain (*), as well as IgSF-like and α-helix and/or β-sheet structures, are shown.
schistosomes. For example, the cathepsin L-like genes of Bu. truncatus have homologs in Bi. glabrata which are known to regulate snail immunity 39 . The expansion of Bu. truncatus chitinase-like proteins (n = 103) is consistent with those in cephalopods 40 and bivalves 41 , in which there has been a marked gene expansion relating to the regulation of immune responses 41 . In Bi. glabrata, there is evidence that a chitinase-like protein is associated with one or more loci that modulate susceptibility to S. mansoni 26 . Also carbohydrate-binding C-type and FREPs/FReD-like lectins could play key roles in the defence of snails against schistosomes 36,[42][43][44][45][46] . We identified a significant expansion of genes encoding FReD-like lectins in Bu. truncatus, similar to that seen in Bi. glabrata. However, the protein members of these species were very distinct, with few homologs clustering together in orthogroups. This finding is consistent with a previous study of Bu. truncatus, in which no canonical Class C and/or B IgSF-FREPs were identified 23 . Detailed comparative analysis of FREPs/ FReD-like lectins, guided by tertiary structure models, allowed the identification/classification in Bu. truncatus of a single IgSF-FREP-like lectin (Class C) and no gene expansion event in this group, in accord with studies of the pond snail, Physella acuta 47 , the common periwinkle, Littorina littorea 48 and the sea hare, Ap. californica 49 . Interestingly, there was a marked expansion in the number and diversity of Class A sFReD-like proteins in Bu. truncatus, with two large groups with a similar diversification to FREPs in Bi glabrata 36 . The presence of additional α-helices and β-sheets in the N-terminal regions of Bu. truncatus sFReDs and signal peptide-like domains in some proteins (Class A α) based on tertiary (but not primary) structure models suggest marked functional diversity, which we propose is central to essential immunobiological processes in Bu. truncatus. Clearly, further work is needed to establish the roles of sFReDs in Bu. truncatus and draw comparisons with information available for Bi. glabrata 46 , Ph. acuta 47 , Mytilus edulis 50 and other molluscs.
The draft genome (Btru.v1) of Bu. truncatus encodes many BIRs/IAPs, calmodulins and Toll-/IL-1-related proteins which have conserved orthologues/paralogues in Bi. glabrata. An expansion of BIRs/IAPs was reported earlier for Bi. glabrata 12 . While the exact role(s) of these molecules is/are not yet understood, the gastropod and bivalve species studied here encode several copies of genes that might relate to a regulatory role in apoptosis and innate immune response, with an observed gene expansion in snails that act as intermediate hosts for schistosomes 12 . Calmodulins transduce signals in response to increases in intracellular Ca 2+ , represent a major component of calcium-dependent signalling pathways 51 and can play a role in pathogen defence 52 . The diversification of calmodulins in molluscs has been reported previously 53 and has been associated with defence against S. mansoni infection in Bi. glabrata 27 , and against bacteria and yeast 54 .
Most invertebrate immune systems include an array of Tolllike receptors (TLRs) that mediate TLR-directed innate immunity to a wide range of pathogen species 55,56 . A much larger number of TLR-like proteins was predicted here from the Btru.v1 genome than reported previously for the transcriptome 23 , and more than two-times the number predicted from the Bi. glabrata genome 12 . Most of the TLRs predicted for Bu. truncatus were categorised (KEGG BRITE) as conserved TLR3/4-like molecules, which are likely to be conserved in invertebrates and mammals 57 . Future work is needed to acurately classify TLRs of planorbid/bulinid snails, to localise their expression in cells and tissues, and to establish which pathogen-associated molecular patterns (PAMPs) they recognise.
The relative conservation of the order of genes in the genomes of the snails studied here indicates that it should be feasible to characterise the genomes and gene orthologues of a range of lymnaeid and physid snails, which are key intermediate hosts (vectors) of socioeconomically important parasitic trematodes other than schistosomes 2using the same technological approach as established here. Such an effort could assist in closing some of the knowledge gaps that exist in the understanding of systematics of these groups 22,33 , and would provide insight into the molecular evolution of molluscan groups.
The availability of the laboratory (BRI) lines for both Bu. truncatus and S. haematobium offers excellent opportunities to now studyunder well-controlled experimental conditionsthe molecular biology of each of these two invertebrates as well as their interactions. Bu. truncatus is an essential intermediate host of S. haematobium. In an aquatic environment, this snail becomes infected by the miracidium of S. haematobium; the ciliated ectoderm of this first larval stage sloughs off upon entry through the snail foot; and the miracidium transforms into a sporocyst, which then undergoes extensive asexual replication within the snail host. Having high-quality genomes and transcriptomes for both Bu. truncatus and S. haematobium now underpins critical molecular investigations of this asexual phase of reproduction, the cross-talk that occurs between the parasite and the snail host during replication, and the mechanisms and/or processes that govern snail susceptibility and/or immunity to the parasite. We propose that the use of a multi-omics approach 58 , involving the use of genome-guided transcriptomic, proteomic, lipidomic and/ or glycomic analyses as well as high-resolution single-cell and spatial transcriptomics 59 , will strongly complement this focus. In addition, explorations of tertiary structure models for all proteins encoded in snail genomes using AlphaFold 60 should allow the identification of distant orthologues and elucidate dark matter in these proteomes.
In conclusion, defining the first nuclear genome (of ≥ 1 Gb) for a well-defined laboratory line of Bu. truncatus opens the door to exploring a range of operational taxonomic units (OTUs) of Bulinus from natural populations as well as other key species of snail hosts of schistosomes, as a basis for future systematic, genetic, epidemiological and ecological investigations. Insights into these areas could significantly assist both fundamental and applied studies of schistosomes and schistosomiases, and enable the development of new interventions for this important neglected tropical disease-complex. Isolation of high molecular weight genomic DNA, library construction and sequencing. High quality genomic DNA was isolated from two adult Bu. truncatus snails using the Circulomics Tissue Kit (Circulomics, Baltimore, MD, USA). The integrity of the DNA was assessed using Genomic DNA ScreenTape and the Agilent 4200 TapeStation (ThermoFisher, MA, USA). Low molecular weight DNA was removed using a 10 kb short-read eliminator (SRE) kit (Circulomics, Baltimore, MD, USA). The high molecular weight DNA from the two individual snails was used to construct Nanopore Rapid Sequencing (SQK-RAPD004; Oxford Nanopore Technologies) and Ligation Sequencing (SQK-LSK109; Oxford Nanopore Technologies) genomic DNA libraries, according to the manufacturer's instructions. Each flow cell used was washed using a Flow Cell Wash Kit (EXP-WSH003; Oxford Nanopore Technologies, Oxford, UK) and re-used to sequence additional SQK-LSK109 libraries. All libraries were sequenced in the MinION sequencer (Oxford Nanopore Technologies). Following sequencing, bases were converted into FASTQ format from raw FAST5 signals using the program Guppy v.4.2.2 (Oxford Nanopore Technologies). Reads with an average quality (Q) value of <7 were removed. A short-insert (500 bp) genomic DNA library was also constructed using the DNA from one snail, and paired-end sequenced (150 base reads) using TruSeq sequencing chemistry and the NovaSeq sequencing platform Isolation of total RNA, Oxford nanopore library construction and sequencing. Total RNA that had been isolated from an adult of Bu. truncatus (Egyptian strain) in an earlier study 23 was used to prepare a long-read library using the Oxford Nanopore PCR-cDNA Sequencing Kit (SQK-PCS109; Oxford Nanopore Technologies, Oxford, UK), as recommended. This library was sequenced on a MinION sequencer (Oxford Nanopore Technologies) using an EXP-FLP002 flow cell priming kit and a R9.4.1 flow cell (FLO-MIN106). Following sequencing, bases were converted into the FASTQ format from raw FAST5 signals using the program Guppy v. 4.2.2. Reads with an average Q value of <7 were removed.

Methods
Assembly of the genome. Long reads from the genomic DNA from two Bu. truncatus snails were used to assemble contigs using FLYE v2.8-b1674 62 with the --nano-raw option and setting a genome size estimate of 900 Mb. Errors in long read sequence data were initially corrected using medaka_consensus in the medaka package v.1.0.3 (https://github.com/nanoporetech/medaka) and nanopore read data. Contigs were polished with the data derived from the short-insert (500 bp) genomic DNA library using pilon v.1.23 63 . Scaffolds were combined with the in situ Hi-C data using 3D-DNA v.180922 64 . Haplotig redundancy was removed by mapping long-read data to the genomic scaffolds, and haplotigs were eliminated employing purge_haplotigs v.1.1.1 65 . Error-corrected long reads were then used to close gaps in scaffolds using TGS-GapCloser v.1.1.1 (https://github.com/BGI-Qingdao/TGS-GapCloser). Following the mapping of short-read data to the genomic scaffolds, haplotig redundancy was eliminated using purge_haplotigs v.1.1.1 65 . The completeness of the genome was assessed (in genome-mode) using BUSCO v 4.0.2 66 .
Assessing genome size, heterozygosity and ploidy. A short-insert (500 bp) genomic DNA library from a single Bu. truncatus snail was used to estimate genome size, heterozygosity and ploidy using the GenomeScope v.2.0 and smudgeplot v.0.2.4 packages 67 . Input into each program was the frequency of 21-mers in the raw short-read data determined using kmc v.3.1.1 68 . Upper and lower frequencies used in smudgeplot were 1300 and 28, respectively. GenomeScope analyses were performed assuming a diploid or tetraploid genome model, according prior evidence for Bulinus from the literature 29 . Reads from the short-insert genomic DNA library were also mapped to the reference genome using bwa v.2 69 to estimate minor allelic frequencies and ploidy using PloiPy https://github.com/ floutt/PloidPy.
Predicting repeat-elements, gene models and protein function. Repeat elements in the genome were predicted using RepeatModeler v. 1.0.8 (http:// www.repeatmasker.org) and EDTA v v.1.9.4 70 . Libraries were combined and redundancy was removed using CD-HIT v.4.8.1 71 . The final repeat element library was used to mask the genome using RepeatMasker v.4.1.1 72 . Gene models were predicted using funannotate v.1.7.4 (https://github.com/nextgenusfs/funannotate), publicly available RNA-seq data (NCBI BioProject PRJNA680620) 23 (Supplementary Table 1) and inferred protein sequence data sets for Bi. glabrata (BB02 strain) 12 76 : 3. The completeness of the gene set was assessed (in protein-mode) using the tool BUSCO v 4.0.2 66 . The annotation of each inferred amino acid sequence was achieved using InterPro v5. 35 77 , EggNOG mapper v.5.0 78 and/or homology (E-value ≤ 10 −5 ) to proteins in the databases Swiss-Prot and TrEMBL within UniProtKB (accessed 20 December 2020) 25 , Kyoto Encyclopedia of Genes and Genomes (KEGG) 79 and/or MEROPS release 12 80 using DIAMOND BLASTp v. 0.9.21 81 . Protein groups and pathways were inferred based on homology to KEGG orthology (KO) terms linked to curated KEGG BRITE and pathway hierarchies. Signal peptide domains and/or transmembrane domains were predicted using phobius v.1.04 82 . The sub-cellular localisation of protein sequences was predicted computationally using the program MultiLoc2 v.2.2.25 83 . Evidence of gene transcription was inferred by mapping short and long RNA-seq data to the genome using HISAT2 v.2.1.0 84 , and the level of transcription per gene (in transcripts per million, TPM) was inferred using StringTie v2.1.2 75 . Gene models were inferred to have transcriptional support if one or more library had a TPM value of >0.2.
Protein-encoding genes were retained based on the features of their gene models. For each gene, the following features were curated: (1) GeneValidator v.2.1.10 85 score estimated using comparisons to proteins in Swiss-Prot within UniProtKB 25 ; (2) evidence of transcription (in TPM); (3) sequence homology (E-value ≤ 10 −5 ) to proteins in the TrEMBL within UniProtKB 25 ; (4) proportion of proteins predicted to be a "low probability subsequence" (LPS) using the program fLPS v.1 86 ; (5) GC content for coding domain; (6) length of inferred mRNA sequences; (7) number of exons in the gene model; (8) presence of one or more Pfam conserved domains inferred using InterProScan; and (9) numbers of genes representing individual groups of orthologous protein between Bu. truncatus (Btru.v1) and Bi. glabrata 12 , established using OrthoFinder v.2.3.11 87 . These features were transformed to normal distributions and subjected to PCA and K-means clustering analyses. Clusters with protein-encoding-like genes were retained for further curation and characterisation. Subsequently, we studied the distribution of repeats in the genome and their association with the 5' and/or 3' untranslated regions (UTRs; 5000 nucleotides for each) of protein-coding genes, employing the Fisher's exact test (pvalue < 0.01) to assess statistical significance of association(s).
Subsequently, single-copy orthologues were inferred from homologous genes shared by Bu. truncatus, Bi. glabrata, E. chlorotica, Ap. californica, P. maximus (= outgroup) and/or Ac. immaculata, and their amino acid sequences conceptually translated. The clusters of single-copy orthologues representing all five or six species were aligned using the program AQUA 92 , employing the programs MUSCLE v3.8.31 93 and MAFFT v.7.271 94 for the alignment and RASCAL v1.34 95 for alignment refinement. Individual clusters of genes with an alignment score of ≥0.8, obtained from the program NorMD 96 , were merged using the program PartitionFinder v2.1.1 97 to assign each merged partition to a predicted amino acid substitution matrix. Partitions that did not contain more than 20 amino acids were removed. Remaining partitions were then subjected to separate phylogenetic analyses using the Bayesian inference (BI) and maximum likelihood (ML) treebuilding methods. BI analysis was conducted using the program MrBayes v3.2.6 98 from four independent Markov chains, run for 1,000,000 metropolis-coupled MCMC iterations, where trees were sampled every 1000 iterations. The resultant tree was inferred by initially discarding 25% of sampled trees as burn-in, and then using the remaining trees to infer tree topology, branch lengths and to calculate Bayesian posterior probabilities (pp). For ML, a partitioned ML tree was constructed using the program RAxML v8.2.6 99selected models for each partition, employing 20 iterations-and the best tree selected for bootstrap analysis (n = 100). A representative tree was prepared using FigTree v.1.31 (http:// tree.bio.ed.ac.uk/software/figtree).
Identification and characterisation of fibrinogen-related domain containing proteins. Within the Bu. truncatus and Bi. glabrata protein data sets, proteins with sequence homology to the fibrinogen beta and gamma chains, C-terminal globular domain (Pfam: PF00147.20) were identified using hmmsearch (HMMER v.3.2.1; http://hmmer.janelia.org/) and using a threshold of E-value ≤ 10 −5 . Proteins containing a conserved fibrinogen domain were then aligned using hmmalign and using the --trim option and sequences with more than 150 amino acid residues across the conserved fibrinogen domain were retained. Trimmed sequences were de-gapped and re-aligned using MUSCLE v3.8.31 93 . Aligned sequences were then subjected to Bayesian inference (BI) analysis using the program MrBayes v3.2.6 98 and using a WAG model with fixed rate matrices, generating 4,000,000 trees and sampling every 400th tree. The resultant tree was inferred by initially discarding 50% of sampled trees as burn-in, and then using the remaining trees to infer tree topology, branch lengths and to calculate Bayesian posterior probabilities (pp). Phylogenetic trees were rendered and annotated using ggtree (v.1.10.5) 100 in R v.3.4.3 (http://www.R-project.org/). Next, mature peptides (excluding predicted signal peptide domains) were subjected to tertiary structure prediction using AlphaFold 60 . The best ranked model was retained and used in subsequent analyses. FReD-like proteins were classified based on known domains described previously 36 , including the presence of a complete or partial C-terminus fibrinogen domain (Class A or A-like, respectively) and with one (Class B; 1-IgSF FREP) or two (Class C) additional N-terminal immunoglobulin superfamily (IgSF) domain. Proteins with FReD domain (Class A) but with additional N-terminal alpha-helices and/or beta-sheets or with complex novel domains were also classified.
Reporting summary. Further information on research design is available in the Nature Research Reporting Summary linked to this article.

Data availability
The nucleotide sequence data linked to the nuclear genome reported in this article is publicly available in the GenBank database and the Sequence Read Archive (SRA; https://www.ncbi.nlm.nih.gov/sra) under the accession numbers SAMN17050146, SAMN16898649 and SAMN16898648 with the NCBI BioProject accession number PRJNA680620. Protein sequences used for sequence homology searches are available from the Swiss-Prot (UniProtKB; https://www.uniprot.org/help/uniprotkb; accessed 20 December 2020) 25 , TrEMBL (UniProtKB; https://www.uniprot.org/help/uniprotkb; accessed 20 December 2020) 25 , Encyclopedia of Genes and Genomes (KEGG; https:// www.genome.jp/kegg/; accessed 20 December 2020) 79 and MEROPS release 12 (https:// www.ebi.ac.uk/merops/) 80 databases. All other data used are referred to in this article and its supplementary files.