Mitochondrial genome of Bulinus truncatus (Gastropoda: Lymnaeoidea): Implications for snail systematics and schistosome epidemiology.

Many freshwater snails of the genus Bulinus act as intermediate hosts in the life-cycles of schistosomes in Africa and adjacent regions. Currently, 37 species of Bulinus representing four groups are recognised. The mitochondrial cytochrome c oxidase subunit 1 (cox1) gene has shown utility for identifying and differentiating Bulinus species and groups, but taxonomic relationships based on genetic data are not entirely consistent with those inferred using morphological and biological features. To underpin future systematic studies of members of the genus, we characterised here the mitochondrial genome of Bulinus truncatus (from a defined laboratory strain) using a combined second- and third-generation sequencing and informatics approach, enabling taxonomic comparisons with other planorbid snails for which mitochondrial (mt) genomes were available. Analyses showed consistency in gene order and length among mitochondrial genomes of representative planorbid snails, with the lowest and highest nucleotide diversities being in the cytochrome c oxidase and nicotinamide dehydrogenase subunit genes, respectively. This first mt genome for a representative of the genus Bulinus should provide a useful resource for future investigations of the systematics, population genetics, epidemiology and/or ecology of Bulinus and related snails. The sequencing and informatic workflow employed here should find broad applicability to a range of other snail intermediate hosts of parasitic trematodes.

Introduction

Human schistosomiasis represents a group of neglected tropical diseases (NTDs) caused by blood flukes of the genus Schistosoma (schistosomes). More than 200 million people worldwide are affected by the hepato-intestinal or urogenital form of schistosomiasis. Notable species within the S. haematobium-group include S. haematobium, Schistosoma intercalatum and Schistosoma guineensis (human pathogens), and Schistosoma bovis, Schistosoma curassoni, Schistosoma mattheei, Schistosoma leiperi and Schistosoma margrebowiei (pathogens infecting domestic and wild ruminants). No vaccines are available to prevent schistosome infections or disease, and the reliance on praziquantel – an anthelmintic drug which is not highly effective against all developmental stages of schistosomes (Gryseels et al., 2006) – for the mass treatment of humans against schistosomiasis, carries a significant risk of treatment failures, and the excessive treatment over time has a risk of resistance against this drug emerging in schistosomes (Greenberg, 2013). For these reasons, the effective prevention and control should rely on a sound understanding of the life-cycle and transmission of schistosomiasis for the implementation of integrated prevention/control, focused on the disruption of transmission to humans (World Health Organization, 2020).

As some species of freshwater snails of the genus Bulinus act as intermediate hosts in the life-cycles of schistosomes in Africa, Madagascar and adjacent regions, knowing which species of Bulinus transmits which species of Schistosoma is crucial. Therefore, the accurate identification of species of Bulinus is central to elucidating and understanding the epidemiology of schistosomiasis. Currently, 37 species of Bulinus represent four “groups”: the “Bu. truncatus/tropicus” complex, “Bu. reticulatus”, “Bu. forskalii” and “Bu. africanus” (Brown, 1994; Rollinson et al., 2001). Traditionally, members of these groups (or operational taxonomic units, OTUs) have been identified using morphological characters (Brown, 1994; Southgate et al., 2000), but specific identification can be unreliable due to a limited number of informative phenotypic characters (e.g. shell morphology) (Stothard et al., 1997). Nonetheless, karyotyping, multilocus enzyme electrophoresis (MEE) and DNA-based methods, particularly those based on PCR, have assisted in discerning species or taxa (Rollinson & Southgate, 1979; Raahauge & Kristensen, 2000; Zein-Eddine et al., 2014). For both fundamental and applied studies, there is a clear need for molecular tools for the accurate assignment of snails to the Bulinus groups, but, more importantly, for the unequivocal identification and differentiation of species within each of the four recognised groups.

Considerable work (Stothard et al., 2002) has been done in this direction, with a particular focus on defining, assessing and/or using genetic markers in nuclear ribosomal and in mitochondrial (mt) DNAs for specific identification, particularly the second internal transcribed spacer (ITS-2) and cytochrome c oxidase subunit 1 gene (cox1), respectively (Stothard et al., 2002; Jorgensen et al., 2013). Indeed, mt cox1 sequences have shown major utility for the identification of species and groups, and to establish the genetic/systematic relationships of species/taxa within the genus Bulinus (Tumwebaze et al., 2019; Chibwana et al., 2020), although, in some instances, there has been limited statistical support for some relationships (Kane et al., 2008). The latter finding indicates that the use of a greater number of phylogenetically-informative genetic markers should enhance the taxonomic classification of Bulinus species or subspecies, and could improve the phylogenetic reconstruction of systematic relationships. A logical extension to this published work is to sequence and annotate a complete mt genome to represent the genus Bulinus, in order to underpin future taxonomic/phylogenetic studies of members of the genus.

Given our recent success with the use of third-generation (long-read or long-range) sequencing technologies to define the complete mt genomes of a number of invertebrate species (parasitic flatworms), including Echinococcus granulosus (genotype G1) and Clonorchis sinensis (Kinkar et al., 2019, 2020), we here used this technology, combined with second-generation (short-read) sequencing and advanced informatics, to define the mt genome of a first, key representative of the genus – Bu. truncatus – to accelerate future mitogenomic and genetic explorations of the snails that act as intermediate hosts of members of the S. haematobium group.

Materials and methods

Isolation and procurement of the snails

Samples of Bu. truncatus originated from a laboratory strain (designated ‘BRI’) which is maintained in the Biomedical Research Institute (BRI), Rockville, Maryland (Lewis et al., 2008); this strain was originally sourced from Egypt (Dr Margaret Mentink-Kane, personal communication, 10 October 2020). Individual snails were washed extensively in phosphate-buffered saline (PBS, pH 7.0) and frozen at −80 °C.

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 an Agilent 4200 TapeStation system (ThermoFisher) and using Genomic DNA ScreenTape (ThermoFisher). Low molecular weight DNA was removed using a 10 kb short-read eliminator kit (Circulomics, Baltimore, MD, USA). The high molecular weight DNA was used to construct Nanopore Rapid Sequencing (SQK-RAD004; Oxford Nanopore Technologies) and Ligation Sequencing (SQK-LSK109; Oxford Nanopore Technologies) genomic DNA libraries, according to the manufacturerʼs instructions. Each flow cell used to sequence the first libraries was washed using a Flow Cell Wash Kit (EXP-WSH003; Oxford Nanopore Technologies) and re-used to sequence additional SQK-LSK109 libraries. All libraries were sequenced using the MinION sequencer (Oxford Nanopore Technologies). Following sequencing, bases were ‘called’ from raw FAST5 reads using the program Guppy v.3.1.5 (Oxford Nanopore Technologies) and stored in the FASTQ format (Cock et al., 2010). A short-insert (500 bp) genomic DNA library was also constructed and paired-end sequenced (150 bp reads) using TruSeq sequencing chemistry and the NovaSeq sequencing platform (Illumina).

Assembly and annotation of the mt genome

Long reads were assembled using FLYE v.2.6 (Kolmogorov et al., 2019) with the --nano-raw option and setting a genome size estimate of 900 megabases (Gregory, 2003). Errors in long-read sequence data were initially corrected using medaka_consensus in the Medaka package v.0.10.0 (https://github.com/nanoporetech/medaka). The assembled genome was then polished with Pilon v.1.23 (Walker et al., 2014) using the short-read Illumina data. Finally, long- and short-reads were mapped to the assembled mt genome using Minimap2 v.2.0 (Li, 2018), and coverage of the genome was determined using mpileup in the SAMtools package v.1.9 (Li et al., 2009).

Annotation of tRNA, rRNA and protein-encoding gene regions was initially undertaken using MITOS webserver (Bernt et al., 2013) and employing the mt genetic code for invertebrates (https://www.ncbi.nlm.nih.gov/Taxonomy/Utils/wprintgc.cgi; translation_table_5). The open reading frame (ORF) of each protein gene was then further curated in the program Geneious v.11.1.5 (Kearse et al., 2012) and using gene regions of published mt genomes of planorbid snails as a guide (Table 1). The complete mt genome sequence was deposited in the GenBank database under the accession no. MT947902; raw data are also available from the Sequence Read Archive (SRA) under the accession no. SAMN17050146 with the NCBI BioProject accession number PRJNA680620.

Table 1

Mitochondrial genomic sequences for snail species/strains used in the present study, with accession numbers and references listed

Accession number	Species – “strain”	Length (bp)	G ​+ ​C content (%)	Reference
MT947902	Bulinus truncatus	13,767	24.3	Present study
NC_026538a	Radix auricularia – ZUE strain	13,745	28.6	Feldmeyer et al. (2015)
MT628577.1	Gyraulus laevis	13,685	25.2	Unpublished
KY514384.1	Planorbella duryi	14,217	27.3	Schultz et al. (2018)
MT628573.1	Bathyomphalus contortus	13,679	24.9	Unpublished
EF433576	Biomphalaria tenagophila – Taim-RS strain	13,722	24.2	Unpublished
MF480756	Biomphalaria straminea voucher Xiaoping Wu:04-BS-2017	13,650	24.7	Zhou et al. (2019)
MG431962	Biomphalaria pfeifferi	13,624	23.4	Zhang et al. (2018)
MG431964	Biomphalaria choanomphala	13,672	23.4	Zhang et al. (2018)
MG431963	Biomphalaria sudanica	13,671	23.4	Zhang et al. (2018)
MG431966	Biomphalaria glabrata – G16 BS90 strain	13,676	24.8	Zhang et al. (2018)
MG431965	Bi. glabrata – G72 M strain	13,667	25.3	Zhang et al. (2018)
AY380531	Bi. glabrata – 1742 strain	13,670	25.4	DeJong et al. (2004)
AY380567	Bi. glabrata – M strain	13,670	25.4	DeJong et al. (2004)

Member of the family Lymnaeidae – used as an outgroup in phylogenetic analyses.

Single loci and whole genome comparative analyses

Available cox1 sequence data for Bulinus spp., selected taxa of the family Planorbidae and Radix auricularia (family Lymnaeidae; outgroup) were downloaded from NCBI (Supplementary Table S1) and aligned as nucleotide sequence using MUSCLE v.3.7 (Edgar, 2004). Sequence alignments were trimmed to the coding region of cox1 (164–1,098 of the 1,542 bp Bu. truncatus cox1 gene) in order to minimise missing data from the 5' and 3' regions of previously sequenced amplicons (Supplementary Table S1). Identical nucleotide sequences were removed using CD-HIT-EST v.4.6 (Fu et al., 2012). The optimal nucleotide substitution model for aligned sequences was then assessed using the program ModelTest-NG v.0.1.6 (Darriba et al., 2020). The aligned sequences were then subjected to phylogenetic analysis using Bayesian inference (BI) employing Monte Carlo Markov Chain analysis in the program MrBayes v.3.2.2 (Ronquist et al., 2012). Posterior probabilities (pp) were calculated using the GTR ​+ ​I ​+ ​G model, generating 1,000,000 trees and sampling every 200th tree until potential scale reduction factors for each parameter approached one. The initial 25% of trees were discarded as “burn-in”, and the others were used to construct a majority rule tree. Phylogenetic trees were rendered and annotated using ggtree v.1.10.5 (Yu et al., 2017) in R v.3.4.3 (http://www.R-project.org/).

The complete mt genome of Bu. truncatus was compared to the available reference mt genomes of other planorbid snails (Table 1) and R. auricularia (family Lymnaeidae; outgroup) using progressiveMauve v.2.4.0 (Darling et al., 2004), and using the following settings: --hmm-identity =0.95 --island-gap-size = 10. Mitochondrial protein-coding genes were then extracted and aligned as individual nucleotide coding or inferred amino acid sequences using MUSCLE. The optimal substitution model for each aligned sequence was then assessed using the program ModelTest-NG v.0.1.6 (Darriba et al., 2020). The aligned sequences were then subjected to phylogenetic analysis using BI or maximum likelihood (ML) methods employing Monte Carlo Markov Chain analysis in the program MrBayes and RaxML v.8.2.12 (Stamatakis, 2014). For BI, posterior probabilities (pp) were calculated using the selected substitution model (Nucleotide sequence: GTR ​+ ​I ​+ ​G; Protein sequence: MtMam for ATP6, COX1-2, CYTB, NAD1-2 and MtRev for COX3 and NAD3-6), generating 1,000,000 trees and sampling every 200th tree until potential scale reduction factors for each parameter approached one. The initial 25% of trees were discarded as “burn-in”, and the others were used to construct a majority rule tree. For ML the same sequence set was used to construct a partitioned ML tree using the selected models for each partition and performing 20 iterations. Nodal support values for the ML tree were inferred from 100 bootstrap replicates. Phylogenetic trees were rendered and annotated using ggtree in R.

Patterns of nucleotide diversity were compared between the aligned mt protein-coding regions of Bu. truncatus and protein-coding regions of available reference mt genomes of planorbid snails and R. auricularia (family Lymnaeidae) as an outgroup (Table 1). A sliding window analysis of nucleotide diversity (steps of 10 bp over 200-bp windows) was performed for each pairwise-alignment of concatenated genes using the PopGenome package (Pfeifer et al., 2014) in R. For each comparison, nucleotide diversity values were plotted using the R package ggplot2 (Wickham, 2009).

Results and discussion

Characteristics of the mitochondrial genome of Bulinus truncatus

From a total of 1,327 long- and 917,758 short-reads (equating to ∼139 Mb), a complete, 13,767 bp circular mt genome was assembled (GenBank: MT947902). The mean coverage of long-reads mapped to the genome was 136 (range: 53–199; Fig. 1). The mean coverage of mapped short-reads was 8,902 (range: 3,892–9,730), with the depth being ∼4,000 reads at the introduced start/stop positions (Fig. 1). We identified 37 genes, including 13 protein-coding (adenosine triphosphatase subunit 6 (atp6), adenosine triphosphatase subunit 8 (atp8), the cytochrome c oxidase subunits 1, 2 and 3 (cox1–cox3), cytochrome b (cytb) and the nicotinamide dehydrogenase subunit 1–6 (nad1–nad6 and nad4l), 22 transfer RNA (two coding for leucine and two coding for serine)) and the small (rrnS) and large (rrnL) subunits of ribosomal RNA genes (Fig. 1 and Table 2). In most instances, start and stop codons were consistent with those of the mt genomes of other invertebrates characterised to date (Lavrov, 2014).

Fig. 1

Reference mitochondrial genome of Bulinus truncatus. Direction of gene transcription is shown with an arrow. Long (rrnL) and short (rrnS) ribosomal RNA subunits are shown in red and protein-encoding genes are shown in black. The graph shows the depth of individual nucleotides of individual long- (LR) and short-read (SR) reads at individual positions in the mt genome (grey dots); indicated is also the mean nucleotide depth across this genome (solid, dark-grey line)

Table 2

Mitochondrial genes of Bulinus truncatus and their locations, GC contents, lengths, direction of transcription and start/stop codons

Gene designations	Location (start/end)	GC (%)	Length (bp)	Transcription direction	Start/stop codons
Protein-coding genes
atp6	5284/5940	22.4	789	Reverse	AAT/TAAa
atp8	5128/5235	18.3	108	Reverse	ATT/TAA
cox1	10517/12058	29.6	1542	Forward	ATT/TAA
cox2	4114/4785	27.7	729	Forward	ATT/TAG
cox3	8666/9463	28.8	882	Reverse	ATG/TAG
cytb	2828/3924	27.3	1107	Forward	TTGa/TTAa
nad1	1647/2540	26.2	894	Forward	ATG/TAA
nad2	9570/10514	23.0	1044	Forward	ATG/TAA
nad3	6819/7184	18.6	378	Reverse	ATA/TAA
nad4	7308/8612	22.3	1404	Forward	ATA/TAA
nad4l	2566/2862	21.2	327	Forward	ATA/TAG
nad5	14/1663	22.8	1650	Forward	ATA/TAG
nad6	13312/13767	18.4	456	Forward	TTG/TAA
Ribosomal RNA genes
rrnS	6064/6776	23.6	713	Reverse	na
rrnL	12104/13153	21.4	1093	Forward	na
Transfer RNA genes
tRNA-D(gtc)	3925/3991		67	Forward	na
tRNA-C(gca)	3987/4048		62	Forward	na
tRNA-F(gaa)	4051/4112		62	Forward	na
tRNA-Y(gta)	4768/4826		59	Forward	na
tRNA-W(tca)	4823/4885		63	Forward	na
tRNA-G(tcc)	4885/4950		66	Forward	na
tRNA-H(gtg)	4942/5003		62	Forward	na
tRNA-Q(ttg)	5006/5066		61	Reverse	na
tRNA-L2(taa)	5056/5125		70	Reverse	na
tRNA-N(gtt)	5236/5303		68	Reverse	na
tRNA-R(tcg)	5944/6005		62	Reverse	na
tRNA-E(ttc)	6001/6064		64	Reverse	na
tRNA-M(cat)	6775/6838		64	Reverse	na
tRNA-S2(tga)	7185/7240		56	Reverse	na
tRNA-S1(gct)	7242/7299		58	Forward	na
tRNA-T(tgt)	8621/8685		65	Reverse	na
tRNA-I(gat)	9504/9570		67	Forward	na
tRNA-K(ttt)	10474/10537		64	Forward	na
tRNA-V(tac)	12062/12125		64	Forward	na
tRNA-L1(tag)	13130/13194		65	Forward	na
tRNA-A(tgc)	13191/13254		64	Forward	na
tRNA-P(tgg)	13255/13315		61	Forward	na

Abbreviation: na, not applicable.

The locations of start and stop codons were identified based on the alignment of protein-coding regions, inferred using published gene models for planorbid snails (Feldmeyer et al., 2015; Schultz et al., 2018; Zhang et al., 2018).

The mt genomic sequence of Bu. truncatus is biased toward A ​+ ​T (75.7%), with T (41.7%) being the most frequent nucleotide, and C (10.9%) the least frequent one, in accord with mt genomes of some other planorbid snails including Bathyomphalus contortus, Biomphalaria spp., Gyraulus laevis and Planorbella duryi (Table 1). A ​+ ​T content was highest in the genes atp8, nad3 and nad6 (81.4–81.7%) and lowest in the genes cox1, cox2 and cox3 (70.4–72.3%) (Table 2). The tRNAs (Table 2) were inferred to have a canonical structure (Supplementary Figure S1). Two copies of a serine and a leucine tRNA were encoded, and all tRNAs were predicted to have DHU and TψC arms, except for tRNA-G(tcc) (without a TψC arm), tRNA-S1(gct) (without a DHU arm) and tRNA-S2(tga) (without a DHU arm).

Verifying the taxonomic status of Bulinus truncatus

As the molecular taxonomy of snails relies on the use of cox1 as the mt genetic marker (Stothard et al., 2002; Pennance et al., 2020), we used extensive, publicly-available (partial) cox1 sequence data sets in a phylogenetic analysis to verify the position of the BRI strain of Bu. truncatus in relation to previously characterised Bulinus groups and taxa (Fig. 2). There was strong nodal support for the present strain grouping with a large number of other Bu. truncatus specimens of the Bu. truncatus/tropicus complex, to the exclusion of members of the three other groups (Bu. reticulatus, Bu. forskalii and Bu. africanus), although there was limited support for grouping with other Bu. truncatus specimens from Egypt (Fig. 2). The BRI strain of Bu. truncatus clustered most closely with isolates KM272998, MG407308, MG407310, MG407312, MG407326 and MG759391, all of which originated from Egypt (Supplementary Table S1; Fig. 2B). The phylogenetic tree constructed using partial cox1 sequence data also revealed limited nodal support for some nodes within some recognised Bulinus groups.

Fig. 2

A Phylogenetic relationships of Bulinus truncatus with other representatives of this species and other (selected) species of planorbid snails (Supplementary Table S1) inferred based on an analysis of the aligned partial cox1 gene sequences by Bayesian inference using Radix auricularia (family Lymnaeidae) as an outgroup. The outermost black line indicates recognised groups of Bu. truncatus (see Brown, 1994). B Part of the phylogenetic tree includes a group comprising Bu. truncatus from the Biomedical Research Institute [BRI] (Lewis et al., 2008) – sequenced herein – and other representatives from Egypt (EGY), Sudan (SUD) and Zimbabwe (ZWE). Posterior probability (pp) support is given at each node of the tree and relates to a heatmap (red: pp = 1; black: pp = 0.5). The partial cox1 gene sequence of the BRI strain of Bu. truncatus sequenced here is denoted in bold-type

Mt genomic comparisons

Having established the characteristics of the mt genome of Bu. truncatus and verified the taxonomic status of our BRI strain, we established levels of nucleotide variability along the whole mt genome between Bu. truncatus, other planorbid snails (Ba. contortus, Biomphalaria spp., G. laevis and P. duryi) using R. auricularia as a distant taxon (outgroup) (Fig. 3). The comparison of the aligned genome sequences revealed 38 conserved, co-linear blocks with ≥ 95% nucleotide identity (comprising 8,582 of the 13,506 aligned nucleotide positions). There was no observed difference in gene order or direction of transcription between Bu. truncatus and other planorbid snails studied here; cox1 (1,563 bp), cox3 (1,106 bp) and cytb (860 bp) were the longest protein-encoding genes (Table 2). A guide-tree, inferred from a whole genome alignment (Fig. 3, left), indicated that Bu. truncatus grouped with G. laevis to the exclusion of Biomphalaria species and other planorbids (Ba. contortus and P. duryi). An assessment of nucleotide diversity in protein-coding genes between Bu. truncatus and other snail species upon pairwise comparison (Fig. 4) revealed low (0.153–0.244) and high (0.215–0.391) mean diversities in the cox and nad genes, respectively. Across all protein-coding genes of Bu. truncatus, most diversity was seen with the lymnaeid R. auricularia (mean nucleotide diversity = 0.304; Fig. 4).

Fig. 3

Alignment of the Bulinus truncatus mt genome with other (selected) species of planorbid snails and an outgroup – Radix auricularia (family Lymnaeidae) (cf. Table 1). Co-linear blocks of nucleotides ≥ 95% identity (maximum of 10 nucleotide gaps allowed within block) are shown. The guide-tree represents the genetic distance between genomes based on shared gene content among aligned genome pairs. The scale-bar indicates estimated evolutionary distance

Fig. 4

Nucleotide diversity in coding genes between Bulinus truncatus and other available planorbid snails as well as Radix auricularia (outgroup) upon pairwise comparison (cf. Table 1). Each solid line represents the average nucleotide diversity across 200, 10-base overlapping nucleotide windows. All species of Biomphalaria are represented in one plot, with each species/strain represented by a distinct grey line. Average nucleotide diversity across all coding regions is indicated by a horizontal, dashed line. Average nucleotide diversities for each gene region are indicated. Asterisk indicates nucleotide diversity values affected by missing data for some species

Relationship of Bulinus truncatus with other planorbids

The phylogenetic/systematic relationship of Bu. truncatus with other planorbid species was constructed using complete data sets publicly available for mt protein-coding genes (Table 1; Fig. 5). Trees constructed using nucleotide or amino acid sequence data provided strong support (posterior probability = 1.00; bootstrap support = 100%) for most nodes, with Bu. truncatus being basal. Using nucleotide sequence data, support was limited (pp = 0.82; bss = 70%) only at nodes grouping Bi. glabrata 1742 with Bi. glabrata M-line as well as Bi. straminea with some other species, such as Bi. tenagophila (pp = 0.89; bss < 50%). Using amino acid sequence data, Bu. truncatus grouped with Ba. contortus and G. laevis (pp = 0.99; bss = 73%), with P. duryi being basal. The basal position of Bu. truncatus in the tree (Fig. 5) agrees with a proposed ancestral position of Bulinus and Indoplanorbis (Albrecht et al., 2007) and supports the previous proposals that Bulinus belong to the family Bulinidae (Albrecht et al., 2007; Bouchet et al., 2017). The characterisation of the mt genomes of other species of Bulinus and Indoplanorbis should help resolve the taxonomic ambiguities within the superfamily Lymnaeoidea.

Fig. 5

Phylogenetic relationship of Bulinus truncatus with selected species of planorbid snails and Radix auricularia (family Lymnaeidae) as an outgroup (cf. Table 1). A Phylogeny inferred from concatenated nucleotide sequences derived from 12 mitochondrial protein-encoding genes using Bayesian inference (BI) or maximum likelihood (ML). B Phylogeny inferred from concatenated amino acid sequences derived from 12 mitochondrial protein-encoding genes using BI or ML. Posterior probability (pp; BI) or percentage bootstrap (bs; ML) support values are indicated at each node of the tree (BI/ML). The scale-bar indicates phylogenetic distance (in substitutions per site)

Implications and future work

Here, we characterised the first complete mt genome representing the genus Bulinus – a complex of at least 37 species presently divided into four main ‘groups’ (Brown, 1994). This mt genome for the BRI strain of Bu. truncatus was assembled using a combination of second- (short-read) and third-generation (long-read) sequence data, enabling taxonomic verification (using cox1 sequence data) as well as comparative analyses with other planorbid snails for which mt genomes were known.

The taxonomic investigation using cox1 sequence data showed that the BRI strain of Bu. truncatus clustered with six samples from Egypt, including those from locations in Tanta (KM272998), Mansouria (MG407308, MG407310 and MG407312), Behera (MG407312) and Giza (MG759391) (Zein-Eddine et al., 2014; Abe et al., 2018). While an exact phylogeographic location was not determined, the clustering of the BRI strain with these Bu. truncatus isolates is consistent with information that this laboratory strain was originally collected from Egypt (Dr Margaret Mentink-Kane, BRI, personal communication, 10 October 2020). The partial cox1 sequence data sets publicly available were not sufficiently phylogenetically informative to yield statistically well-supported trees/phylogeographic relationships (Fig. 2) – a finding that is consistent with evidence provided in previous studies (Zein-Eddine et al., 2014; Abe et al., 2018; Tumwebaze et al., 2019) – although identification and differentiation should be achieved if cox1 sequences of at least 1 kb are employed (Kane et al., 2008). The molecular characterisation of Bulinus species has usually relied on the use of single mt and/or nuclear (e.g. ITS-2) genetic markers (cf. Kane et al., 2008; Chibwana et al., 2020; Pennance et al., 2020). In most cases, identification or differentiation of members of the four major Bulinus groups was achieved (Kane et al., 2008; Zein-Eddine et al., 2014; Pennance et al., 2020), but, in some instances, incongruences were observed due to limited signal (saturation of mutational variation) or nodal support (Stothard et al., 2001; Jorgensen et al., 2013; Zein-Eddine et al., 2014; Abe et al., 2018; Tumwebaze et al., 2019).

Defining the first complete reference mt genome for Bu. truncatus (BRI strain) using short- and long-read data sets paves the way to characterising such genomes for the 37 known members of the genus Bulinus and then, using complete protein-coding gene data sets, to test the validity of the four currently-recognised/-proposed taxonomic groups (“Bu. truncatus/tropicus” complex, “Bu. reticulatus”, “Bu. forskalii” and “Bu. africanus”). In the future, it will be important to conduct mt genome studies of type (museum) specimens of Bulinus species using refined methods for the isolation and sequencing of ancient DNA (e.g. Spyrou et al., 2019). Using all of these mt genomic data sets, it should be possible to define genetic markers for distinct OTUs as a basis to conduct molecular epidemiological investigations of the prevalence and geographical distribution of particular Bulinus taxa involved in the transmission of urogenital schistosomiasis, explorations of OTU-schistosome affiliations as well as large-scale population genetic studies of the snail vectors of members of the S. haematobium-group (cf. Rollinson et al., 2001). These efforts would assist in providing an informed position regarding the control of urogenital schistosomiasis and the feasibility of eliminating urogenital schistosomiasis in endemic regions in the context of mass drug administration regimens (Knopp et al, 2012, 2019; Allan et al., 2020).

The relatively conserved gene order in the mt genomes and the phylogenetic relationships of the snails studied here (Fig. 3, Fig. 5) show the feasibility of rapidly characterising the mt genomes of a range of members of the family Planorbidae – some of which are vectors of socioeconomically important parasitic trematodes other than schistosomes – using the technological approach established here. Such an effort could assist significantly in closing some of the knowledge gaps that exist in the understanding of systematics of these groups, and would allow some insights into molecular evolution and codon usage (cf. Clary & Wolstenholme, 1985; Batuecas et al., 1988; Hoffmann et al., 1992; Boore & Brown, 1994). Complementing such an effort, it would be imperative to start to sequence the nuclear genomes of respective snail vectors of key parasitic trematodes.

Given the prime importance of Bi. glabrata as a vector of Schistosoma mansoni, draft nuclear genomes have been published for this species (Adema et al., 2017; Lu et al., 2020), although there have been challenges associated with karyotyping (Goldman et al., 1984) and ploidy determination (Wheeler et al., 2018) for this snail. Nonetheless, through major advances in sequencing and informatics, a recent study (Kenny et al., 2020) has succeeded in yielding a chromosome-level assembly of a large and complex genome of the mollusc Pecten maximus – great scallop (family Pectinidae) – which bodes well for the application of a similar workflow to sequence the genomes of molluscan vectors of trematodes. As a key intermediate host for S. haematobium, Bu. truncatus would be an ideal candidate for future nuclear genome characterisation and comparative genomic explorations. We anticipate that the ploidy/karyotype(s) of Bulinus (Wu, 1972) will present unique challenges for genome assembly efforts, but, if successful, the availability of both mt and nuclear genomes for this snail species would represent major resources for studies of its ecology, population genetics, origin and evolution as well as the epidemiology of urogenital schistosomiasis and snail host-parasite affiliation(s).

In conclusion, the present study emphasises the relevance of genomic work on Bu. truncatus and other key snail vectors of parasitic trematodes, as a basis for future systematic, population genetic, epidemiological, ecological and biological investigations. The sequencing and bioinformatic platform established here provides exciting prospects for large-scale comparative and population genetic studies of such snail vectors. Fundamental insights into snail genetics and biology could assist in better understanding the transmission patterns of trematodiases and in finding ways of blocking disease transmission to humans and other vertebrates.

Funding

This research project was supported by grants from the (NHMRC) and the (ARC) to RBG, NDY and PKK.

CRediT author statement

Neil D. Young: Methodology, Writing - Original Draft, Conceptualisation, Writing - Review & Editing, Supervision, Project administration, Funding acquisition. Liina Kinkar: Methodology, Writing - Review & Editing. Andreas J. Stroehlein: Methodology, Writing - Review & Editing. Pasi K. Korhonen: Writing - Review & Editing, Funding acquisition. David Rollinson: Writing - Review & Editing. J. Russell Stothard: Writing - Review & Editing. Robin B. Gasser: Conceptualisation, Writing - Original Draft, Writing - Review & Editing, Supervision, Project administration, Funding acquisition. All authors read and approved the final manuscript.

Data availability

The nucleotide sequence of the mitochondrial (mt) genome reported in this article is publicly available in the GenBank database under accession no. MT947902. Raw data are also available from the Sequence Read Archive (SRA) under the accession no. SAMN17050146 with the NCBI BioProject accession number PRJNA680620. All other data analysed in this study are referred to in this article and its supplementary files.

Declaration of competing interests

The authors declare that they have no competing interests.