Large-scale and small-scale population genetic structure of the medically important gastropod species Bulinus truncatus (Gastropoda, Heterobranchia)
Parasites & Vectors volume 15, Article number: 328 (2022)
Gastropod snails remain strongly understudied, despite their important role in transmitting parasitic diseases. Knowledge of their distribution and population dynamics increases our understanding of the processes driving disease transmission. We report the first study to use high-throughput sequencing (HTS) to elucidate the population genetic structure of the hermaphroditic snail Bulinus truncatus (Gastropoda, Heterobranchia) on a regional (17–150 km) and inter-regional (1000–5400 km) scale. This snail species acts as an intermediate host of Schistosoma haematobium and Schistosoma bovis, which cause human and animal schistosomiasis respectively.
Bulinus truncatus snails were collected in Senegal, Cameroon, Egypt and France and identified through DNA barcoding. A single-end genotyping-by-sequencing (GBS) library, comprising 87 snail specimens from the respective countries, was built and sequenced on an Illumina HiSeq 2000 platform. Reads were mapped against S. bovis and S. haematobium reference genomes to identify schistosome infections, and single nucleotide polymorphisms (SNPs) were scored using the Stacks pipeline. These SNPs were used to estimate genetic diversity, assess population structure and construct phylogenetic trees of B. truncatus.
A total of 10,750 SNPs were scored and used in downstream analyses. The phylogenetic analysis identified five clades, each consisting of snails from a single country but with two distinct clades within Senegal. Genetic diversity was low in all populations, reflecting high selfing rates, but varied between locations due to habitat variability. Significant genetic differentiation and isolation by distance patterns were observed at both spatial scales, indicating that gene flow is not strong enough to counteract the effects of population bottlenecks, high selfing rates and genetic drift. Remarkably, the population genetic differentiation on a regional scale (i.e. within Senegal) was as large as that between populations on an inter-regional scale. The blind GBS technique was able to pick up parasite DNA in snail tissue, demonstrating the potential of HTS techniques to further elucidate the role of snail species in parasite transmission.
HTS techniques offer a valuable toolbox to further investigate the population genetic patterns of intermediate schistosome host snails and the role of snail species in parasite transmission.
In 2014, the WHO highlighted vector-borne diseases as a global public health priority. These diseases are either actively transmitted through vectors, such as mosquitos, sand flies or ticks, or passively by snail intermediate hosts. Despite their medical importance, snail intermediate hosts remain largely understudied . Bulinus truncatus (Audouin, 1827) is one of the main intermediate hosts of Schistosoma haematobium (Bilharz, 1852), the causative agent of human urinary schistosomiasis. It is also an intermediate host of the bovine parasite Schistosoma bovis (Bilharz, 1852) and different Paramphistoma and Echinostoma species (Brown, 1994). Bulinus truncatus is a tetraploid species belonging to the monophyletic B. truncatus/tropicus complex [2,3,4,5]. It is a typical r-species characterised by a high reproduction potential, the ability to self-fertilise in new environments, wide tolerance ranges and good dispersal capacities, facilitated by passive spreading through river currents or on the feet and feathers of waterbirds [5,6,7,8,9].
The distribution and abundance of B. truncatus, and of the associated schistosomiasis risk, are heavily impacted by anthropogenic activities. Dam construction has been shown to alter snail abundances differently in various areas [10,11,12,13], and human-induced climate change is predicted to shift snail host distributions poleward while reducing habitat suitability in some parts of sub-Saharan Africa [14, 15]. Furthermore, increasing human mobility in combination with rising temperatures could (re)introduce schistosome parasites in areas where a suitable intermediate host is present, as exemplified by the recent foci in Corsica [15,16,17,18,19]. These anthropogenic influences will likely lead to a schistosomiasis risk expansion to more temperate regions and a risk reduction in some tropical regions .
Little is known about the population genetic structure of B. truncatus and the role of snail genetics in determining parasite susceptibility [21,22,23,24,25]. Furthermore, a reference genome has only recently become available for Biomphalaria glabrata (Say, 1818)  but remains absent for other intermediate host species. These knowledge gaps need to be addressed in order to fully understand schistosome transmission dynamics.
This study aims to provide a first characterisation of the genetic structuring of B. truncatus populations on a regional and inter-regional scale, ranging from 17 to 5400 km, using genotyping-by-sequencing (GBS) . To the best of our knowledge, this is the first study where high-throughput sequencing (HTS) has been used to assess the population genetic structure of schistosome intermediate snail hosts on multiple spatial scales.
Samples of B. truncatus (range in numbers of specimens collected per locality: 4–16) were collected from 2011 to 2014 in Cameroon (1 site: Cam_BAK), Senegal (4 sites: Sen_DIAM, Sen_GUEC, Sen_MBO, Sen_NDO), Egypt (2 sites: Egy_BEK, Egy_MAN) and France (1 site: Cor_SUT) (Table 1; Fig. 1). The morphological identification of specimens was verified using DNA barcoding . DNA was extracted using the DNeasy Blood and Tissue Kit (Qiagen, Hilden, Germany) following the manufacturer's protocol, and a 658-bp fragment of cytochrome c oxidase subunit I (COI) was amplified using primers LCO1480 and HCO2198  and/or primers BulCox6 and BulCox12 . PCR assays were performed in a 25-μl reaction volume containing 1.5 mM MgCl2, 0.2 mM dNTPs, 0.4 µM of each primer, 0.75 U of Platinum Taq DNA Polymerase (Invitrogen™, Thermo Fisher Scientific, Waltham, MA, USA), 2.5 µl 10× reaction buffer (Invitrogen™) and 1.5 µl of DNA extract (with DNA concentrations ranging from 10 to 100 ng/l). The PCR cycling parameters consisted of an initial denaturation at 94 °C for 3 min, followed by 35 cycles of denaturation at 94 °C for 60 s, annealing at 55 °C for 60 s and extension at 72 °C for 90 s, with a final extension at 72 °C for 7 min. PCR products were purified by passage through GFX purification columns (GE Healthcare, Chicago, IL, USA), subjected to sequencing reactions using the 3.1 BigDye Cycle Sequencing Kit (Applied Biosystems, Thermo Fisher Scientific, Foster City, CA, USA) and sequenced in both directions with an ABI Prism 3100 Genetic Analyzer (Applied Biosystems). Nucleotide sequences were aligned using the default parameters of the IUB scoring matrix of ClustalW, as implemented in Geneious R9 . Primers were trimmed, and the fragment was translated into amino acids to check for the absence of internal stop codons. The R package Adhoc  was used to implement molecular identifications using all reference sequences publicly available for the genus Bulinus in the Barcode Of Life Data Systems (BOLD, http://www.boldsystems.org/). Identification was based on the Best Close Match criterion with 1% distance threshold (corresponding to K2P ≤ 0.01) . Genetic distances between the BOLD reference sequences and the 87 B. truncatus queries from this study (plus 8 other queries from specimens that were eventually discarded due to their low genotyping-by-sequencing [GBS] output) were visualised through a Neighbour Joining (NJ) tree (see Additional file 3: Figure S1) based on K2P genetic distance  as calculated in Mega 6 .
A single-end GBS reduced representation library (RRL) was built and sequenced on an Illumina HiSeq 2000 platform (Illumina, Inc., San Diego, CA, USA) at the Cornell University Biotechnology Resource Center, Ithaca, NY, USA. The library was prepared by digesting DNA with EcoT22I, ligating P1 and P2 Illumina sequencing primers (the former including molecular identifiers for sequence demultiplexing) and Illumina adapters, before proceeding with titration and library enrichment as described in Elshire et al. . The choice of the restriction enzyme followed preliminary comparisons between EcoT22I and PstI on representatives of B. truncatus.
Processing of HTS data
The HiSeq run produced 257 × 106 single-end raw reads of 101 bp. FastQC (http://www.bioinformatics.babraham.ac.uk/projects/fastqc/) was used for sequence quality check before trimming the last 7 bp of raw sequences, as implemented in FASTX-Toolkit (http://hannonlab.cshl.edu/fastx_toolkit/). We explored if, and to what extent, parasite infection could be detected using GBS data by adding two snail samples from an additional site of Senegal (Sen_PAK) that were known to be infected with S. bovis (i.e. cercariae collected after shedding in Boon et al. ). To remove possible contaminant sequences originating from Schistosoma spp. flukes, reads were mapped against the S. bovis and S. haematobium reference genomes (GenBank assembly GCA_003958945.1 and GCA_000699445.1, respectively) using DeconSeq 0.4.3 applying default settings. Decontaminated data were processed using the Stacks 1.21 pipeline  and implementing: (i) process_radtags for demultiplexing and quality filtering and (ii) denovo_map.pl to call single nucleotide polymorphisms (SNPs). After some pilot runs with various parameter settings, denovo_map.pl considered a minimum threshold of 10 raw reads per stack (m), maximum two mismatches between loci from a single individual (M), maximum two mismatches when aligning secondary reads to primary stacks (N) and maximum two mismatches between loci when building the catalogue (n). SNPs were scored when occurring in at least four samples. VCFtools v0.1.14  was used to filter out specimens with > 60% missing data, and SNPs with minimum allele frequencies < 0.05. To eliminate putative paralogues we conducted a heterozygosity excess test for each SNP using VCFtools and applied a threshold level of significance of 0.01. The data were further filtered by selecting one single SNP per restriction site associated DNA marker (RAD tag). This generated a main dataset of 10,750 polymorphic sites scored in 87 specimens (following discarding of 8 specimens with low-quality GBS profiles from an initial set of 95 vouchers; see Additional file 1: Table S1). When not specified otherwise (see below), this dataset was used for downstream analyses.
Population structure and diversity
Allelic richness (ar) and observed and expected heterozygosity (Hobs, Hexp) were calculated using the R package ‘diveRsity’ . Arlequin 3.5  was used to test significant departures from Hardy–Weinberg equilibrium (HWE; following  and to calculate and test inbreeding coefficients (Fis) , as well as pairwise fixation index (Fst) values  per population (10,000 permutations). Probability values of repeated tests were corrected for Type I errors using the false discovery rate procedure  with P = 0.001. The function glPCA of the R package ‘adegenet’  was used to perform exploratory principal component analyses (PCAs) of raw data. Levels of genetic admixture among samples were quantified via STRUCTURE 2.3.4  after inferring the optimal number of clusters (K) following the method of Evanno et al.  in Structure Harvester 0.6.94 . For each value of K, 10 iterations were run for 200,000 generations (admixture model, burn-in of 100,000 generations), and the posterior estimates of cluster memberships of the three runs with the highest estimated log probability of the data were summarised using CLUMPP 1.1.2  and visualised in DISTRUCT 1.1 .
Isolation by distance (IBD) was verified by a Mantel test between Euclidean genetic and Euclidean geographic distances. The occurrence of discontinuities in the distribution of inter-individual distances was visualised through two-dimensional kernel density as implemented by the function kde2d of “adegenet” . Geographic patterns of genetic structuring were further investigated via the multivariate spatial autocorrelation method of Smouse and Peakalll [51, 52] implemented in GENALEX 6.5 [53, 54]. We randomly selected a subset of 917 SNPs by setting the ‘–thin’ option of VCFtools to 15,000, and pairwise geographic distances were calculated as kilometric distance between x- and y-UTM individual coordinates. The spatial autocorrelation coefficient r assesses the pairwise genetic similarity between individuals sampled within a specified distance class. We defined multiple distance class sizes (‘MultipleDclass’) to visualise the extent of spatial autocorrelation as a function of cumulative geographic distance. Statistical significance of r (H0 = lack of spatial genetic structure) was verified using 1000 permutations, and distance classes at which autocorrelation coefficients no longer remained significant were considered to approximate the true extent of identifiable genetic structure. Analyses were performed at smaller scale (distance range: 17–180 km), across samples from the Senegal River (Sen_DIAM, Sen_GUEC, Sen_MBO, Sen_NDO), and at larger scale, considering one sample for each country (Cam_BAK, Egy_BEK, Sen_DIAM, Cor_SUT; distance range: 2600–5400 km).
Maximum Likelihood (ML) tree reconstructions were performed on the un-partitioned data using the GTRCAT approximation of RAxML 8  and rapid bootstrap analysis based on 100 runs. In order to verify the consistency of tree reconstructions  ML analyses were repeated on three different datasets obtained by considering only polymorphic sites occurring in at least 30, 95 and 100% of specimens and including 6577, 2209 and 763 SNPs, respectively.
The DNA barcoding identification of specimens from this study relied on a reference library of 181 BOLD public records from 14 Bulinus species [B. truncatus; B. africanus (Krauss, 1848); B. barthi Jelnes, 1979; B. camerunensis Mandahl-Barth, 1957; B. cernicus (Morelet, 1867); B. forskalii (Ehrenberg, 1831); B. globosus (Morelet, 1866); B. nasutus (Martens, 1879); B. natalensis (Küster, 1841); B. nyassanus (Smith, 1877); B. senegalensis Müller, 1781; B. tropicus (Krauss, 1848); B. umbilicatus Mandahl-Barth, 1973; B. wrighti Mandahl-Barth 1965]. All of the B. truncatus queries but four (see Additional file 2: Table S2) could be unambiguously identified as B. truncatus. We consider that these four ambiguous identifications (Cam_BAK_1, Cam_BAK_7, Cam_BAK_9, Cam_BAK_10) were due to the presence of a misidentified/mislabelled sequence of B. tropicus in the reference dataset (GenBank accession KJ157497). We assume that this reference sequence was misidentified as: (i) it was the only B. tropicus sequence in an otherwise mono-specific clade of B. truncatus (83% bootstrap support) and (ii) all other B. tropicus reference sequences were recovered in two clades, one including three specimens from Cameroon (60% bootstrap support) and the other also including B. natalensis and B. nyassanus (76% bootstrap support; see Additional file 3: Figure S1).
Overall, low levels of parasite sequences were identified, as across all samples the median percentage of reads assigned to S. haematobium or S. bovis was as low as 1% (see Additional file 4: Figure S2). Levels of contamination were three- to ninefold higher in the two positive controls (2.8% and 9.3%, respectively, of the reads mapped to the parasite genomes). In addition, two samples of Sen_GUEC showed similar high incidence of contamination (5.2% and 6.9%, respectively, contaminant reads).
Allelic richness measured for Bulinus (Table 1) did not show marked variations across populations, with values ranging from 1.056 (Sen_GUEC) to 1.140 (Sen_DIAM). Inbreeding coefficients (Fis) were significantly different from zero in all samples, except in the population of Cameroon (CAM_BAK), notably the population exhibiting the smallest sample size (Table 1). Pairwise Fst values (Fig. 2) ranged from 0.018 (between the two samples from Egypt: Egy_MAN and Egy_BEK) to 0.876 (between the samples from Sen_GUEC and Cor_SUT) and were significantly different from zero in all pairwise comparisons except between Egy_MAN and Egy_BEK. The PCA of axes 1 and 2 (explaining 22.0% and 17.21% of the variance, respectively; Fig. 3) allowed resolving of the four main sample groups corresponding to specimens from: (i) Egypt (sites Egy_BEK, Egy_MAN); (ii) Cameroon and Corsica (sites Cam_BAK and Cor_Sut); (iii) Senegal (sites Sen_DIA, Sen_MBO, Sen_NDO); and (iv) Senegal (site Sen_GUEC). PCA of axes 1 and 3 (explaining 22.01% and 13.2% of the variance, respectively) allowed further resolving of the specimens from Cameroon and Corsica and grouped together all specimens from Senegal (Fig. 3). The STRUCTURE analysis showed ∆K values  peaking at K = 4 (see Additional file 5: Figure S3) and indicates that the main hierarchical level of the population structure is based on four groups. As this method detects the uppermost level of population structure when different hierarchical levels exist , we further explored the genetic sub-structuring of our dataset by arbitrarily considering K = 5 and K = 6. At K = 4 (Fig. 4), samples from Senegal were assigned to two separate groups, one including Sen_DIAM, Sen_MBO and Sen_NDO, and the other containing Sen_GUEC. The analysis could also resolve the two samples from Egypt (Egy_BEK, Egy_MAN) on one hand, and the sample from France (Cor_SUT) on the other hand. Conversely, specimens from Cameroon (Cam_BAK) did not show clear assignment patterns at K = 4, while they could only be resolved when considering K = 5. Increasing K to 6 did not provide further sample resolution.
The Mantel test showed that individual geographic and genetic distances were highly correlated (observed r = 0.679, p < 0.001). A discontinuous distribution of individual geographic and genetic distances was observed due to the two different spatial scales being sampled (small scale within the Senegal River Basin (SRB), large scale across countries) (Fig. 5). Positive spatial autocorrelation was detected across samples from the Senegal River at all geographic distances investigated (from 0–10 to 0–150 km); hence even at this scale gene flow between adjacent populations was not high enough to counter genetic drift. Similarly, at the larger spatial scale (across countries) positive autocorrelation could be observed at all distance classes, starting from 0 to 2500 km (Fig. 5).
The midpoint-rooted ML tree reconstructions (Fig. 6; Additional file 6: Figure S4, 7: Figure S5) consistently recovered five main and highly supported clades (bootstrap support = 100%), including specimens from: (i) Sen_DIAM, Sen_MBO, Sen_NDO; (ii) Sen_GUEC; (iii) Cor_SUT; (iv) Cam_BAK; and (v) Egy_BEK, Egy_MAN. The two samples from Egypt were resolved as a single group. Samples from the Senegal River were recovered as two sister groups, with Sen_GUEC in one group and Sen_DIAM, Sen_MBO, Sen_NDO in the other group. The sample from France (Cor_SUT) appeared to be more closely related to samples from Senegal than to samples from Egypt or Cameroon.
Population genetic structure of B. truncatus
The phylogenetic analysis based on 2209 SNPs identified five distinct clades, with samples from each country forming well-supported monophyletic groups (Fig. 6). In addition, all individual specimens could be unambiguously assigned to their country of origin (Fig. 4), indicating little gene flow between countries. Remarkably, the Senegalese populations form two distinct groups where the pairwise genetic distance between the Sen_GUEC population and the other Senegalese populations is as large as the pairwise genetic distance between clades from different countries (Fig. 2).
Consistent with previous studies [24, 58, 59], a significant genetic differentiation between populations could be observed at both a regional (within countries) and an inter-regional (between countries) scale. The genetic differentiation between the two populations from Egypt was low (0.022, not significant), and this was also true between three of the Senegalese populations (Sen_DIAM, Sen_MBO and Sen_NDO; Fst 0.140–0.256, significant difference), indicating relatively high amounts of gene flow on a regional scale of 100 km (cf. Gow et al. [23, 24] in Cameroonian B. truncatus populations). In Senegal, the low differentiation might also be a relic of the Diama dam construction (see below) that was followed by a rapid expansion of snail populations, leading to highly inbred and genetically similar snail populations in the lower valley of the SRB, as was found for Biomphalaria pfeifferi in the same area . However, there was a significant genetic differentiation (Fst 0.709–0.894) between these three populations and the Sen_GUEC population only 100 km further upstream, of a magnitude comparable to the genetic differentiation between the Senegalese and Egyptian populations (Fst 0.689–0.751), which are almost 5000 km apart. To interpret this differentiation of the Sen_GUEC population, a few factors influencing gene flow need to be considered. First, high stream velocities in rivers might flush away snails, either through detachment from their substrate or through exceedance of the snail’s tolerance limits, giving rise to gene flow from upstream to downstream populations [61, 62]. Secondly, gastropods disperse over large distances on the feet and feathers of birds [63, 64] or through humans . The combination of these phenomena can lead to the biotic homogenisation of snail populations over large distances [23, 61, 64, 65]. On the contrary, genetic differentiation on very small scales can also occur due to genetic drift. Populations with a small effective population size, either as a result of seasonal extinction events or frequent selfing, will experience more drift. These phenomena are typical of Bulinus species .
Considering these phenomena, the genetic dissimilarity of the Sen_GUEC population could be explained by its location along the Doué River, which discharges in the Senegal River approximately 60 km downstream and originates from the same river approximately 160 km upstream. It is possible that the current from the Doué River is an effective barrier for snails from the lower valley of the SRB to migrate upstream, while snails coming from the higher part of the Senegal River can migrate downstream easily. Furthermore, because it is a tributary of the Senegal River, the Doué River populations contribute relatively little to the gene pool in the lower SRB compared to the populations coming from the upper parts of the Senegal River. A similar pattern has been observed in Egypt, where populations in the lower Nile River are genetically similar to each other but are distinct from populations upstream . Finally, the sampling locality was a closed site, a bit further away from the river. Seasonal desiccation events might have induced strong bottleneck events as opposed to the other three sites that were situated on the riverbank.
A clear IBD pattern and positive spatial autocorrelation is present, both on a regional scale and on an inter-regional scale, even after omitting the Sen_GUEC population from the regional analysis. This finding is consistent with the results from microsatellite studies by Viard et al.  and Nyakaana et al.  on B. truncatus and B. globosus, respectively, and indicates that (regional) gene flow is not strong enough to completely counteract the effects of high selfing rates and genetic drift in B. truncatus [58, 68].
Genetic diversity of B. truncatus
The hermaphroditic nature of B. truncatus strongly influences the population genetic diversity. In this study we found an excess of homozygotes and low levels of genetic diversity in all populations, which is consistent with previous studies on highly selfing freshwater snails (e.g. [23, 24, 65,66,67, 69, 70]). However, variable inbreeding coefficients between geographic locations have been observed in our study (e.g. high in Senegal vs low in Egypt) and in previous studies (e.g. [58, 59, 65, 71, 72]). This variation might be attributed to different frequencies of population bottlenecks (e.g. in temporary habitats or after molluscicidal treatments), temperature-dependent phally polymorphism favouring aphallic individuals at high temperature (resulting in higher selfing rates in the population) and the higher relative importance of selfing in permanent habitats [24, 58, 66, 73,74,75]. In Senegal, high selfing rates might have facilitated the rapid expansion of B. truncatus after the completion of the Diama dam in 1985 and the Manantali dam on the Bafing River in 1989, which led to the creation of large open water bodies [76, 77], followed by a rapid expansion of the parasite. Such a scenario has been described for Bi. pfeifferi, where a microsatellite study showed a low genetic diversity for all snail populations from the lower valley of the SRB. According to the authors, this low diversity, in combination with the higher mobility of the parasite (mediated through humans) compared to the snail, resulted in the massive outbreak of intestinal schistosomiasis .
Viard et al. [58, 66] showed that B. truncatus populations in temporary habitats have a larger genetic diversity than populations in permanent habitats, despite the high selfing rates in both habitats. This is also true in our study where the slightly higher allelic richness of the Egyptian and Senegalese populations (except for the Sen_GUEC population) might be attributed to the temporary habitats in both countries. The Egyptian snails were sampled inside irrigation canals that dry up frequently, while the Senegalese snails originate from the Senegal, Lampsar and Doué River that are prone to large debit fluctuations  and from the irrigation canals linked with these same rivers. The good aestivation capacities of B. truncatus and the ability to quickly (re)colonise previously dried-out habitats [9, 79] result in genetic mixing of snails coming from different aestivation spots or refugia and the maintenance of a higher genetic diversity in temporary habitats [75, 80, 81]. In permanent habitats, the combination of high selfing rates and limited relative gene flow erodes the genetic diversity of B. truncatus populations .
The allelic richness values observed in this study were far below values observed in studies on Cameroonian and Nigerian B. truncatus populations ([57, 70] respectively). However, both of these studies made use of microsatellite markers that are characterised by a high level of polymorphism . SNP markers are much less polymorphic, with usually only two different alleles per locus. Another point of attention in our study is the low number of snails genotyped in some locations (e.g. the Cam_BAK population), which might result in an under- or overestimation of genetic diversity.
A note on snail population genetic diversity and parasite susceptibility
The population genetic structure of snails might affect their susceptibility to parasite infections. The low genetic diversity observed in all B. truncatus populations in this study in combination with the rapid adaptation potential of parasites might indicate that a parasite genotype able to infect the local snail genotype can easily cause a high infection prevalence in these snail populations [83, 84]. It has been shown that, in the SRB, the genetic differentiation pattern of Bi. pfeifferi populations corresponds to the pattern observed in S. mansoni populations [60, 85] and that these S. mansoni populations are highly adapted to the local Bi. pfeifferi populations [60, 86, 87]. These observations might explain the high S. mansoni transmission rates in the SRB . Similar to these findings, we found a clear genetic differentiation of the B. truncatus population from Sen_GUEC from the other Senegalese populations that corresponds well to the differentiation between the S. haematobium populations from the middle and lower valley of the SRB, as shown by Boon et al. . This genetic differentiation of both B. truncatus and S. haematobium possibly explains why B. truncatus acts as an intermediate host for S. haematobium in the upper valley while B. truncatus is not susceptible to S. haematobium in the lower valley . The middle valley (where we sampled Sen_GUEC) might be a transition zone between both regions, with the presence of both susceptible and non-susceptible B. truncations populations .
The GBS technique was able to confirm snail infections in samples that were previously identified to be infected during shedding experiments in the field, but it could also detect prepatent infections as it detected schistosome DNA in two samples that were negative after shedding. Although Preston et al.  showed that trematode biomass within Ramshorn snails (Helisoma trivolvis Say, 1817) comprised only about 3–4% of the total snail biomass, in our samples up to 9% of total reads could be attributed to schistosome parasites. Although the percentage of parasite biomass in a snail cannot directly be related to the percentage of parasite reads in that sample, this still shows that even blind HTS techniques can detect parasite DNA in snail tissue without prior amplification rounds, opening up new avenues for future research on the importance of genetics in determining snail-parasite compatibility.
In this study, five distinct B. trunatus clades could be identified, each comprising snails from a single country but with two distinct clades within Senegal. The genetic diversity within populations was low because of high selfing rates, but varied between geographic locations, probably due to habitat variability. Considerable genetic differentiation at both a regional and an inter-regional scale was observed, mainly due to the clear differentiation of the Sen_GUEC population from the other Senegalese populations. However, a clear IBD pattern remains present on both spatial scales, even after the Sen_GUEC population was omitted from the analysis. This indicates that regional gene flow is not strong enough to counteract the effects of bottlenecks, high selfing rates and genetic drift. The genetic differentiation patterns observed in B. truncatus are remarkably close to the patterns observed in S. haematobium in the area, which could explain the differences in parasite susceptibility of snails from the lower and middle valley of the SRB. HTS techniques offer a valuable toolbox to further investigate these patterns.
Availability of data and materials
The genotype data supporting the conclusions of this article are available on DataDryad (DOI:https://doi.org/10.5061/dryad.37pvmcvpc).
Barcode of life database
Cytochrome c oxidase subunit I
Isolation by distance
Principal coordinate analysis
Reduced representation library
Single nucleotide polymorphism;
Senegal River Basin
Adema CM, Bayne CJ, Bridger JM, Knight M, Loker ES, Yoshino TP, et al. Will all scientists working on snails and the diseases they transmit please stand up? PLoS Negl Trop Dis. 2012;6:5–6.
Brown DS, Shaw KM. Freshwater snails of the Bulinus truncatus/tropicus complex in Kenya: tetraploid species. J Molluscan Stud. 1989;55:509–32.
Jorgensen A, Madsen H, Nalugwa A, Nyakaana S, Rollinson D, Stothard JR, et al. A molecular phylogenetic analysis of Bulinus (Gastropoda: Planorbidae) with conserved nuclear genes. Zool Scr. 2011;40:126–36.
Jørgensen A, Jørgensen LVG, Kristensen TK, Madsen H, Stothard JR. Molecular phylogenetic investigations of Bulinus (Gastropoda: Planorbidae) in Lake Malawi with comments on the topological incongruence between DNA loci. Zool Scr. 2007;36:577–85.
Brown D. Freshwater snails of Africa and their medical importance. Boca Raton: Taylor & Francis; 1994.
Cowie RH, Dillon RT, Robinson DG, Smith JW. Alien non-marine snails and slugs of priority quarantine importance in the United States: a preliminary risk assessment. Am Malacol Bull. 2009;27:113–32.
Maes T, Hammoud C, Volckaert FAM, Huyse T. A call for standardised snail ecological studies to support schistosomiasis risk assessment and snail control efforts. Hydrobiologia. 2021;848:1773–93.
Najarian HH. Biological studies on the snail, Bulinus truncatus, in central Iraq. Bull World Health Organ. 1961;25:435–46.
Watson JM. Ecology and distribution of Bulinus truncatus in the Middle East; with comments on the effect of some human activities in their relationship to the snail host on the incidence of Bilharziasis haematobia in the Middle East and Africa. Bull World Health Organ. 1958;18:833–94.
Picquet M, Ernould JC, Vercruysse J, Southgate VR, Mbaye A, Sambou B, et al. The epidemiology of human schistosomiasis in the Senegal river basin. Trans R Soc Trop Med Hyg. 1990;90:340–6.
Abdel-Wahab MF, Yosery A, Narooz S, Esmat G, El Hak S, Nasif S, et al. Is Schistosoma mansoni replacing Schistosoma haematobium in the Fayoum? Am J Trop Med Hyg. 1993;49:697–700.
Barakat R, El Morshedy H, Farghaly A. Human schistosomiasis in the Middle East and North Africa region. In: Mary AM, Sima R, editors. Neglected tropical diseases—Middle East and North Africa. Neglected tropical diseases series. Vienna: Springer; 2014. p. 23–57.
Hotez PJ, Savioli L, Fenwick A. Neglected tropical diseases of the Middle East and North Africa: review of their prevalence, distribution, and opportunities for control. PLoS Negl Trop Dis. 2012;6:1475.
Stensgaard AS, Vounatsou P, Sengupta ME, Utzinger J. Schistosomes, snails and climate change: current trends and future expectations. Acta Trop. 2019;190:257–68.
Boissier J, Grech-Angelini S, Webster BL, Allienne J-F, Huyse T, Mas-Coma S, et al. Outbreak of urogenital schistosomiasis in Corsica (France): an epidemiological case study. Lancet Infect Dis. 2016;16:971–9.
Githeko A, Lindsay S, Confalonieri U, Patz J. Climate change and vector-borne diseases: a regional analysis. Bull World Health Organ. 2000;78:1136–47.
Berry A, Moné H, Iriart X, Mouahid G, Aboo O, Boissier J, et al. Schistosomiasis haematobium, Corsica. Fr Emerg Infect Dis. 2014;20:1595–7.
Martínez-Ortí A, Vilavella D, Bargues MD, Mas-Coma S. Risk map of transmission of urogenital schistosomiasis by Bulinus truncatus (Audouin, 1827) (Mollusca Gastropoda, Bulinidae) in Spain and Portugal. Anim Biodivers Conserv. 2019;42:257–66.
Mulero S, Rey O, Arancibia N, Mas-Coma S, Boissier J. Persistent establishment of a tropical disease in Europe: the preadaptation of schistosomes to overwinter. Parasit Vectors. 2019;12:379.
Stensgaard A-S, Booth M, Nikulin G, McCreesh N. Combining process-based and correlative models improves predictions of climate change effects on Schistosoma mansoni transmission in eastern Africa. Geospat Health. 2016;11:406.
Coustau C, Gourbal B, Duval D, Yoshino TP, Adema CM, Mitta G. Advances in gastropod immunity from the study of the interaction between the snail Biomphalaria glabrata and its parasites: a review of research progress over the last decade. Fish Shellfish Immunol. 2015;46:5–16.
Rollinson D, Webster JP, Webster B, Nyakaana S, Jrgensen A, Stothard JR. Genetic diversity of schistosomes and snails: implications for control. Parasitology. 2009;136:1801–11.
Gow JL, Noble LR, Rollinson D, Tchuenté LAT, Jones CS. Contrasting temporal dynamics and spatial patterns of population genetic structure correlate with differences in demography and habitat between two closely-related African freshwater snails. Biol J Linn Soc. 2007;90:747–60.
Gow JL, Noble LR, Rollinson D, Mimpfoundi R, Jones CS. Breeding system and demography shape population genetic structure across ecological and climatic zones in the African freshwater snail, Bulinus forskalii (Gastropoda, Pulmonata), intermediate host for schistosomes. Mol Ecol. 2004;13:3561–73.
Jarne P, Théron A. Genetic structure in natural populations of flukes and snails: a practical approach and review. Parasitology. 2001;123:27–40.
Adema CM, Hillier LDW, Jones CS, Loker ES, Knight M, Minx P, et al. Whole genome analysis of a schistosomiasis-transmitting freshwater snail. Nat Commun. 2017;8:1–12.
Elshire RJ, Glaubitz JC, Sun Q, Poland JA, Kawamoto K, Buckler ES, et al. A robust, simple genotyping-by-sequencing (GBS) approach for high diversity species. PLoS ONE. 2011;6:1–10.
Hebert PDN, Cywinska A, Ball SL, DeWaard JR. Biological identifications through DNA barcodes. Proc R Soc B Biol Sci. 2003;270:313–21.
Folmer O, Black M, Hoeh W, Lutz R, Vrijenhoek R. DNA primers for amplification of mitochondrial cytochrome c oxidase subunit I from diverse metazoan invertebrates. Mol Mar Biol Biotechnol. 1994;3:294–9.
Kane RA, Stothard JR, Emery AM, Rollinson D. Molecular characterization of freshwater snails in the genus Bulinus: a role for barcodes? Parasit Vectors. 2008;1:15.
Kearse M, Moir R, Wilson A, Stones-Havas S, Cheung M, Sturrock S, et al. Geneious basic: an integrated and extendable desktop software platform for the organization and analysis of sequence data. Bioinformatics. 2012;28:1647–9.
Sonet G, Jordaens K, Nagy ZT, Breman FC, De Meyer M, Backeljau T, et al. Adhoc: an R package to calculate ad hoc distance thresholds for DNA barcoding identification. Zookeys. 2013;365:329–35.
Meier R, Shiyang K, Vaidya G, Ng PKL. DNA barcoding and taxonomy in diptera: a tale of high intraspecific variability and low identification success. Syst Biol. 2006;55:715–28.
Kimura M. A simple method for estimating evolutionary rates of base substitutions through comparative studies of nucleotide sequences. J Mol Evol. 1980;16:111–20.
Tamura K, Stecher G, Peterson D, Filipski A, Kumar S. MEGA6: molecular evolutionary genetics analysis version 6.0. Mol Biol Evol. 2013;30:2725–9.
Boon NAM, Mbow M, Paredis L, Moris P, Sy I, Maes T, et al. No barrier breakdown between human and cattle schistosome species in the Senegal River Basin in the face of hybridisation. Int J Parasitol. 2019;49:1039–48.
Catchen JM, Amores A, Hohenlohe P, Cresko W, Postlethwait JH. Stacks: building and genotyping loci de novo from short-read sequences. G3 Genes Genomes Genet. 2011;1:171–82.
Danecek P, Auton A, Abecasis G, Albers CA, Banks E, DePristo MA, et al. The variant call format and VCFtools. Bioinformatics. 2011;27:2156–8.
Keenan K, McGinnity P, Cross TF, Crozier WW. Diversity: an R package for the estimation and exploration of population genetics parameters and their associated errors. Methods Ecol Evol. 2013;4:782–8.
Excoffier L, Lischer HEL. Arlequin suite ver 3.5: a new series of programs to perform population genetics analyses under linux and windows. Mol Ecol Resour. 2010;10:564–7.
Guo SW, Thompson EA. Performing the exact test of Hardy–Weinberg proportion for multiple alleles. Biometrics. 1992;48:361.
Wright S. Coefficients of inbreeding and relationship. Am Nat. 1922;56:330–8.
Weir BS, Cockerham CC. Estimating F-statistics for the analysis of population structure. Evolution. 1984;38:1358–70.
Benjamini Y, Hochberg Y. Controlling the false discovery rate: a practical and powerful approach to multiple testing. J R Stat Soc Ser B. 1995;57:289–300.
Jombart T, Devillard S, Dufour A-B, Pontier D. Revealing cryptic spatial patterns in genetic variability by a new multivariate method. Heredity. 2008;101:92–103.
Pritchard JK, Stephens M, Donnelly P. Inference of population structure using multilocus genotype data. Genetics. 2000;155:945–59.
Evanno G, Regnaut S, Goudet J. Detecting the number of clusters of individuals using the software structure: a simulation study. Mol Ecol. 2005;14:2611–20.
Earl DA, VonHoldt BM. Structure harvester: a website and program for visualizing Structure output and implementing the Evanno method. Conserv Genet Resour. 2012;4:359–61.
Jakobsson M, Rosenberg NA. CLUMPP: a cluster matching and permutation program for dealing with label switching and multimodality in analysis of population structure. Bioinformatics. 2007;23:1801–6.
Rosenberg NA. Distruct a program for the graphical display of population structure. Mol Ecol Notes. 2004;4:137–8.
Smouse PE, Peakall R. Spatial autocorrelation analysis of individual multiallele and multilocus genetic structure. Heredity. 1999;82:561–73.
Peakall R, Ruibal M, Lindenmayer DB. Spatial autocorrelation analysis offers new insights into gene flow in the Australian bush rat, Rattus fuscipes. Evolution. 2003;57:1182–95.
Peakall ROD, Smouse PE. genalex 6: genetic analysis in excel population genetic software for teaching and research. Mol Ecol Notes. 2006;6:288–95.
Peakall R, Smouse PE. GenAlEx 6.5: genetic analysis in excel population genetic software for teaching and research-an update. Bioinformatics. 2012;28:2537–9.
Stamatakis A. RAxML version 8: a tool for phylogenetic analysis and post-analysis of large phylogenies. Bioinformatics. 2014;30:1312–3.
Rivers DM, Darwell CT, Althoff DM. Phylogenetic analysis of RAD-seq data: examining the influence of gene genealogy conflict on analysis of concatenated data. Cladistics. 2016;32:672–81.
Coulon A, Fitzpatrick JW, Bowman R, Stith BM, Makarewich CA, Stenzler LM, et al. Congruent population structure inferred from dispersal behaviour and intensive genetic surveys of the threatened Florida scrub-jay (Aphelocoma cœrulescens). Mol Ecol. 2008;17:1685–701.
Viard F, Bremond P, Labbo R, Justy F, Delay B, Jarne P. Microsatellites and the genetics of highly selfing populations in the freshwater snail Bulinus truncatus. Genetics. 1996;142:1237–47.
Nalugwa A, Jorgensen A, Nyakaana S, Kristensen TK. Genetic variation within and between populations of hermaphroditic Bulinus truncatus tetraploid freshwater snails of the Albertine Rift, East Africa. Hydrobiology. 2011;673:53–61.
Campbell G, Noble LR, Rollinson D, Southgate VR, Webster JP, Jones CS. Low genetic diversity in a snail intermediate host (Biomphalaria pfeifferi Krass, 1848) and schistosomiasis transmission in the Senegal River Basin. Mol Ecol. 2010;19:241–56.
Habib MR, Lv S, Rollinson D, Zhou X. Invasion and dispersal of Biomphalaria species : increased vigilance needed to prevent the introduction and spread of schistosomiasis. Front Med. 2021;8:1–16.
Dazo BC, Hairston NG, Dawood IK. The ecology of Bulinus truncatus and Biomphalaria alexandrina and its implications for the control of bilharziasis in the Egypt-49 project area. Bull World Health Organ. 1966;35:339–56.
Neubauer TA, Mandic O, Harzhauser M, Jovanović G. The discovery of Bulinus (Pulmonata: Planorbidae) in a Miocene palaeolake in the Balkan Peninsula. J Molluscan Stud. 2017;83:295–303.
Pfenninger M, Salinger M, Haun T, Feldmeyer B. Factors and processes shaping the population structure and distribution of genetic variation across the species range of the freshwater snail Radix balthica (Pulmonata, Basommatophora). BMC Evol Biol. 2011;11:135. https://doi.org/10.1186/1471-2148-11-135.
Zein-Eddine R, Djuikwo-Teukeng FF, Dar Y, Dreyfuss G, Van den Broeck F. Population genetics of the Schistosoma snail host Bulinus truncatus in Egypt. Acta Trop. 2017;172:36–43.
Viard F, Justy F, Jarne P. Population dynamics inferred from temporal variation at microsatellite loci in the selfing snail Bulinus truncatus. Genetics. 1997;146:973–82.
Nyakaana S, Stothard JR, Nalugwa A, Webster BL, Lange CN, Jorgensen A, et al. Bulinus globosus (Planorbidae; Gastropoda) populations in the Lake Victoria basin and coastal Kenya show extreme nuclear genetic differentiation. Acta Trop. 2013;128:226–33.
Rousset F. Genetic differentiation and estimation of gene flow from F-statistics under isolation by distance. Genetics. 1997;145:1219–28.
Lounnas M, Correa AC, Vazquez AA, Dia A, Escobar JS, Nicot A, et al. Self-fertilization, long-distance flash invasion and biogeography shape the population structure of Pseudosuccinea columella at the worldwide scale. Mol Ecol. 2017;26:887–903.
Jarne P, Viard F, Delay B, Cuny G. Variable microsatellites in the highly selfing snail Bulinus truncatus (Basommatophora: Planorbidae). Mol Ecol. 1994;3:527–8.
Djuikwo-Teukeng FF, Njiokou F, Nkengazong L, De Meeus T, Ekobo AS, Dreyfuss G. Strong genetic structure in Cameroonian populations of Bulinus truncatus (Gastropoda: Planorbidae), intermediate host of Schistosoma haematobium. Infect Genet Evol. 2011;11:17–22.
Chlyeh G, Henry PY, Sourrouille P, Delay B, Khallaayoune K, Jarne P. Population genetics and dynamics at short spatial scale in Bulinus truncatus, the intermediate host of Schistosoma haematobium, in Morocco. Parasitology. 2002;125:349–57.
Betterton C. Spatiotemporal distributional patterns of Bulinus rohlfsi (Clessin), Bulinus forskali (Ehrenberg) and Bulinus senegalensis (Muller) in newly-irrigated areas in northern Nigeria. J Molluscan Stud. 1984;50:137–52.
Schrag SJ, Read AF. Temperature determination of male outcrossing ability in a simultaneous hermaphrodite. Evolution. 1992;46:1698–707.
Jarne P. Mating system, bottlenecks and genetic polymorphism in hermaphroditic animals. Genet Res. 1995;65:193–207.
Baker HG. Self-compatibility and establishment after “long-distance” dispersal. Evolution. 1955;9:347–9.
Pannell JR, Auld JR, Brandvain Y, Burd M, Busch JW, Cheptou P, et al. The scope of Baker’s law. New Phytol. 2015;208:656–67.
Sakho I, Dupont JP, Cisse MT, El JS, Loum S. Hydrological responses to rainfall variability and dam construction: a case study of the upper Senegal River basin. Environ Earth Sci. 2017;76:1–19.
Chu KY, Bijan H, Massoud J. The ability of Bulinus truncatus, Biomphalaria alexandrina and Lymnaea gedrosiana to survive out of water in the laboratory. Ann Trop Med Parasitol. 1967;61:1–5.
Mavárez J, Amarista M, Pointier JP, Jarne P. Fine-scale population structure and dispersal in Biomphalaria glabrata, the intermediate snail host of Schistosoma mansoni, in Venezuela. Mol Ecol. 2002;11:879–89.
Mukaratirwa S, Siegismund HR, Kristensen TK, Chandiwana SK. Population genetics and genetic variability of Bulinus globosus (Gastropoda: Planorbidae) from the two main river systems in Zimbabwe. J Hered. 1996;87:288–94.
Zane L, Bargelloni L, Patarnello T. Strategies for microsatellite isolation: A review. Mol Ecol. 2002;11:1–16.
King KC, Lively CM. Does genetic diversity limit disease spread in natural host populations. Heredity. 2012;109:199–203.
Gandon S, Capowiez Y, Dubois Y, Michalakis Y, Olivieri I. Local adaptation and gene-for-gene coevolution in a metapopulation model. Proc R Soc B Biol Sci. 1996;263:1003–9.
Van den Broeck F, Maes GE, Larmuseau MHD, Rollinson D, Sy I, Faye D, et al. Reconstructing colonization dynamics of the human parasite Schistosoma mansoni following anthropogenic environmental changes in northwest Senegal. PLoS Negl Trop Dis. 2015;9:e0004090.
Southgate VR, Tchuem Tchuenté LA, Théron A, Jourdane J, Ly A, Moncrieff CB, et al. Compatibility of Schistosoma mansoni Cameroon and Biomphalaria pfeifferi Senegal. Parasitology. 2000;121:501–5.
Diaw OT, Vassiliades G, Seye M, Sarr Y. Epidemiology of intestinal schistosomiasis with Schistosoma mansoni in Richard-Toll (delta of the Senegal River). Malacological study. Bull Soc Pathol Exot. 1991;84:174–83 (in French).
Talla I, Kongs A, Verlé P. Preliminary study of the prevalence of human schistosomiasis in Richard-Toll (the Senegal river basin). Trans R Soc Trop Med Hyg. 1992;86:182.
Sène M, Southgate VR, Vercruysse J. Bulinus truncatus, hôte intermédiaire de Schistosoma haematobium dans le bassin du fleuve Sénégal. Bull Soc Pathol Exot. 2004;97:29–32.
Southgate VR, De Clercq D, Sène M, Rollinson D, Ly A, Vercruysse J. Observations on the compatibility between Bulinus spp. and Schistosoma haematobium in the Senegal river basin. Ann Trop Med Parasitol. 2000;94:157–64.
Preston DL, Orlofske SA, Lambden JP, Johnson PTJ. Biomass and productivity of trematode parasites in pond ecosystems. J Anim Ecol. 2013;82:509–17.
The authors wish to sincerely thank the field team (M. Wade and F. Van den Broeck) who helped with snail collections in Senegal, and Prof. Dr. Jerome Boissier of the University of Perpignan–Via Domitia for providing the snail samples from Corsica.
We would like to acknowledge the financial support of the Belgian Science Policy (BELSPO) through JEMU (Joint experimental molecular unit). Tim Maes benefited from an FWO fellowship grant (ref. no. 1S86319N) of the Research Foundation—Flanders.
Ethics approval and consent to participate
Consent for publication
The authors have no conflicts of interests to declare.
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Sample list and information.
Best close match DNA barcoding identification (1% distance threshold) of vouchers considered in this study.
Neighbour Joining (NJ) tree representing K2P genetic distances among vouchers considered in this study (blue dots) and 181 public reference sequences (http://www.boldsystems.org/) from 14 Bulinus genera (B. truncatus, B. africanus, B. barthi, B. camerunensis, B. cernicus, B. forskalii, B. globosus, B. nasutus, B. natalensis, B. nyassanus, B. senegalensis, B. tropicus, B. umbilicatus, B. wrighti).
Percentage of reads assigned to Schistosoma haematobium contamination per snail specimen.
Log-likelihood probability values (LnP(D)) and ∆K (according to Evanno et al. 2005) as obtained in STRUCTURE with K ranging from 1 to 8 (value obtained by averaging the posterior probabilities of three independent runs).
Maximum likelihood tree reconstruction based on 763 SNPs recovered in the 100% of specimens considered in this study (see text for explanations).
Maximum likelihood tree reconstruction based on 6577 SNPs recovered in at least 70% of specimens considered in this study (see text for explanations).
About this article
Cite this article
Maes, T., De Corte, Z., Vangestel, C. et al. Large-scale and small-scale population genetic structure of the medically important gastropod species Bulinus truncatus (Gastropoda, Heterobranchia). Parasites Vectors 15, 328 (2022). https://doi.org/10.1186/s13071-022-05445-x
- Bulinus truncatus
- Population structure
- Isolation by distance