Evidences That Host Genetic Background More Than the Environment Shapes the Microbiota of the Snail Bulinus truncatus, an Intermediate Host of Schistosoma Species

TL;DR

This study shows that the genetic background of snails has a stronger influence on their microbiota than the environment, in the context of Schistosoma parasites.

Contribution

The study introduces new microsatellite markers and provides evidence for host genetics as a key driver of microbiota composition in wild snail populations.

Findings

Snail population genetics and geographic distribution strongly influence microbiota composition.

Environmental bacterial communities have a weaker but significant effect on snail microbiota.

Shell size and trematode infection status do not significantly impact microbiota structure.

Abstract

Microbiota have emerged as fundamental regulators of host physiology, shaping both ecological interactions and evolutionary trajectories. Yet, the determinants of microbiota diversity and structure in wild populations—particularly the respective roles of host genetics and environmental context—are still poorly understood. In this study, we investigated these influences in the freshwater snail Bulinus truncatus, a key intermediate host for human and animal Schistosoma parasites, using a multifactorial approach. We developed 31 new microsatellite markers to resolve population genetic structure across nine sites in Senegal. Metabarcoding methods were then employed to profile the bacterial microbiota of individual snails and to characterize environmental bacterial assemblages from each location via environmental DNA. Shell measurements and molecular diagnostics for trematode infection status were included to assess additional potential contributors. Employing multiple regression on distance matrices (MRM), we quantified how snail population genetics, site‐specific environmental bacterial communities, spatial patterns, and infection status shape microbiota composition. Our analyses reveal that snail geographic distribution and population genetic structure drive the composition of Bulinus truncatus microbiota, with environmental bacterial communities exerting a weaker but still significant effect. In contrast, neither shell size nor trematode infection status impacted microbiota structure significantly. Notably, a considerable fraction of variation remains unexplained, indicating the likely involvement of other ecological or intrinsic factors. These results advance understanding of microbiota determinants in natural populations and underscore the intricate interplay between host genetics, environment, and microbial communities.

Figures

Locations of the studied aquatic transmission sites (black dots) nearby the closest established villages, on a satellite map of northern Senegal.

(A) Relative abundance of the ten most abundant phyla in the microbiota of B. truncatus from the nine different sites. (B) Relative abundance of the ten most abundant phyla in the environmental samples from the nine different sites. (C) Bacterial composition of the core microbiota of B. truncatus snails at the family level across the 9 studied sites.

Principal coordinate analysis (PCoA) plot showing environmental samples (triangles) and B. truncatus snail associated (dots) bacterial communities distributed along the two first PCoA axes computed based on pairwise Jaccard distances obtained between each pair of samples.

Principal coordinate analysis (PCoA) plot showing the bacterial communities of B. truncatus according to their original sampling sites (see colour legend), and to their infection status by trematodes based on a molecular diagnostic (dot size) distributed along the two first PCoA axes based on pairwise Jaccard similarity indices computed between each sample. Each coloured small dot represents one snail host from a given sampling site. Large dots are snail hosts positive to trematode molecular diagnostic. Medium dots are the centroids of the bacterial communities of snails originating from each sampling site. Each line represents the connection between the sample and the centroids of the bacterial communities of snails originating from each sampling site.

Distance‐decay scatterplots showing bacterial communities distance (determined by the pairwise Jaccard similarity index) as a function of geographic distance across sampling sites for environmental samples (A), and B. truncatus snails (B). Scatterplot of B. truncatus snails microbiota dissimilarity (determined by the Jaccard index) as a function of Nei's genetic distance between B. truncatus individuals computed based on 31 microsatellite loci (C). Blue lines of best fit are included to highlight trends.

Tables

TABLE 1. Genetic diversity statistics for B. truncatus at each sampling site based on the 31 microsatellite markers developed in this study. N a : Number of alleles; A r : Allelic richness; H e : Gene diversity (Nei 1978); H o : Observed heterozygosity; F is: Inbreeding coefficient; p‐values of F is after 10,000 randomizations of gene copies among individuals.

Site	Sample size	N a	A r	H e	H o	F is	p (F is ≠ 0)
Ndiawara	7	2.52	2.42	0.405	0.382	0.063	0.028
Ouali Diala	11	2.42	2.27	0.375	0.389	−0.042	0.054
Guia	16	2.58	2.21	0.384	0.406	−0.06	0.001
Dioundou	5	2.13	2.13	0.369	0.396	−0.087	0.011
Fonde Ass	15	3.06	2.53	0.417	0.39	0.069	0.001
Khodit	15	2.97	2.63	0.453	0.396	0.133	0
Mbane	12	3.23	2.75	0.401	0.41	−0.023	0.239
Saneinte	16	3.26	2.68	0.383	0.382	0.004	0.806
Lampsar	13	2.9	2.57	0.421	0.392	0.072	0
Overall sites	110	5.52	3.18	0.499	0.394	0.211	0

Full text

Introduction

1

The ‘microbiota’, consisting of all microorganisms associated with a host, is now broadly acknowledged as a crucial factor in various aspects of organismal biology, encompassing developmental, physiological, reproductive, behavioural and immunological phenotypes (Liberti and Engel 2020; Daisley et al. 2020; Yuan et al. 2021; Davidson et al. 2024). The microbiota therefore partly influences the ecology and the evolution of their hosts and that of biotic interactions between host species (Sharon et al. 2010; Henry et al. 2021; Lange et al. 2023). Despite recent efforts to understand the factors shaping host microbiota, identifying and understanding the ecological and evolutionary determinants of microbiota composition and structure in natural host populations remains a major scientific challenge.

The composition and structure of the microbiota of hosts in natural environments are particularly complex and depend on multiple environmental and host genetic factors (Benson et al. 2010). On one hand, the microbiota is partly composed of some either obligatory or facultative symbionts that are transmitted, generally via maternal inheritance, from one host generation to another (Funkhouser and Bordenstein 2013). Their evolutionary history hence follows that of their hosts with which they have co‐evolved at both macroevolutionary and microevolutionary timescales leading to phylosymbiosis (Kohl 2020; Dapa et al. 2023). These microorganisms constitute the ‘core microbiota’ which consists of species consistently associated with the host regardless of environmental conditions (Risely 2020). Accordingly, the composition and structure of hosts' microbiota is expected to be partly shaped by a combination of stochastic (e.g., dispersal, drift) and deterministic (e.g., selection) processes associated with hosts ecology and evolution (Benson et al. 2010; Furman et al. 2020; Hayashi et al. 2024).

On the other hand, most microorganisms that constitute hosts' microbiota are facultative and are acquired from the environment throughout the host's lifelong development. This is well illustrated by some species of Lepidopteran that lack a resident microbiota and generally acquire bacterial communities from their host plants or from the environmental microbial communities (Montagna et al. 2016; Phalnikar et al. 2018; Minard et al. 2019; Liu et al. 2020). Moreover, several environmental biotic factors including the communities of free‐living microorganisms present in the environment (Díaz‐Sánchez et al. 2018) and abiotic factors (e.g., temperature, pH) can influence the composition of the hosts' microbiota, which can in turn modify hosts' biology (Bernardo‐Cravo et al. 2020). In other words, the ecology of hosts (e.g., behaviour, diet) also influences their microbiota (Archie and Tung 2015; Kennedy et al. 2020). Thus, it is expected that organisms from a common genetic pool and established in different environments harbour different communities of microorganisms (Berg et al. 2016). Our comprehension of factors shaping hosts' microbiota hence requires accounting for the ecological context in which hosts are established and the genetic background of host populations.

Surprisingly however, few studies have yet specifically investigated the relative contribution of organisms' genetic diversity and that of their environment in shaping natural hosts' microbiota and lead to distinct conclusions. For instance, Suzuki and collaborators found that the structure of gut microbiota of wild house mouse ( Mus musculus domesticus ) populations is well predicted by genetic distances between host populations computed from an extensive exome genomic dataset and not by different environmental conditions, including temperature and diet (Suzuki et al. 2019). Conversely, the structure of the gut microbiota associated with California voles ( Microtus californicus ) populations across a contact zone between two recently diverged lineages was best explained by the spatial distribution of hosts and not by lineage divergence, hence suggesting a strong influence of the environment on the structure of voles' gut microbiota (Lin et al. 2020). In the same vein, Rothschild et al. (2018) found that the genetic background among 1046 healthy human individuals has a minor role in determining individuals' gut microbiota while environmental factors such as housing, diet and anthropometric measurements can explain some of the inter‐individual variability in microbiota (Rothschild et al. 2018). Finally, in the Nematostella vectensis sea anemone model, Fraune et al. have shown that both the environment and genetics play a role in shaping the microbiota, and that these elements are complexly interconnected (Fraune et al. 2016). In the light of these studies, it appears that the relative weight of the genetic or environmental determinants of microbiota are potentially species and/or context dependent.

Bulinus truncatus constitutes the main intermediate host for several trematode species including Schistosoma haematobium, the etiological agent for urogenital bilharziasis in human populations in Africa (Toledo and Fried 2014). It is also a major intermediate host for Schistosoma bovis, a widespread livestock parasite responsible for bovine schistosomiasis across Africa, with important veterinary and economic impacts (Demlew and Tessma 2020). This snail species is tolerant to a wide range of abiotic conditions, has a high passive dispersal capacity, and is capable of self‐fertilization (Jarne et al. 1993). As a result, this species is widely distributed across Africa, sometimes at high abundance, and can establish in a wide variety of freshwater ecosystems including rivers, lakes, (temporal) ponds, and reservoirs (Tumwebaze et al. 2022). Due to its ability to proliferate across diverse environments, this species serves as a valuable model for exploring the relative impact of snail genetics and environmental factors on the diversity and structure of its microbiota.

Ecological studies aimed at characterizing the communities of microorganisms associated with natural snail populations have recently flourished, especially in freshwater snails that are associated with disease transmission (Li et al. 2023). Based on culture‐independent molecular methods, Van horn et al. have shown highly diverse intestinal microbial communities among three freshwater planorbid snails collected from the field. Two of these snail species serve has intermediate hosts of several parasites including digenetic trematodes of the genus Schistosoma, responsible for human Schistosomiasis (Van Horn et al. 2012). Similarly, Huot et al. demonstrated that seven species and strains within the Planorbidae exhibited highly specific core bacterial composition that are closely related to the host strain. This bacterial composition show high congruence with the host phylogeny, hence revealing a phylosymbiotic pattern (Huot et al. 2020). More recently, McCann et al. showed that the microbiota of Galba truncatula , the main intermediate host for the zoonotic trematode Fasciola hepatica , vary between natural populations established at two geographically and ecologically distinct sites, hence highlighting the existence of natural diversity in the composition of bacterial communities associated with this snail species (McCann et al. 2024).

However, whether such microbiota diversity in natural populations is driven by snails' genetic background or by environmental factors remains to be fully elucidated. If the freshwater snails' microbiota is important in the circulation of associated pathogens, as has been suggested (Portet et al. 2021; Clec'h et al. 2022), it is essential to study the natural diversity, the structure and the factors structuring these communities in the field.

Here, we took advantage of high‐throughput sequencing technologies to identify and assess the relative contribution of the ecological and host genetic factors structuring the microbiota associated with natural populations of the freshwater snail B. truncatus .

We characterized the diversity and structure of the whole‐body microbiota of B. truncatus across nine freshwater aquatic sites in Northern Senegal, a well‐documented region endemic for S. haematobium and S. bovis (Léger et al. 2020). The genetic background of each snail host was characterized based on 31 newly developed microsatellite markers using a genotyping‐by‐sequencing approach. We also measured snails' shell size and checked for the presence of trematodes developing within each snail using a 16S trematode metabarcoding assay (Douchet et al. 2022). Aquatic bacterial communities associated with each of the nine studied freshwater sites were also characterized and used as a proxy of the ecological context in which B. truncatus populations are established. We used these complementary datasets to identify the main factors structuring the microbiota associated with B. truncatus in natural populations.

Materials and Methods

2

Field Collection of Samples

2.1

Field sampling was conducted in February 2022 during the dry season in the Senegal River basin region (Figure 1). We focused on nine natural sites of S. haematobium transmission, where urogenital schistosomiasis' prevalence among school‐aged children was previously reported (Senghor et al. 2022). These sites differ in terms of ecological contexts and in B. truncatus abundance (Douchet et al. 2024). Six transmission sites are located near the villages of Ndiawara, Ouali Diala, Dioundou, Fonde Ass and Khodit, all of which being located along river “le Doue”, a tributary of the Senegal River in the middle valley. Sites along the river “le Doue” varied in terms of aquatic plant cover, substrate size and quantity of artificial substrates (e.g., plastics), although we did not quantify these variables. One transmission site, near the village of Guia, is also located in the middle valley and consists of an irrigation canal that drains water from the river “le Doue”. In this canal the water flow was low and the vegetation cover (both bottom and surface) high. Two transmission sites, near the villages of Mbane and Saneinte are located along the east shore of Guiers lake. The Mbane site had intermediate vegetation cover and rapidly increasing depth while the Saneinte site had shallow depth, little vegetation, and numerous (semi)‐emerged anthropogenic substrates (e.g., plastic containers, tires). One transmission site near the village of Lampsar is located in the lower valley and consists of an inlet of the Senegal River delta (Table S1). This site consists of a gently sloping beach surrounded by dense aquatic vegetation. All these nine sites are frequented by nearby human populations (e.g., for washing clothes, dishes, and bathing) and/or considered as transmission sites for S. haematobium. Pictures of these sites are provided in Figure S1.

Locations of the studied aquatic transmission sites (black dots) nearby the closest established villages, on a satellite map of northern Senegal.

Field Collection of Environmental Samples

2.1.1

At each site, we first filtered water from the surface to the water–sediment interface along the shores (i.e., max depth 20 cm). We used filtration capsules of 0.45 μM mesh size (Waterra USA Inc.) connected to an electric water pump as previously described in Douchet et al. (2022). Filtration was performed until the filtration capsule was clogged and the final filtration volume was retrieved. We here consider that filtering until capsule clogging allows standardizing the same amount of particles collected across all samples and maximizing the DNA detection. Once the filtration completed, the filtration capsule was depleted from its water content, filled with 50 mL of Longmire solution (100 mM Tris, 100 mM EDTA, 10 mM NaCl, 0.5% SDS, 0.2% sodium azide, (Sahu et al. 2025)), vigorously shaken, and preserved in the dark at ambient temperature until used for environmental DNA (eDNA) extractions. At each site but Khodit (due to logistical supply issue), a technical field negative control was obtained by filtering 1.5 L of commercial spring water using the same 0,45 μM filter capsule as for eDNA sampling, following the same protocol for preservation (i.e., using 50 mL of Longmire solution once the filter capsule depleted from any remaining water).

Field Collection and Preservation of Snails' Samples

2.1.2

Following eDNA sampling, we manually collected all encountered aquatic snails by scooping the aquatic vegetation using a colander for about 30 min to 1 h, along a 10 to 30‐m long transect along the shore of the targeted waterbody. After taxonomic identification (following the keys of (Mandahl‐Barth 1962)), up to 17 non‐emitting B. truncatus individuals per site were extracted from their shell after a brief heating step at 70°C for 1 min, using decontaminated forceps, and individually preserved in sterile 1.5 mL tubes containing 70% ethanol. Forceps were decontaminated between each sample using a DNA AWAY solution. The extraction of snail bodies from their shells was achieved to avoid DNA contamination with microorganisms, including bacteria, established on snail shells. The length and width of each individual shell were measured to the nearest tenth of a centimetre using a calliper (see Table S2 and Figure S2). Overall, 124 snails were prepared as above.

DNA Extraction From eDNA Samples and

B. truncatus Snails

2.2

All pre‐PCR steps including DNA extraction and the preparation of PCR mixes were conducted under a sterile hood decontaminated before and after each use with 10% bleach, 70% ethanol, a DNA AWAY solution followed by UV light exposure for 20 min.

Total eDNA from water filtrations were extracted following Douchet et al. (Douchet et al. 2022). Briefly, the Longmire buffer contained within each filtration capsule was equally split into three 50 mL tubes. For the field negative controls (i.e., spring water filtrates), the Longmire solution from each filtration capsule was recovered in one single 50 mL tube. The tubes were centrifuged for 20 min at 16,000 × g and the supernatant was discarded. We then collected a similar amount of sediment from the pellet from each tube (less than 1 g as recommended in the protocol). Thus, total eDNA from each capsule was extracted in triplicate. This strategy allowed us to capture sediment originating from the full range of microhabitats relevant to the species while ensuring replication within the molecular workflow. No pellets were observed for filtration negative controls after centrifugation. For these samples, 500 μL of Longmire were retained and resuspended after removing the supernatant and were processed as the other samples. This step led to the processing of 35 samples (i.e., 3 extraction replicates for each of the nine eDNA samples and one field negative control per site except for Khodit, see Figure S3 for a graphical vision of the processed samples). The total eDNA of these 35 samples was extracted using the DNeasy PowerSoil Pro kit (Qiagen) according to manufacturer's protocol performing the physical lysis with a MagNA Lyser at a speed of 7000 × g for 30s.

Total genomic DNA from each of the 124 snails removed from their shell was extracted using DNeasy 96 Blood and Tissue Kit (Qiagen), according to the manufacturer's protocol. Two extraction controls were also added in the process. Snail assignation to B. truncatus was confirmed via a molecular diagnostic tool based on specific lamp amplification first developed to detect B. truncatus eDNA from water samples (Blin et al. 2023) and recently applied to snails (Douchet et al. 2024).

Infection Status and Genotyping of

B. truncatus Snails

2.3

Infection Status of

B. truncatus Snails

2.3.1

The presence of developing trematodes within each of the 124 snails was diagnosed using the Trem_16S_F1 (GACGGAAAGACCCCRAGA) and Trem‐16S_R2 (CRCCGGTYTTAACTCARYTCAT) 16S trematode metabarcode (Douchet et al. 2022). DNA extracts were diluted to 1/100th to minimize PCR inhibitors and hence limit PCR false negatives. PCR reactions were prepared and run under standard conditions (Douchet et al. 2022). After this first screening step, a new PCR was achieved using DNA templates from positive samples and sent to the Bio‐Environment platform (University of Perpignan Via Domitia, France) for metabarcoding sequencing to identify the trematodes potentially developing within snail hosts and possible coinfections. At the platform, individual libraries were performed using the Illumina Nextera index kit and Q5 high fidelity DNA polymerase (New England Biolabs). Indexed PCR products were then normalized with SequalPrep plates (ThermoFisher) and paired‐end sequenced on a MiSeq instrument using v2 chemistry (2 × 250 bp). The obtained sequences were processed following Douchet et al. (2022) and blasted on the NCBI database for molecular identification as in Douchet et al. (2022, 2024). OTUs were assigned to species when they showed ≥ 97% coverage and ≥ 99% pairwise identity with a reference sequence. OTUs showing ≥ 97% coverage and ≥ 96% identity were assigned to genus. OTUs that did not meet these thresholds were assigned to higher taxonomic ranks using a clustering approach based on pairwise genetic distances with a set of 174 reference sequences.

Microsatellite Development and Sequence‐Based Microsatellite Genotyping

2.3.2

The B. truncatus genome (GenBank accession GCA_021962125.1, (Young et al. 2022)) was used for microsatellite discovery, leading to the selection of 96 primer pairs optimized for multiplex PCR (details in Appendix S1). Multiplex PCR amplification and sequencing libraries were constructed, with detailed methods provided in Appendix S2.

Genetic Analyses

2.3.3

Preliminary analyses led to the selection of 31 among the 96 initially developed microsatellite markers to avoid markers of bad quality and have no missing genetic data overall individuals. Descriptive statistics of genetic diversity including the mean number of alleles, the mean allelic richness (computed based on the minimal number of genotypes obtained over sampling sites, N = 5), the mean expected (He) and observed (Ho) heterozygosity, and the mean Fis were computed at each sampling site using the SPAGeDi software (v. 1.5d; (Hardy and Vekemans 2002)) accounting for the allotetraploidy nature of B. truncatus (Njiokou et al. 1993).

We next computed pairwise Nei's genetic standard distance Ds (Nei 1978) and F st values between each sampling sites using the SPAGeDi software (v. 1.5d; (Hardy and Vekemans 2002), Table S3). From the obtained D _ s _ values, we ran a Mantel test to test for possible isolation by distance pattern among sampling sites using matrices of pairwise D _ s _ values and pairwise linear geographical distances obtained between each sampling site. The statistical significance of the Mantel coefficient r was tested based on 10,000 permutations as implemented in the ‘adegenet’ R package (Jombart 2008). A similar Mantel test analysis was run using pairwise (F _ st _/(1—F _ st _)) values obtained between sampling sites.

Sequencing of Bacterial Communities and Metabarcoding Analyses

2.4

A total of 198 bacterial 16S metabarcoding libraries were prepared following the Illumina two‐step PCR protocol, including the DNA extracts from the 124 snails, 70 (35 × 2) eDNA replicates consisting of 54 eDNA samples (i.e., each extraction triplicate from each of the nine sampled sites was amplified in duplicate) and 16 technical field negative controls (each control was amplified in duplicate), and 4 negative PCR controls (milliρ water) (Table S4). We targeted the variable V3‐V4 loops region of the 16S rDNA gene using the 341F (5′‐CCTACGGGNGGCWGCAG‐3′) and 805R (5′‐GACTACHVGGGTATCTAATCC‐3′) primers (Klindworth et al. 2013) combined with universal Illumina adapters.

The first PCRs were performed using the Q5 High‐Fidelity 2× Master Mix (New England BioLabs), in a 25 μL final volume containing 2 μL of template DNA and using a PCR program consisting of an initial denaturation step of 30 s at 98°C followed by 32 cycles containing a denaturation step of 6 s at 98°C, an annealing step of 30 s at 55°C, and an elongation step of 8 s at 72°C and ending with a final elongation step of 60 s at 65°C. We used 5 μL of the PCR products to check the quality and integrity of amplicons through electrophoresis on agarose gels. The 20 μL left were sent to the Bio‐Environment platform (University of Perpignan Via Domitia, France). Libraries were performed using the Illumina Nextera index kit and Q5 high fidelity DNA polymerase (New England BioLabs). Indexed PCR products were then normalized with SequalPrep plates (ThermoFisher) and paired‐end sequenced on a MiSeq instrument using v2 chemistry (2 × 250 bp).

Metabarcoding raw data obtained from sequencing were processed using R (version 2023.03.0 + 386) using a pipeline based on the ‘dada2’ package (v. 1.28.0) (Callahan et al. 2016) on the IFB (French Institute of Bioinformatics) cloud. Briefly, this pipeline consists of removing primers from the obtained sequenced reads (max.mismatch = 1), trimming and quality filtering reads (minLen = 150, maxN = 0, maxEE = c (3, 3), truncQ = 2), denoising, dereplicating, merging paired‐end reads (maxMismatch = 0), building the Amplicon Sequence Variant (ASV) table and removing potential chimeric sequences produced during the process (method = ‘consensus’). The taxonomic assignment (multithread = TRUE, minBoot = 60) of the resulting filtered and cleaned ASVs was performed using the Ribosomal Database Project (RDP) classifier (Wang et al. 2007) and based on the Silva_train_set and silva_species_assignment 138.1 datasets (Quast et al. 2012). Eukaryotic, mitochondrial, and chloroplast sequences were removed, as well as ASVs not affiliated to ‘Bacteria’ at the Kingdom taxonomic levels. The obtained ASV sequence table was finally merged with the taxonomy database and the sample metadata matrix for subsequent analyses using the ‘phyloseq’ R package (v. 1.22.3) (McMurdie and Holmes 2013).

ASV present at low abundance (< 0.005%; (Bokulich et al. 2013)) were filtered from the whole dataset to account for possible sequencing errors. Briefly, the whole phyloseq object was transformed in relative abundance with the ‘transform_sample_count’ function of the ‘phyloseq’ R package. Then, the low abundant taxa were removed through the ‘prune_taxa’ function of the ‘phyloseq’ R package. To account for the effect of potential contamination, a relative abundance filtering was performed by setting a filter to remove ASVs with an abundance inferior to 0.1% from each sample (Karstens et al. 2019). To further assess potential contamination, we compared the ASVs identified in PCR negative controls with those detected in biological samples. Two Phyloseq objects were built: one for the four PCR negative controls and another with B. truncatus samples. All ASVs present in the negative controls (with taxonomy details provided in Table S5) were screened for co‐occurrence in B. truncatus samples. For each shared ASV, we compared its read abundance in biological versus control samples to evaluate whether its presence was consistent with low‐level background contamination (see Figure S4A). The same comparative approach was applied to environmental samples (see Figure S4B).

To assess whether sequencing depth was satisfactory across samples, rarefaction curves were generated using a ‘ggrare’ custom function. The abundance of each ASV was normalized based on the lowest sample in terms of ASV counts (i.e., 5000, see Table S6) by random sub‐sampling using the ‘rarefy_even_depth’ function of the ‘phyloseq’ package (rngseed = 1000) (McMurdie and Holmes 2013). We chose to rarefy our data in order to provide a balanced view of the most abundant microbial communities *across B. truncatus

• populations. This approach aimed to minimize the influence of low‐abundance taxa that could introduce excessive noise in the data. Consequently, our analysis highlights the core microbiota shared among populations, while potentially underrepresenting rarer taxa present at lower abundances. To quantify read counts per sample during cleaning steps, refer to Table S6.

Statistical Analyses

2.5

Diversity in Environmental and Snail‐Associated Bacterial Communities and Characterization of

B. truncatus Core Microbiota

2.5.1

The diversity of bacterial communities associated with environmental matrices and B. truncatus snails was assessed at the sampling site level based on the absolute number of ASVs, the Shannon and Simpson's diversity indices as implemented in the ‘microbiome’ R package (v. 1.22.0) (Lahti and Shetty 2012). Kruskal‐Wallis nonparametric tests were computed to compare each of these alpha‐diversity indices between sample type (i.e., snails versus eDNA samples) and between sites using the ‘stats’ R package (v. 4.3.0) (R Core Team 2023). Those tests were followed by pairwise multiple comparisons (i.e., Dunn test, (Dunn 1964)) using the ‘dunnTest’ of the ‘FSA’ R package (v. 0.9.5) and adjusting the resulting p‐values according to Benjamini‐Hochberg (Benjamini and Hochberg 1995). To check whether snails' size could influence the alpha diversity of their microbiota we ran a linear mixed‐effect model (‘lmer function’; lme4 R package, (Bates et al. 2015)) using the number of ASV as the dependent variable and the length (or width) of snails' shell as an explanatory variable and fixing the sampling site as a random factor. To visualize the top ten most abundant phyla associated with B. truncatus , we aggregated taxonomic ranks at the phylum level and selected the 10 most abundant taxa based on overall relative read abundance. Relative abundances were normalized using the ‘transform_sample_counts’ function in the ‘phyloseq’ package (McMurdie and Holmes 2013). Barplots representing the mean composition across sampling sites were generated using the ‘plot_composition’ function from the ‘microbiome’ package (Lahti and Shetty 2012). Phylum identities were manually curated and mapped to the corresponding ASVs for figure display. The same process was applied to environmental samples.

To characterize the core microbiota of B. truncatus in natura, which consists in the most stable part of the bacterial communities shared among all individuals (Risely 2020), we followed established guidelines (Neu et al. 2021). The core microbiota was assessed at the family level using the ‘core_members’ function of the ‘microbiome’ R package (Lahti and Shetty 2012) and setting the detection parameter to 0 and the prevalence parameter to 90%.

Identification of Factors Structuring Environmental and Snail Bacterial Communities

2.5.2

To investigate for possible structuration of bacterial communities in environmental samples and B. truncatus snails we conducted a Principal Coordinate Analysis (PCoA) using the ‘pco’ function of the ‘ecodist’ package (v. 2.1.3) (Goslee and Urban 2007). Pairwise Jaccard and Bray‐Curtis indices were conducted using the ‘distance’ function of the ‘phyloseq’ R package. Based on these beta diversity indices, we ran permutational multivariate analyses of variance (PERMANOVA) to test for: 1‐ Differences between aquatic bacterial communities (i.e., from eDNA samples) and snails associated bacterial communities (i.e., microbiota). We accounted for the sampling site as an additional explanatory variable and for possible interaction between the sampling site and the nature of the sample (i.e., environment versus snail). 2‐ An effect of trematode infection status (all trematode species combined) on B. truncatus microbiota Beta‐diversity. Here we accounted for the sampling site as an additional explanatory variable and an interaction term between the sampling site and the infection status. PERMANOVA were conducted using the ‘adonis2’ function implemented in the ‘vegan’ R package (Oksanen et al. 2022), and setting the number of permutations to 10,000.

We investigated whether differences in snail shell size influence the beta diversity of the snails' microbiota at the site level. To this end, we conducted a series of partial Mantel tests to assess the correlation between pairwise matrices of absolute differences in shell size and microbiota beta diversity (Jaccard and Bray‐Curtis indices) while accounting for genetic distance among individual snails from each site independently. Those partial Mantel tests were conducted using the ‘mantel. partial’ function from the ‘vegan’ R package (Oksanen et al. 2022). p‐values were adjusted based on the Benjamini and Hochberg method (Benjamini and Hochberg 1995) (Table S7).

To assess the relative effect of the environmental bacterial communities and the genetic background of snails on whole‐body snail microbiota, we used four pairwise matrices, including the geographical distance matrix between each sampling site, the Jaccard indices computed between environmental bacterial communities at each sampling site, the Jaccard indices computed between bacterial communities associated with B. truncatus individuals and the Nei's genetic distance between snail hosts. To ensure comparable dimensions for the four matrices, 15 samples from the original matrices were discarded using the ‘subset_samples’ function from the ‘phyloseq’ R package. This step led to two matrices of 102*102 dimensions. Based on those matrices, we pseudoreplicated the geographic distance and the Jaccard indices computed between environmental bacterial communities at each sampling site matrices to transform them in the same dimensions. We conducted a multiple regression on distance matrices (MRM) fixing the Jaccard distance matrix obtained from the bacterial communities between B. truncatus individuals as the dependent variable and using three matrices as explanatory variables: the Jaccard distance between bacterial communities at each site (pseudo‐replicated at the individual level), the Nei's genetic distance computed between individuals and the Euclidean geographic distance between sites (also pseudo‐replicated at the individual level). This MRM analysis was ran using the ‘MRM’ function of the ‘ecodist’ R package and the significance of this multiple regression test was assessed based on 9999 permutations. This approach was complemented by a series of Mantel tests between these four matrices using the ‘mantel’ function of the ‘vegan’ R package that were visually represented using a custom function to plot matrices against one another.

Results

3

Abundance, Size and Infection Status of

B. truncatus Among Sites

3.1

Overall, 124 B. truncatus snails were collected among the nine sampled sites (mean = 14, min = 5, max = 17, Table S1). Mean snail shell size per site ranged from 3.30 mm in Dioundou to 8.01 mm in Saneinte (Table S2). The total trematode infection prevalence ranged from 0% (Dioundou, Mbane, Saneinte) to 29.4% (Fonde Ass) (Table S1). Nine species of trematodes were identified. Aside from S. haematobium and Schistosoma bovis, we also found Orientocreadium batrachoides, two other unidentified species of the genus Haematolochus and Petasiger, and two species of Paramphistomoidea and two species of Diplostomoidea. Among the 18 infected snails from the overall sites, 7 were co‐infected (1 in Ndiawara; 3 in Ouali Diala; 1 in Guia; 1 in Khodit; 1 in Lampsar).

Genetic Diversity and Structure of

B. truncatus Among Sampling Sites

3.2

A total of 110 snails were fully genotyped at 31 newly developed microsatellite markers (Table 1). The remaining samples did not yield DNA of sufficient concentration or quality for reliable sequencing‐based genotyping. Based on these individual genotypes we found little genetic diversity with little variation in genetic diversity among sites. Allelic richness (Ar) ranged from 2.13 (Dioundou) to 2.75 (Mbane) and observed heterozygosity (Ho) ranged from 0.38 (Ndiawara and Saneinte) to 0.41 (Guia and Mbane) (Table 1). Inbreeding coefficients (F is) were different from zero for all but the two sites Mbane and Saneinte located on the Guiers lake. Pairwise Ds genetic distances computed between sites ranged from 0.0019 between Ndiawara and Dioundou, both sites being located along the Doué river approximately 4.8 km one from each other; and 0.4354 between Ouali Diala (Doué River) and Lampsar, distant by approximately 160.3 km. The observed genetic differentiation between sites followed an isolation by distance pattern (Mantel r = 0.76; p‐value = 0.0001, Figure S5).

**TABLE 1: Genetic diversity statistics for B. truncatus at each sampling site based on the 31 microsatellite markers developed in this study. N

a : Number of alleles; A

r : Allelic richness; H

e : Gene diversity (Nei 1978); H

o : Observed heterozygosity; F is: Inbreeding coefficient; p‐values of F is after 10,000 randomizations of gene copies among individuals.**

Diversity in Environmental and Snail‐Associated Bacterial Communities and Characterization of

B. truncatus Snails' Core Microbiota

3.3

The bacterial communities from the environment displayed equal alpha diversity among sites except for Lampsar that displayed lower alpha diversity whatever the index used (i.e., Species Richness, Shannon, Simpson) (Figure S6). Moreover, B. truncatus microbiota displayed lower diversity than aquatic bacterial communities, except for the Lampsar site, when considering the ASV richness (Kruskall‐Wallis chi squared = 91.13; p‐value < 0.001), the Shannon diversity index (Kruskall‐Wallis chi squared = 95.32; p‐value < 0.001) and the Simpson index (Kruskall‐Wallis chi squared = 97.41; p‐value < 0.001) (Figure S6). We found significant differences in the mean alpha diversity indices in the microbiota of B. truncatus from the different sample sites (all p‐value < 0.05, Table S8 and Figure S7A). Pairwise post hoc tests showed site‐specific differences in alpha diversity (see Tables S9–S11).

We found no evidence that snails' shell size influences the alpha diversity (number of ASVs) of their microbiota according to our linear mixed‐effect model when accounting for length (t value = −1.598; p‐value = 0.113) or width (t value = −0.757; p‐value = 0.451). Similar results were obtained when using the Shannon index. Moreover, we found no significant differences in alpha diversity between infected and uninfected snails (all p‐value > 0.05, Table S12 and Figure S7B).

A total of 30 phyla of bacteria were found in overall samples including environmental DNA and snails, 18 of which (i.e., 60%) were common to eDNA and snail samples, and 12 (i.e., 40%) specifically found in the environment (Figure S8). Among the 10 most represented phyla identified in the environmental samples, seven were common to those found in the snail microbiota (Figure 2). The three phyla found in B. truncatus are Patescibacteria, Deinococcota and Fusobacteriota which ranked fourth, fifth and nineth in terms of relative abundance overall snails (Figure 2A). The three phyla present only in the environmental samples are the Verrucomicrobiota, Desulfobacterota and Acidobacteriota which were found in all environmental samples and ranked sixth, nineth and tenth overall sites in terms of relative abundance (Figure 2B).

(A) Relative abundance of the ten most abundant phyla in the microbiota of B. truncatus from the nine different sites. (B) Relative abundance of the ten most abundant phyla in the environmental samples from the nine different sites. (C) Bacterial composition of the core microbiota of B. truncatus snails at the family level across the 9 studied sites.

Proteobacteria were predominant in the microbiota of B. truncatus across all sampling sites, followed by Bacteroidota and Firmicutes, which relative abundance differed between sites (Figure 2A).

The core microbiota of B. truncatus was characterized at the family level by considering bacterial families present in at least 90% of the overall snails. The core microbiota consisted of 6 families out of the 134 overall families identified in the snails although their relative abundances varied (Figure 2C). These families included Chitinophagaceae (73.3% of core microbiota reads), Oxalobacteraceae (3.3%), Sphingomonadaceae (2.4%), Weeksellaceae (14.1%), Moraxellaceae (3.7%), and Flavobacteriaceae (3.2%).

Identification of Factors Structuring Environmental and Snail Bacterial Communities

3.4

The Microbiota of

B. truncatus Snails Differ From Environmental Bacterial Communities

3.4.1

Based on the microbial communities characterized from the 54 environmental samples and the 118 B. truncatus individuals (initially 124, with sample loss during rarefaction), we found that the total microbiota associated with whole‐body B. truncatus differs from bacterial communities present in the environment. Bulinus truncatus and environmental samples cluster into two distinct groups along the first axis which captures 18.27% of the total variance (Figure 3). For each matrix (i.e., snail and environment), samples then cluster according to the sampling site along the second axis which captures 6.36% of the total variance. (Figure 3). Such differences are confirmed by the PERMANOVA analysis which indicates significant effect of the type of sample (i.e., B. truncatus snail and environmental sample) (F‐value = 50.5; p‐value < 0.001) and of the sampling site (F‐value = 6.81; p‐value < 0.001). Similar results were obtained from the PERMANOVA analysis based on the Bray‐Curtis distances (nature of the sample: F‐value = 106.3; p‐value < 0.001; site: F‐value = 11.17; p‐value < 0.001; Figure S9).

Principal coordinate analysis (PCoA) plot showing environmental samples (triangles) and B. truncatus snail associated (dots) bacterial communities distributed along the two first PCoA axes computed based on pairwise Jaccard distances obtained between each pair of samples.

The Microbiota of

B. truncatus Snails Is Geographically Structured, Partly Driven by Their Genetic Background and the Geographic Distance and to a Lesser Extent by the Environment

3.4.2

The first two axes of the Principal Coordinates Analysis (PCoA) performed based on the Jaccard distance matrix obtained from B. truncatus microbiota explained 17.38% (9.85% and 7.53%) of the total variance (Figure 4). Our PERMANOVA analysis supports an effect of the sampling site with the structure of B. truncatus microbiota (F‐value = 6.09; p‐value < 0.001, Figure 4). The same results were observed using pairwise Bray‐Curtis distances (F‐value = 9.71; p‐value < 0.001, Figure S10). Conversely, we found no effect of the infection status of B. truncatus on the beta diversity of snails' microbiota, neither based on the Jaccard distances (Figure 4) nor based on the Bray‐Curtis distances (all p‐values > 0.05, Table S13 and Figure S10) even when accounting for the effect of the sampling site. Based on our series of partial Mantel tests conducted at the site level, we also found no effect of size difference between B. truncatus and beta diversity between snails' microbiota (Table S7).

Principal coordinate analysis (PCoA) plot showing the bacterial communities of B. truncatus according to their original sampling sites (see colour legend), and to their infection status by trematodes based on a molecular diagnostic (dot size) distributed along the two first PCoA axes based on pairwise Jaccard similarity indices computed between each sample. Each coloured small dot represents one snail host from a given sampling site. Large dots are snail hosts positive to trematode molecular diagnostic. Medium dots are the centroids of the bacterial communities of snails originating from each sampling site. Each line represents the connection between the sample and the centroids of the bacterial communities of snails originating from each sampling site.

Based on the multiple regression on distance matrices (MRM) and on our series of correlation tests (Figure 5), the variables that explain the mean Jaccard indices calculated between B. truncatus microbiota at the individual level are, in order of importance, the geographical distance between sites pseudo‐replicated at the individual level (r = 0.37; p‐value = 0.0001), the Nei's genetic distance between individuals (r = 0.17; p‐value = 0.0001), and the mean Jaccard indices calculated between the aquatic bacterial communities at the site level (r = 0.07; p‐value = 0.008) pseudo‐replicated at the individual level.

Distance‐decay scatterplots showing bacterial communities distance (determined by the pairwise Jaccard similarity index) as a function of geographic distance across sampling sites for environmental samples (A), and B. truncatus snails (B). Scatterplot of B. truncatus snails microbiota dissimilarity (determined by the Jaccard index) as a function of Nei's genetic distance between B. truncatus individuals computed based on 31 microsatellite loci (C). Blue lines of best fit are included to highlight trends.

Discussion

4

Composition of the Microbiota in Natural Populations of

B. truncatus

4.1

Bacterial Communities Display Lower Alpha‐Diversity in

B. truncatus Than in Environmental Samples

4.1.1

Our results provide new insights on the composition of the microbiota associated with B. truncatus natural populations. Overall, B. truncatus display microbiota less diversified than bacterial communities present in the surrounding freshwater environment. This was somewhat expected because bacterial communities collected from the environment encompass free‐living bacteria, bacteria that are released by all organisms associated with freshwater habitats and bacteria associated with freshwater microorganisms that could have been trapped during eDNA sampling (e.g., algae, invertebrates) (Mikhailov et al. 2019; Gendron et al. 2019; Sadeghi et al. 2021). Moreover, host microbiota are generally considered as a subsample of environmental bacterial communities resulting from a filter induced by the intrinsic properties of the hosts (Reese and Dunn 2018; Mazel et al. 2018; Suzzi et al. 2023). In particular, the host constitutes an ecological niche, which suitability for bacteria is largely influenced by hosts immune system that acts as a selective pressure, promoting the growth of some beneficial bacteria while inhibiting potential pathogens and maintaining the microbiota homeostasis (Belkaid and Harrison 2017). In line with this, seven among the 10 most represented bacterium phyla present in B. truncatus overall sites, are also found in the surrounding aquatic environment. The lower alpha‐diversity in snail microbiota compared to environmental bacterial communities was also described in recent studies conducted on other freshwater snails including Ampullaceana balthica and Galba truncatula (Herlemann et al. 2024; McCann et al. 2024). In our case, the observed difference in bacterial community richness between environmental samples and B. truncatus might have been sharpened here since eDNA was collected along the water columns down to the water–sediment interface hence providing a broad vision of environmental bacterial communities at the sampling sites.

Although global test indicated significant inter‐site differences, post hoc comparisons showed that only the Mbane site differed significantly from three other sites (Guia, Saneinte and Lampsar), with Mbane exhibiting lower alpha diversity than these sites. The particularly low alpha diversity observed in Mbane compared to Saneinte is intriguing, given that these two sites are separated by only 3.2 km and lie along the same shore of Guiers Lake. In our study, no obvious technical (e.g., difference in sequencing coverage), or measured biological factors (e.g., host size, host genetic diversity), can satisfactorily explain this pattern either at an individual or site scale. Ecological differences between these two sites, for instance in terms of food resources available for local B. truncatus populations, could be one explanatory factor explaining such a difference, although we did not collect this information in the field. Indeed, studies have suggested that, additionally to host physiological traits, dietary diversity (exogenous microbial diversity) also impact gut microbiota alpha diversity (Reese and Dunn 2018). Fine‐scale geographical studies at the scale of the lake accounting for more detailed ecological factors could be useful to better understand such results.

Bulinus truncatus Microbiota Composition in Natural Populations

4.1.2

We observed that Proteobacteria, Bacteroidota, and Firmicutes were the dominant phyla in the microbiota of B. truncatus across all sampling sites, which is consistent with previous studies on many different animals (McCann et al. 2024; Shi et al. 2024; Yuan et al. 2025). Among those phyla, the Bacteroidota and Firmicutes were, to some extent, enriched in snails compared to their surrounding aquatic environment. Those phyla are widely recognized for their roles in nutrient metabolism, fermentation, and host immunity (Huot et al. 2020; McCann et al. 2024; Yuan et al. 2025). The prevalence of Bacteroidota, in particular, may reflect host responses to physiological stress such as parasitic infection, as previously suggested for freshwater gastropods (McCann et al. 2024). Conversely, the low representation of Cyanobacteria and Actinobacteria, despite their abundance in surrounding habitats, highlights the selective pressures of the snail gut ecosystem (Herlemann et al. 2024). These patterns support the idea that microbiota composition reflects a combination of ecological exposure, host physiology, and historical association, echoing findings in B. straminea where microbial communities appear to shift in response to habitat changes (Yuan et al. 2025).

At the family level, the ‘common’ core microbiota of B. truncatus included Chitinophagaceae, Weekselaceae, and Flavobacteriaceae that all belong to the Bacteroidota phylum, as well as Oxalobacteraceae, Sphingomonadaceae, and Moraxellaceae that belong to the Proteobacteria phylum. Interestingly Proteobacteria, Bacteroidota and Firmicutes were previously detected as major bacterial taxa in a pool of B. truncatus initially originating from Spain and maintained in the laboratory for several generations (Huot et al. 2020). In particular, Flavobacteriaceae (Bacteroidota) and Sphingomonadaceae (Proteobacteria) are families that are shared by this laboratory‐reared Spanish B. truncatus strain and all natural B. truncatus populations from northern Senegal studied here. These bacteria are hence likely to display important functions for the development of B. truncatus .

Conversely, some bacterial families and in particular the Chitinophagaceae were found in high relative abundance in B. truncatus natural populations studied here but were absent or present at low abundance in the microbiota of the previously studied Spanish B. truncatus laboratory strain (Huot et al. 2020). Interestingly, the Chitinophagaceae family, as their name implies, can hydrolyze chitin from the environment and are found primarily in soils and aquatic sediments (Lim et al. 2009; Madhaiyan et al. 2015). Under natural conditions, B. truncatus feed on recalcitrant food resources including detritus, decaying macrophytes, diatoms and filamentous algae, and thus potentially require a highly diverse microbial communities for digestion (Madsen 1992; Van Horn et al. 2012). In this regard, a reduction in the complexity of the food resource available for the laboratory population (i.e., organic washed lettuce) may also explain a relaxation of selection pressure and a loss of some bacteria taxa involved in snails' digestion. Alternatively, although nonexclusively, the differences observed in the core microbiota between our B. truncatus field populations from Senegal and that of the Spanish laboratory strain may be a result of co‐evolutionary processes driven in part by the evolutionary history of B. truncatus populations; the population established in Spain having probably diverged from the populations of Senegal for a long time. This latest hypothesis implies that the genetic structure of B. truncatus populations resulting from their evolutionary history influences the bacterial communities associated with these populations. This is precisely what our results show at the scale of northern Senegal, as we discuss in the next section.

Structuring Factors of the Microbiota of Natural

B. truncatus Snails

4.2

Our result show that B. truncatus microbiotas are well structured geographically although they display a less pronounced ‘distance–decay’ pattern than that of aquatic bacterial communities. Communities of free‐living bacteria generally display such a pattern, although the robustness of the relationship may vary depending on the ecological context and the geographical scale (Clark et al. 2021). While the factors that lead to this ‘distance–decay’ pattern are still debated, it is accepted that it results at least from a limitation in the dispersion of bacterial communities and the heterogeneity of the environment (Green and Bohannan 2006; Sadeghi et al. 2021). We believe that, the strength of the ‘distance–decay’ pattern observed among aquatic bacterial communities in our study result from the fact that ecological heterogeneity and geographical distances between sites are tightly linked, with closer sites being associated to the same ecological system (e.g., lake de Guiers, River ‘Le Doué’). The less pronounced ‘distance‐decay’ pattern of B. truncatus microbiota observed here is in line with the recent study from Herlemann et al. conducted on the freshwater snail A. balthica in Northern Europe at a similar geographical scale (Herlemann et al. 2024). As the authors suggest in their model, we can here hypothesize that B. truncatus create more uniform environments for their associated bacteria, compared to the heterogeneity of the aquatic environments where B. truncatus and aquatic bacterial communities are found. Another non‐exclusive hypothesis would be that the structure of B. truncatus microbiota is determined by other predominant factors, and particularly some evolutionary factors specific to B. truncatus populations. In this respect, the MRM analysis shows that the genetic distance between B. truncatus hosts also explains the beta diversity of their associated microbiota (albeit slightly lesser) than the geographical distance. This suggests that the evolutionary history of B. truncatus populations, which is strongly influenced by migration and drift as shown by the strong isolation by distance (IBD) pattern observed, is also an important element in the structuring of B. truncatus microbiota. Thus, while phylosymbiosis between snail species and their microbiota is already documented at the inter‐specific level (Huot et al. 2020; Schols et al. 2023), we here argue that phylosymbiosis could also occur at the intra‐specific level.

The strong IBD pattern observed among B. truncatus populations was expected and fit well the results from previous genetic studies conducted on this species in different countries of Africa, including Senegal. Early studies based on allozymes or a limited number of microsatellites have shown that self‐fertilization is frequent in natural populations of B. truncatus (Njiokou et al. 1993; Jarne et al. 1994). On the other hand, populations of B. truncatus are often subjected to fluctuating episodes of extinction and recolonization that follow the alternation of wet and dry seasons that greatly shape the (sometimes temporary) habitats of this species. Episodes of extinction and recolonization by a small number of self‐fertilizing individuals generally explain the low intra‐site genetic diversity and genetic differentiation between populations, even if they are geographically close (Njiokou et al. 1993; Maes et al. 2022). Based on the 31 microsatellite markers newly developed in this study, we found no clear evidence of self‐fertilization, and the low values of genetic differentiation observed between geographically close populations suggest a non‐negligible gene flow between these populations that tends to fade as the geographical distance between populations increases. We explain these results by the fact that all sites included in this study are permanent habitats with water, and possibly B. truncatus populations, present all through the year. This limits (but does not fully exclude) temporary population extinctions and allows populations to develop demographically, thus reducing the effect of genetic drift and potentially encouraging cross‐fertilization. Moreover, water flow in the river Le Doué certainly promotes gene flow between geographically close populations and explains the little genetic differentiation observed between these populations. In this respect, it is interesting to note the higher genetic differentiation between the population established in the Guia canal, which is fed and connected to the Le Doué River, compared to populations established along the Le Doué River at equal or greater geographical distances. Similarly, the sampling sites on the banks of Lake Guiers are only 3.2 km apart, and it is likely that the sole effect of wind‐generated surface currents facilitates gene flow between these two populations. Conversely, high values of genetic differentiation were observed between geographically distant populations and are in line with the genetic patterns recently documented on B. truncatus populations in the same region based on 10,750 SNPs (Maes et al. 2022).

Combined together, the geographical distance between B. truncatus populations, the genetic structure of B. truncatus populations and to a lower extent the local environmental bacterial communities do not fully explain the structure of the microbiota of B. truncatus . This could be at least partly due to a sampling bias due to the unequal number of snails collected per site. Indeed, a smaller sample size can reduce the accuracy of allele frequency estimates and thereby increase the variance of Nei's genetic distances. However, this effect appears limited in our dataset, as illustrated by the smallest genetic distance observed between Ndiawara and Dioundou, two geographically close populations despite different sample sizes. Regarding the microbiota, a reduced number of snails per site could in theory underestimate alpha diversity and consequently inter‐site dissimilarities. Yet, no such systematic bias was detected in our data. We therefore believe that the impact of unequal sampling effort on our conclusions is negligible. Another explanation would be that other factors, probably ecological, which we missed in this study, are certainly at work in shaping the microbiota of B. truncatus . For example, low genetic differentiation and little geographical distance clearly fail to explain the huge difference observed between the microbiota of B. truncatus established at the two sites along Lake Guiers (Saneinte and Mbane). Moreover, since environmental bacterial communities are very similar at these two sites, this factor is also unlikely to explain the differentiation of the microbiota between these two B. truncatus populations. More generally, when controlling for B. truncatus population genetics and inter‐site geographic distances, the aquatic bacterial community structure explains little variation in snail microbiota structure. This was somewhat unexpected for at least two reasons. First, we might expect that some environmental bacteria are transferred to B. truncatus through several processes including ingestion during snail feeding or colonization of external snail tissues. Indeed, most freshwater snails are detritivorous or feed upon biofilms hence providing opportunity for different environmental bacteria to colonize and eventually establish within the gastrointestinal tract of snails established in different habitats (Madsen 1992; Van Horn et al. 2012; Kivistik et al. 2023). The little influence from aquatic bacterial communities on B. truncatus microbiota further supports the idea that B. truncatus acts as a filter for specific bacteria. This filtering effect, that might partly depend on the genetic background of B. truncatus lineages, limits the differences between the environmental bacterial communities that occasionally associate with B. truncatus . It is worth stressing that we here focused on B. truncatus whole microbiota and not specifically on the gut microbiota. By focusing only on these tissues, it is likely that we would have found a better congruence between environmental bacterial communities and B. truncatus gut microbiota. In fact, we did find several bacteria phyla that are present both in the environment and in B. truncatus as previously mentioned.

Finally, as part of the biotic environment, the presence of developing trematodes inside the snails did not appear to strongly influence the microbiota of B. truncatus . Most infected individuals harboured bacterial communities similar to those of non‐infected sympatric conspecifics. Because we did not determine infection stage (patent vs. pre‐patent), we cannot assess how long infected snails had been carrying trematodes. However, Portet et al. empirically showed that trematode experimental infections can induce transient dysbiosis in laboratory‐reared snails, followed by a gradual return to baseline within days or weeks (Portet et al. 2021). This suggests that some infected B. truncatus in our study may have been infected long enough for their microbiota to recover. Conversely, the few infected individuals whose microbiota differed from that of non‐infected snails may correspond to recent infections, before recovery could occur. Alternatively, these individuals with different microbiota may have been exposed to a trematode but were incompatible from an immune perspective, potentially triggering an immune response that led to a temporary dysbiosis (Schols et al. 2025). Because infection stage and compatibility were not assessed, these interpretations remain hypothetical.

Lastly, considering that we analysed whole snail microbiota, we potentially sequenced the microbiota of the B. truncatus snail and its potential infecting parasites so we could find different microbial signatures depending on the infecting trematode species as observed by Salloum and collaborators in the mud snail Zeacumantus subacarinatus infected by four different trematode species (Salloum et al. 2023). However, due to the low observed prevalence of infection in our study, we might have insufficient results to robustly explore the link between the infection status and the microbiota composition of B. truncatus , which means that we may have missed such signals.

In addition to these biological and methodological considerations, we also evaluated the potential impact of low‐level contamination in our sequencing data. A small number of ASVs detected in the negative controls also appeared in snail and environmental samples (Figure S4). Their relative abundances varied across samples and controls, with some ASVs being more abundant in controls and others more abundant in biological samples—a pattern consistent with low‐level background contamination commonly observed in amplicon‐based microbiome studies. Because these ASVs did not show a consistent contamination signature and their inclusion did not alter the major community patterns, they were retained in the dataset.

Conclusions

5

Our study highlights a well‐structured geographical pattern in the microbiota of natural populations of B. truncatus that is best explained by the geographic and genetic structure of B. truncatus populations. This suggests an important role of the evolutionary trajectory of B. truncatus populations in shaping their microbiota. This calls for further studies to investigate for potential genomic determinants that could influence the composition of hosts microbiota. Importantly however, some strong differences in the composition of B. truncatus microbiota from genetically and geographically close populations remain unexplained and suggest that some unexplored ecological factors or hosts' physiological traits (e.g., age or fitness) could also influence the microbiota of this species. Among these ecological factors we found that neither the bacterial communities from the surrounding aquatic environment nor the presence of trematode developing within hosts influence the overall structure of hosts microbiota. To identify potential environmental factors at play in shaping the microbiota of natural snail populations at the ecological scale, the additional temporal monitoring of the microbiota associated with hosts populations and that of the environment in which the hosts are established could provide a valuable approach. This would require a deep characterization of the environmental parameters associated with snails' aquatic ecosystems. In the case of B. truncatus , populations from habitats with high temporal environmental fluctuations such as temporary ponds should be targeted. More generally, such temporal monitoring studies in more temperate regions with important seasonal environmental changes would help better understand the potential environmental factors influencing the microbiota of freshwater gastropods. Finally, further research on the interactions between genetic and environmental variables would be necessary to specifically assess the plasticity of hosts microbiota, hence contributing to a better understanding of the relationship between the environment, hosts and their microbiota at the intra‐specific level.

Author Contributions

O.R., P.D., B.G., M.J.J. were responsible for research design. O.R., P.D., B.S., J.B. conducted the field work. M.J.J. and P.D. conducted the lab work. O.L. and E.C. developed the microsatellite dataset. M.J.J., P.D., O.R., E.T., and T.L. contributed to the data analysis. M.J.J., O.R. and P.D. wrote the initial draft of the manuscript and all authors contributed to revisions. O.R. and B.G. supervised the project.

Funding

This work was funded by the MICROVECT Project (défi clé RIVOC Occitanie Region, University of Montpellier); the European and Developing Countries Clinical Trials Partnership (EDCTP) program and by the French Agency for Food, Environmental and Occupational Health & Safety. This study is set within the framework of the “Laboratoire d'Excellence (LabEx)” TULIP (ANR‐10‐LABX‐41) and LabEx CeMEB (ANR‐10‐LABX‐04‐01). The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.

Disclosure

Benefit‐Sharing Statement: This project obtained from the Nagoya office in Senegal an exemption from authorization of access and use of genetic resources (number: 001339 of November 15, 2021, reference: V/L du 28 octobre 2021) by the Competent National Authority (Directorate of National Parks of Senegal).

Ethics Statement

The project has received approval from the National Ethical Committee (CNERS) of Senegal (agreement number: 00061/MSAS/CNERS/SP).

Conflicts of Interest

The authors declare no conflicts of interest.

Supporting information

Appendix S1: Microsatellite development. Appendix S2: Sequence‐Based Microsatellite Genotyping. Table S1: Abundance of B. truncatus , trematode prevalence among sites. Table S3: Paiwise genetic distances and F st values computed between sites from B. truncatus individual genotypes based on 31 newly developed microsatellite markers. Table S7: Partial Mantel test—Effect of morphological difference on beta diversity among the different sites. Table S8: Results of Kruskall‐Wall test comparing the B. truncatus alphadiversity depending on the sampling site of individuals based on ASV richness, diversity Shannon and Gini‐Simpson indices. Table S9: Results of pairwise post hoc comparisons between sampling sites based on ASV richness. Table S10: Results of pairwise post hoc comparisons between sampling sites based on Shannon diversity. Table S11: Results of pairwise post hoc comparisons between sampling sites based on Gini Simpson diversity. Table S12: Results of Kruskall‐Wall test comparing the B. truncatus alphadiversity depending on the infectious status of individuals based on ASV richness, diversity Shannon, and Gini‐Simpson indices. Table S13: Results of PERMANOVA analysis comparing the B. truncatus betadiversity based on Bray‐Curtis distance. Figure S1: Pictures of the nine sampled sites in Northern Senegal. Figure S2: Boxplot representation of B. truncatus snail sizes per sampling site. Figure S3: Graphical vision of the processed environmental samples. Figure S4: Common ASV between negative controls and B. truncatus samples (A) or environmental samples (B) and their abundance in each type of sample. Figure S5: Scatterplot displaying Nei's genetic distance between individuals on 31 loci as a function of the geographic distance across sampling sites. Figure S6: Alpha diversity computed on the bacterial communities from the environmental water‐sediment samples and isolated from B. truncatus snail hosts at each sampling site. Figure S7: Alpha diversity computed on the bacterial communities isolated from B. truncatus snail hosts at each sampling site. Figure S8: Venn diagram illustrating the number of shared phyla between B. truncatus snails and environmental samples. Figure S9: PCoA plot showing environmental samples and B. truncatus snails associated bacterial communities based on pairwise Bray‐Curtis distance obtained between each sample. Figure S10: PCoA plot showing the bacterial communities of B. truncatus snails according to their original sampling sites, and to their infectious based on pairwise Bray‐Curtis distance computed between each sample.

Table S2: Shell length and width of individual Bulinus snails across sampling sites.

Table S4: asv table before any cleaning step. See file.

Table S5: ASV found in the 4 PCR negative controls and their corresponding taxonomy.

Table S6: Read counts per sample during cleaning steps. See file.