The population genetic structure of Biomphalaria choanomphala in Lake Victoria, East Africa: implications for schistosomiasis transmission

Sampling methodology and population selection

Snails were collected from 29 sites around the perimeter of Lake Victoria, using hand held scoops (Figure 1). Sampling was semi-quantified, with two collectors each surveying a 50 m length of shoreline for approximately 20 minutes and collecting all Biomphalaria found. Sites were selected to ensure an even geographical spread along the lakeshore, in all three countries and in both marsh and lake habitats; five localities were specifically included because they incorporated both habitat types in close proximity. All sites showed signs of human activity near the shoreline. In order to compare sites of similar population density, all sites selected had a high abundance of Biomphalaria, defined as 30 or more individuals collected, which accounted for the paucity of sites along the eastern shoreline of the lake in Tanzania, where Biomphalaria densities were generally low.

Snails were putatively identified as B. choanomphala or B. sudanica, based on shell morphology, as these are the species thought primarily to inhabit Lake Victoria [15]. However, recent taxonomic work has suggested that all the Biomphalaria in Lake Victoria can be considered a single species of which ecophenotypy of shell form gives rise to the previously considered taxonomic distinctions [20], and thus for the purposes of this study, individuals were simply classed as being one of the two morphotypes of B. choanomphala.

All collected snails (in total >30 per site, given the site selection criteria) were exposed, as a group, to sunlight for 2–4 hours to check for infection with Schistosoma mansoni or other parasite cercariae [21]. A random sample of snails was examined in this way on subsequent days, to ensure maximum reliability of infection testing. After shedding, snails were placed in glass tubes containing 95% ethanol solution.

DNA extraction, amplification and sequencing

Genomic DNA was extracted from between 10 and 12 individuals from each of the selected 29 populations, using a standard CTAB extraction [22] with a final re-suspension in pure water. The mitochondrial cytochrome oxidase sub-unit one (COI) gene and the mitochondrial 16S sub-unit of the ribosomal RNA gene were amplified using the Folmer ‘universal’ primers [23] and modified 16ar and 16br primers [24], respectively. The four microsatellite loci were Bpf1, Bpf2, Bpf3 and Bpf10 [25] - all nine of the primers described in referenced paper were tested, but the above four were the only ones which amplified successfully and were polymorphic in the Lake Victoria Biomphalaria choanomphala under examination here. Mitochondrial (COI and 16S) amplifications were carried out in 25 μl total volume, with 2.5 μl MgCl₂ (20 mM concentration), 2.5 μl 5 × buffer, 2.5 μl pre-mixed dNTPs (20 mM concentration), 1 μl each of forward and reverse primer (10 pmol concentration) and one unit of TAQ per reaction. Microsatellite reactions were amplified in 12.5 μl total volume, with half of the volumes stated above, and with forward primer concentration of 50 pmol and reverse (fluorescently labeled) primer concentration of 20 pmol. Cycling conditions in all cases followed published methods [23]-[25].

Positively amplified products were purified using Millipore PCR₉₆ Cleanup kits on a vacuum manifold (Millipore, Billerica, USA) as per manufacturer’s instructions, using pure water for washing and re-suspension. Product concentration was quantified on a Nanodrop ND-1000 Spectrophotometer (Nanodrop Technologies Inc., Willington, USA), and sequencing reactions were performed on mitochondrial purified PCR products using an Applied Biosystems Big Dye Kit (version 1.1) and run on an Applied Biosystems 3730 DNA Analyzer (Applied Biosystems, Carlsbad, USA). Microsatellites were diluted 1:10 in HiDiformamide (Applied Biosystems, as above) and analyzed using an ABI3730 automated sequencer, using a GeneScan 500 LIZ size standard (both Applied Biosystems, as above).

Population analyses

The overall framework of the population genetics analyses were designed to examine the population structure, diversity, demography, phylogeography and spatial correlations of the B. choanomphala populations surveyed from Lake Victoria [2].

COI and 16S mitochondrial sequences were aligned using MUSCLE on the web (http://www.ebi.ac.uk/Tools/msa/muscle/) and edited visually in MacClade v 4 [26]. The aligned sequences were then compared and reduced to unique haplotypes, and the proportion of each haplotype per population was recorded. The 16S data contained many gap-rich regions, for which it was impossible to be assured of homology; any region with a gap was therefore removed from the alignment before calculating genetic distances amongst sequences. There was only one indel region in the COI alignment, of a single codon repeat, in three haplotypes, which was therefore easily alignable and retained for haplotype analysis.

Measures of population differentiation were calculated from mitochondrial sequence data using analysis of molecular variance (AMOVA), F-statistics, mismatch analysis, Tajima and Fu’s tests for neutrality and a Mantel z-test for correlation by distance. Each AMOVA was conducted using “country” as the primary grouping, to investigate the level of variation partitioned between each country, versus between populations within each country and within each population. “country” was deemed to be a potentially important geographical factor given that each country in the region undertakes different levels of schistosomiasis control initiatives, thus potentially influencing the level of parasitism experienced by local snails. All tests were executed in Arlequin v 3.2 [27].

A distance-based, neighbor-joining tree was created in Paup* 4.0 [28] for the unique haplotypes that remained when indels were removed from the alignment from each gene, using the GTR + gamma model of sequence evolution, determined to be the most appropriate model through likelihood testing, also in Paup* 4.0. Median-joining networks were constructed using Network 4.6 [29] based on the gap-free unique haplotypes.

The number and proportion of private microsatellite alleles per population was calculated using the web-based version of Genepop [30]. Calculations of F_ST and tests for Hardy-Weinberg equilibrium (through estimation of F_IS values per population and per locus) and linkage disequilibrium amongst loci were also done using Genepop. The loci were also tested for evidence of previous bottlenecking events by comparing the allelic diversity against the observed and expected heterozygosity, using the program Bottleneck [31].

Mitochondrial data

The COI and 16S sequences were very genetically variable with 127/308 (CO1) and 181/300(16S) unique haplotypes/total number of sequences respectively. Haplotypes were not highly divergent (0.002-0.058 pairwise distance for COI haplotypes and 0.000-0.047 for 16S). GenBank accession numbers for these datasets are HM769132-HM769258 for COI and HM768950-HM768980 and HM768982-HM769131 for 16S. When unalignable regions with gaps were removed from the analysis of the 16S sequences, the number of unique haplotypes was reduced from 181 to 64. Many common haplotypes were spread throughout the lake, but others were locally restricted and there was high variation in the haplotype richness of the sites: overall, intra-site haplotype diversity estimates varied from 0.200-0.982 for COI and 0.378-1.000 for 16S, where zero indicates complete homogeneity and one indicates maximum diversity.

Figure 2 shows the geographical distribution of the five most abundant COI and 16S (“A” and “B” maps in Figure 2, respectively) haplotypes (without removal of indels) per population across Lake Victoria, as compared to the frequency of other shared haplotypes and “private” haplotypes (i.e. haplotypes only observed at one site). The map reveals patterns of higher and lower haplotype sharing versus uniqueness within and between localities and countries; for example, Kenyan sites tended to share haplotypes whereas Ugandan sites, especially those in the central and western regions of the lakeshore, had high proportions of private haplotypes. Maps showing the distribution of all the haplotypes for COI and 16S, plus a full list of haplotype frequencies per site, can be found in the supplementary information (See Additional file 1: Figure S1, Additional file 2: Figure S2, Additional file 3: Table S1 and Additional file 4: Table S2 for COI and 16S respectively).

The phylogenetic trees for COI and 16S (Figure 3) were consistent, as would be expected from linked loci, and showed few deep divisions, but several more recent processes of divergence. These supported the hypothesis of geographical structuring, as divergent clades tended to be made up of haplotypes from the same country, or limited to two out of the three countries. The networks for each mitochondrial marker were difficult to analyze visually due to the large number of haplotypes (See Additional file 5: Figure S3 for COI and Additional file 6: Figure S4 for 16S). However there appeared to be evidence of geographical clustering by country; for both markers, the most abundant and geographical widespread haplotype (H1 for COI and H2 for 16S; the network software utilized the gap-free dataset) was a central node within each respective network.

When pair-wise fixation index (F_ST) values were calculated between each of the populations in turn, the majority demonstrated significant levels of differentiation, for both types of molecular marker. The range in values was 0.000-0.921 for COI and 0.000-0.845 for 16S. For the 16S data, all but one of the non-significant F_ST values were between sites in the same country, once again demonstrating geographically defined structuring of the populations for this gene. All but one paired site showed significant differentiation between the populations from lake versus marsh habitat, indicating even small-scale geographical segregation. Full tables of the pairwise F_ST values and their significance can be found in the supplementary information (See Additional file 7: Tables S3 and Additional file 8: Tables S4 respectively for COI and 16S).

All of the populations conformed to the assumptions of the neutral hypothesis, as tested for by Tajima’s D statistic for the 16S data; one population, T001, deviated significantly (p <0.001) from assumptions of normality when the COI data was analyzed. However, this population also had a significant result for the sum of squares mismatch analysis (p <0.001), suggesting that demographic change, such as population expansion, may have biased the neutrality test. No other populations, for either marker, came out as significant in the mismatch analysis.

The relationship with geography was elucidated through the significance of the Mantel z-test for both the COI and 16S data (p =0.001, for both), meaning that the genetic and geographical distances between the populations was positively correlated; in other words, sites closer together were more likely to be genetically similar, even if the composition of distinct haplotypes differed per site (as seen in Figure 3 and Additional file 1: Figure S1). This demonstrated an overall relationship between genetic differentiation and geographic distance, on top of the localized differentiation already described, such as between marsh and lake habitats (Additional file 5: Figure S3 and Additional file 6: Figure S4).

Microsatellite data

The microsatellite data were tested for departures from Hardy-Weinberg equilibrium (HWE), linkage disequilibrium and for an excess of heterozygosity relative to allelic diversity, which would be consistent with populations having recently passed through a bottleneck.

There was also no evidence for linkage disequilibrium between the loci, suggesting that the alleles considered are not closely associated in the genome and are mixing randomly during gamete production. Four populations (K013b, T001, T027b and U028), when tested for evidence of a past bottleneck, showed heterozygote excess (p = 0.031 for all four), but this value was not significant once the Bonferroni correction for multiple tests was applied.

Population structure of Biomphalaria choanomphala in Lake Victoria

The populations studied here were characterized by high intrapopulation diversity as well as high levels of population structure, with low levels of gene flow as inferred from molecular data and negligible inbreeding. High diversity is consistent with what would be expected if parasitism were strongly influencing population structure [3],[7], and the low F_IS values and apparent random mating observed here further supports this hypothesis. The snails collected for this research were examined for infection with schistosomes and other parasites in the field as part of a different study; however, without more accurate and consistent methods such as molecular probes [32] for the detection of infected snails (and indeed, infection with potentially more than one parasite species [33]) and infected snail populations, it is difficult to determine whether parasitism is indeed a key driver of population structure in this system. Similarly, future studies should seek to use genetic markers which are known to correlate to parasite resistance or susceptibility, which were not available when this study was conducted, but have since been identified in several studies [34],[35]. Even so, the challenges of modifying markers for use on field populations of less well-studied African Biomphalaria remain.

We would also expect parasite populations to adapt to local snails, thus mirroring to some extent their population structure. We used a Mantel test [36] to compare the population structure of B. choanomphala from the six sites from which there were published data on schistosome population structure [19], and found no correlation. Notwithstanding the limitation of the small number of sites, this conforms to other studies where parasite populations were much less differentiated than those of the intermediate host, suggesting local compatibility in this context is being overshadowed by other factors such as high parasite migration [36],[37]. Definitive host migration has been put forward as a factor that might homogenize parasite populations across a wide area, and this is compatible with the Lake Victoria context, where human populations are known to move widely throughout the region [18].

Role of other factors influencing population structure

Given the possibility of parasitism as a driver for population structure in Lake Victoria snail populations, it is important to consider what other factors might be influential in this system, and particularly how these differ from previous studies of Biomphalaria in other settings. Such considerations link back to parasitism, and particularly control of schistosomiasis in human populations, in terms of gaining a better understanding of drivers for distribution and abundance of host snails under different conditions. Studies of B. pfeifferi in Madagascar, for example, revealed high levels of interpopulation variation but, in contrast to what we observed here, consistently low levels of intrapopulation genetic differentiation [11]. This was explained as being due to the combination of low levels of migration combined with habitat stochasticity; B. pfeifferi are also known to be frequent in-breeders, which would assist in maintaining a metapopulation with seasonal local extinction events [38].

In contrast, Lake Victoria possesses very different environmental conditions; as a very large, permanent lake, it rarely experiences the kind of environmental perturbations that characterize temporally transient ponds and streams. Such homogeneity of habitat, with associated infrequent extinctions and opportunity for genetic drift, could account for the high levels of intrapopulation diversity that are maintained in Lake Victorian B. choanomphala[2]. This applies mainly to the lake habitats; indeed, this study showed significantly lower levels of gene and nucleotide diversity in the marsh habitats. However, in both lake and marsh populations, the tests for bottleneck events were uniformly negative, suggesting that sudden demographical events are not overwhelmingly contributing to the patterns of diversity of Lake Victoria Biomphalaria.

Effect on transmission of intestinal schistosomiasis

Although the low population structure in S. mansoni observed previously in this region counters the hypothesis of highly specialized local parasite compatibility [19] for both host and parasite, certain B. choanomphala haplotypes were widespread around the lake. This suggests that, as seen in B. pfeifferi from elsewhere in Africa, there may be a risk of transmissive genotypes becoming widely successful, if there is selection for traits independent of parasite susceptibility [39].

Another factor to consider is human-mediated changes to parasite populations: all three countries bordering Lake Victoria have initiated efforts to control schistosomiasis in human populations through mass treatment campaigns, although these initiatives are not synchronized and differ in their degree of coverage. The low levels of population structuring seen in the parasite [19] are likely a function of high terminal host migration and rapid dispersal of the parasite into novel localities [36],[40]; if some populations of S. mansoni experience more intense selection pressure based on frequent mass drug administration campaigns, and given the generally transmissive snail populations, there is a risk that these highly adapted forms of S. mansoni will be able to spread rapidly throughout the region. This emphasizes the need for coordinated control strategies, both at the level of human treatment as well as improved water access and sanitation, throughout lakeshore communities in the region.