Mathematical model

To specifically assess the effect of snail host abundance and the effect of possible trematode-trematode antagonistic interactions on the production and emission of cercariae of Schistosoma species we built a mathematical model. We constructed a dynamical model of an aquatic ecosystem in which a snail population is infected by n trematode species, one of them being the Schistosoma species. Trematodes enter the system as miracidia, which can infect susceptible snails. Infected snails shed cercariae, which are taken up by mammal hosts. This leads to a model with four types of dynamical variables: the density of miracidia for the different trematode species, the density of susceptible snails, the density of snails infected by one or more (up to a maximum of n) trematode species, and the density of cercariae for the different trematode species. The general model becomes rapidly overwhelming with increasing trematode diversity. In the Additional file 1 we show that the number of variables increases exponentially with the number of trematode species n. To keep the model manageable, we applied the following assumptions:

1. We required that infected snails have a smaller growth rate and a larger death rate than non-infected snails. For simplicity, we implemented an extreme version of this requirement, and assumed that the growth rate of infected snails and the death rate of non-infected snails are both zero.

2. We assumed that snails infected by multiple trematode species only shed cercariae of a single trematode species. We also required that the Schistosoma species is less competitive, and translated this difference in competitivity as a lower shedding rate of Schistosoma cercariae. For simplicity we implemented an extreme version of this requirement, and assumed that snails that are co-infected by the Schistosoma species and one or more other trematode species never shed cercariae of the Schistosoma species.

3. We assumed that all trematode species, including the Schistosoma species, are equivalent, in the sense that they have the same trait values apart from the shedding rate (see previous assumption). Specifically, the input rate of miracidia is the same for all trematode species, as are the miracidia loss rate, the death rate of infected snails, the cercariae loss rate and the uptake rate of cercariae by mammal hosts.

In the Additional file 1 we express these assumptions in terms of the model parameters, and use them to derive a drastic simplification of the general model. Specifically, we show that the model equilibrium can be computed from a system with a number of variables that increases only linearly in the number of trematode species n.

We studied the equilibrium properties as a function of two model parameters: the number of trematode species and the carrying capacity of the snail population. When varying the number of trematode species, we assumed that the miracidia input rate per trematode species is independent of diversity, so that the total input rate (i.e., summed over all trematode species) increases proportionally with diversity. In particular, the input rate of Schistosoma miracidia is unaffected by the presence of other trematode species. We also assumed that the total shedding rate (i.e., the shedding rate per infected snail and summed over all shedded trematode species) is independent of the number of species infecting the snail. This implies that the shedding rate per snail infected by a particular trematode species decreases with trematode diversity, as some of these snails are co-infected with other trematode species. This holds especially for snails infected by Schistosoma, as co-infection by another trematode species causes the Schistosoma shedding rate to drop to zero (see third assumption).

We reported the simulation outcome using three variables: the co-infestation rate, defined as the density of infected snails by two or more trematode species divided by the density of infected snails by one or more trematode species; the Schistosoma prevalence, defined as the density of snails infected by the Schistosoma species divided by the total density of snails (infected or not); and the Schistosoma cercariae ratio, defined as the density of Schistosoma cercariae divided by the total density of cercariae.

In the Additional file 1 we provide an in-depth description of the model equations and our simulations. We also provide the Matlab R2024a code (MathWorks, inc., Natick, USA) for all reported model simulations (Additional file 2 & 3).

Empirical field workStudy sites description and ethical consideration

An empirical fieldwork study was carried out in the region of the Senegal river basin during the dry season in February 2022 (Fig. 1). We targeted nine aquatic sites previously identified as potential S. haematobium transmission sites that differ from an ecological point of view (i.e. river, irrigation canal and lake) and hence possibly harbouring different trematode communities. Indeed, each of these transmission sites is close to a village where the prevalence and intensity of urogenital schistosomiasis among children was previously reported [19].

Fig. 1

Satellite map of the studied area in Northern Senegal. This map indicates the location of the nine targeted sites (black dots) and the name of the nearby nine villages associated to these Schistosoma haematobium transmission sites

The transmission sites close to the villages of Ndiawara (16°35′04″N / 14°50′58″W), Ouali Diala (16°35′56″N / 14°56′10″W), Dioundou (16°35′50″N / 14°53′22″W), Fonde Ass (16°36′20″N / 14°57′38″W) and Khodit (16°35′45″N / 14°56′42″W) are located in the middle valley along a tributary of the Senegal River (“le Doue” river). The transmission site close to the village of Guia (16°35′51″N / 14°55′31″W) is located along an irrigation canal that drains water from the river “le Doue”. The transmission sites nearby the villages of Mbane (16°16′15"N / 15°48′7"W) and Saneinte (16°14′32"N / 15°48′6"W) are located along the east shore of the lake de Guiers. The transmission site close to the village of Lampsar (16°6′34"N / 16°20′58"W) is in an inlet in the lower valley of the Senegal river delta (Fig. 1).

eDNA, snails and trematodes field sampling

To assess the abundance of B. truncatus and to identify all trematode species present at each of the nine studied sites as free-living stages, aquatic environmental DNA samples (eDNA) were collected using water filter capsules according to Douchet et al. (2022) [20]. Filtrations were carried out from the entire water column along the banks until the filter’s membrane clogged, and the volume filtered at each site was recorded. At each site but Khodit, 1.5 L of commercial spring water was filtered following the same protocol as a technical field negative control. Once the water filtrations completed, capsules were drained, filled with 50 ml of Longmire buffer solution to preserve eDNAs [21], vigorously shaken, and stored at room temperature and at dark until subsequent DNA extraction.

Once the eDNA sampling achieved, all freshwater snails found at each sampling site were systematically collected manually or by scooping the grass on the water bench using a colander for about 30 min to one hour. Snails were morphologically identified using a taxonomic key [22] and individualized on well plates filled with dechlorinated water and left to emit trematodes for two hours in the afternoon under natural sunlight. Cercariae from each emitting snail were collected and transferred onto FTA cards and stored at room temperature until DNA extraction and molecular identification. Coupled with eDNA samples, this approach allows a better characterisation of the trematode community at each transmission site and a quantification of the prevalence of each emitted trematode species among the locally established snail populations including B. truncatus. All snails morphologically identified as B. truncatus (either emitting or non-emitting) were preserved individually in alcohol until DNA extraction. Subsequent molecular analyses aimed at validating the species attribution of each B. truncatus individual, check for the presence and identification of potentially developing trematodes within each snail, refine measures of prevalence of each trematode species, and identify possible coinfections.

Experimental infection

Twenty-six non-emitting snails (although possibly naturally infected by trematodes) sampled at sites Guia, Fonde Ass, Khodit and Saneinte were kept alive and individually exposed to three S. haematobium miracidia for two hours in 24 well plates filled with dechlorinated water. The miracidia used for this experimental mollusc exposition were from a homogenised pool of S. haematobium eggs of eight children from the Lampsar village provided by the Sen_Hybrid_Invasion project (see acknowledgment section). The 104 exposed snails (i.e., 26 × 4) were maintained at 24 °C in dechlorinated water under natural light and fed ad libitum for a period of one month. During this period, survival was checked daily. At the end of this period, the infection status of each B. truncatus was monitored by emission of cercariae once a week until the snails died. To this end, B. truncatus were individualized on clean well plates filled with dechlorinated water and stimulated to emit trematodes for 4 h in the afternoon under artificial light. Cercariae released from emitting snails were collected and stored at −20 °C until DNA extraction for subsequent molecular identification. All along the experiment (i.e., 4 months), each emitting snail was kept individualised to identify possible trematode species emission shifts through time. At the end of the experiment, all B. truncatus individuals (either emitting or non-emitting) were preserved in alcohol individually until DNA extraction.

Molecular approachesDNA extractions from samples

To extract eDNAs from water filtrates, the Longmire solution contained in each capsule was poured into three 50 ml tubes as technical replicates. For the field negative controls (i.e., Spring water filtrates), each capsule content was recovered in one 50 ml tube only. All tubes were centrifugated at 16,000 × g for 20 min and the supernatant was removed. We then collected 0.25 g to 0.5 g of sediment from the pellet from each tube or 500 µl of Longmire remaining at the base of the tube when not enough material was observed. For negative controls, 500 µl of Longmire were systematically retained when discarding the supernatant. This pre-extraction step led to the processing of 35 samples (i.e., eight negative controls and three extraction replicates for each of the nine environmental samples). Total environmental genomic DNA was extracted from each triplicate and negative controls using the Qiagen’s Dneasy PowerSoil Pro Kit following supplier recommendation performing the physical lysis with a MagNA Lyser at a speed of 7000 × g for 30 s.

DNAs from all cercariae obtained from the field survey, from the exposure experiment and from all B. truncatus snails were extracted using the Qiagen DNeasy Blood & Tissue kit (Qiagen, Hilden, Germany). Regarding cercariae preserved onto FTA cards, we manually isolated 0.5 cm in diameter pieces containing the biological material using a sterile punch and placed them in 1.5 ml tubes. Four samples consisting in 0.5 cm in diameter pieces from cards on which ultrapure water had been deposited were also isolated as DNA extraction negative controls. For cercariae preserved in 1.5 ml tubes stored at −20 °C, we thawed the tubes and centrifuged them at 20,000 × g for 10 min before removing the supernatant. To obtain extraction and PCR negative controls, we also centrifuged two tubes that contained ultrapure water following the same protocol. For B. truncatus, we first rinsed the snails with tap water and ground them entirely and individually with a sterile pestle in 1.5 ml tubes. To obtain extraction and PCR negative controls, we also processed seven assuredly non-infected B. truncatus from our laboratory collections in the same way. Moreover, three extraction/PCR positive controls were prepared from three additional non-infected laboratory B. truncatus individuals, to which one S. haematobium cercariae from our laboratory collections was artificially added before DNA extraction processing. These pre-extraction steps led to the processing of 328 samples and 16 controls (i.e., 29 FTA card pieces and four associated negative controls, 28 cercaria preserved in alcohol and two associated negative controls and 281 B. truncatus including 173 B. truncatus from the field and 98 from the experimental infection and three and seven associated positive and negative controls). From these samples, we then followed the Tissue protocol as recommended by the supplier applying a two-hour lysis for cercariae and an overnight lysis for snails.

Taxonomical assignment of trematode cercariae by barcoding

To taxonomically assign each cercariae emitted from snail hosts to a trematode species, DNAs extracted from cercariae of each emitting snail were SANGER sequenced at the 28S D2 rDNA gene domain (28S) and at the 16S rDNA gene (16S) [20]. PCRs were run using the obtained 57 DNAs extracted from cercariae and the six negative controls on both markers. PCRs were performed using the GoTaq® G2 Hot Start Polymerase kit of Promega (Promega, France). Each PCR reaction contained Colorless Buffer at 1 × , MgCl2 at 1.5 mmol/L, dNTPs at 0.2 mmol/L, primers at 0.4 µmol/L, 1.25 units of GoTaq G2 Hot Start, 2 µl of DNA sample, and ultrapure water for a total PCR reaction volume of 35 µl. For the 28S marker, the PCR program was used as follows: An initial denaturation step at 94 °C for 3 min followed by 40 cycles with a denaturation step at 95 °C for 30 s, a hybridization step at 56 °C for 30 s and an elongation step at 72 °C for 30 s. We finally performed a final elongation step at 72 °C for 5 min. For the 16S marker, the same PCR program was used, except that 35 cycles were performed, the hybridization temperature was 54 °C and the elongation time was 15 s. Ten microliters of the resulting PCR products were migrated on a 2% agarose gel for 20 min at 135 V and revealed using a Vilber Infinity 1000 imaging system (Vilber, France). Each individual PCR product from the 57 cercariae DNA extracts that displayed an expected theoretical size was then sequenced in both the forward and reverse directions on an ABI 3730xl sequencer (Thermo Fisher scientific, USA) at the GenoScreen platform (GenoScreen, Lille, France). The sequences generated were aligned, trimmed, and taxonomically assigned with a MEGABLAST analysis [23] for the 28S and the 16S markers. For the taxonomical assignment of each cercariae, the common taxonomic rank to all hits above 96% of identity over 97% of coverage were kept.

Molecular taxonomic validation of B. truncatus and infection diagnostic

The 271 Bulinus spp. sampled in the field from which total DNA was extracted were diagnosed to verify that they belonged to the species B. truncatus using a molecular diagnosis by loop-mediated isothermal amplification (LAMP) based on the internal transcribed spacer (ITS2) B. truncatus species-specific [24]. DNAs extracted from Bulinus spp. were first diluted at 1/100e and a negative control consisting in a 1/100e diluted DNA extraction from Bulinus globosus was used. Each LAMP reaction contained an Isothermal Amplification Buffer II reaction buffer [20 mmol/L Tris-HCl, 10 mmol/L (NH4)2SO4, 150 mmol/L KCl, 2 mmol/L MgSO4, 0.1% Tween 20, pH 8.8 at 25 °C (New England Biolabs, UK)] at 1 X, additional MgSO4 at 3 mmol/L, dNTP at 1.0 mmol/L, internal primers FIP and BIP at 1.2 µmol/L, external primers F3 and B3 at 0.2 μmol/L, LOOP primers LB and LF at 0.4 µmol/L, 1 U of Bst 2.0 WarmStart DNA polymerase (New England Biolabs, UK), 1 µl of DNA sample, and ultrapure water for a total reaction volume of 10 µl as described in Blin et al., 2023 [24]. Reaction was performed in a thermocycler at 63 °C for 45 min followed by an enzyme inactivation phase at 80 °C for 5 min. Result visualizations were done using a final point visual detection of fluorescence after adding 1 μl of 1:50 diluted 10,000 × g concentration of SYBR Green (Invitrogen, USA) (green: positive = B. truncatus species; orange: negative = other than Bulinus species).

To detect trematodes within the tissues of each B. truncatus, we used the 16S rDNA gene metabarcode initially developed to characterize trematode communities from eDNA [20] as a molecular diagnostic tool. PCRs were performed using the GoTaq® G2 Hot Start Polymerase kit of Promega (Promega, France) following the same PCR condition as described previously, except that reactions were performed in a final volume of 10 µl and using 2 µl of DNA extracted from Bulinus spp. diluted 1:100. Amplification success was assessed visually by migrating the ten microliters of the resulting PCR products on a 2% agarose gel for 20′ at 135 V and revealed using a Vilber Infinity 1000 imaging system (Vilber, France).

Trematodes-specific MiSeq sequencing on the infected B. truncatus individuals and on the eDNA to characterize trematode communities

A total of 137 positive 16S metabarcoding NGS libraries were prepared following the Illumina two-step PCR protocol, using the Trem_16S_F1 and the Trem_16S_R2 primers set up with Illumina adapters (Illumina, USA) using 2 µl of eDNA or DNA extracted from B. truncatus diluted at 1:100e as previously described [20]. Libraries were sent and paired end sequenced (2 × 250 bp) on an Illumina MiSeq™ at the BioEnvironnement platform (University of Perpignan Via Domitia, France). Three samples were eliminated from the analyses of the MiSeq sequencing because they did not meet the defined threshold of 25,000 sequences required to normalise samples by rarefaction (i.e., one non emitting B. truncatus from the site Ouali Diala, and two PCR duplicates from water eDNA).

Abundance of B. truncatus determination by ddPCR on eDNA

Abundance of B. truncatus at each study site was assessed by digital droplet PCR (ddPCR) from eDNA sampled by water filtration according to [25]. We first pooled the eDNA extraction triplicates from each of the nine sampling sites. As negative control, we also pooled DNAs extracts from one individual of each mollusc species sampled during the field work in addition to DNAs extracts from three B. truncatus sister species (i.e., B. globosus, B. senegalensis and B. umbilicatus) from our laboratory collection. As positive control, we used a DNA extract at 0.01 ng/µl from one B. truncatus from our laboratory collection. We ran ddPCRs using the TaqMan technology on a QX200 AutoDG Droplet Digital System (BioRad, USA) at the Bio-environnement platform (Perpignan, France). We used the Btco2F (5′-ATTTTGACTTTTACCACCAT-3′) and Btco2R (5′-GATATCCCAGCTAAATGAAG-3′) primers combined with the FAM-labelled probe Btco2P (5′-TCGAAGGAGGGGTTGGAACAGG-FAM-3′). ddPCR reactions were performed using the BioRad ddPCR Supermix for Probes (No dUTP) (BioRad, USA). Each ddPCR reaction contained the BioRad MIX at 1 X, primers at 0.5 µmol/L, probes at 0.25 µmol/L, 7 µl of eDNA template or 2 µl of control DNAs templates, and ultrapure water for a total ddPCR reaction volume of 20 µl. After the droplets generation, the ddPCR program was used as follows: An enzyme activation step at 95 °C for 10 min followed by 40 cycles with a denaturation step at 94 °C for 30 s and an annealing/extension step at 60 °C for 60 s, with a ramp setting to 2 °C/s. We finally performed an enzyme deactivation step at 98 °C for 10′. The resulting fluorescence signal were analysed using QuantaSoft software V1.7 (BioRad, USA). A signal was considered positive (i.e., with the presence of B. truncatus DNA) if at least one positive droplet was detected and if the positive droplet displayed the same order of fluorescence magnitude as the positive droplet obtained from a ddPCR positive control.

Data analysisCharacterisation of trematode communities present in the water and exploiting B. truncatus populations

The resulting amplicon sequence datasets from the MiSeq sequencing was processed using the Find Rapidly OTUs with Galaxy Solution (FROGS) [26] according to reference [20]. Briefly, the produced datasets were pre-processed by filtering out the sequences to keep amplicon sizes from 150 to 400 nucleotides. The remaining sequences were next clustered into operational taxonomic units (OTUs) using the swarm algorithm and using denoising and an aggregation distance of three [27]. The resulting dataset was filtered out for chimeras using VSEARCH [28]. Singletons and underrepresented clusters (i.e., clusters whose number of sequences were < 0.1% of the total number of sequences) were removed. Lastly, we conservatively considered that a given OTU was present in a library if its number of sequences was > 0.1% of the total number of sequences in this library.

Each OTU was next assigned to a taxonomic level (either a species or a genus) using a two-step BLAST affiliation process. The first BLAST analysis was computed using the standalone blastn program (NCBI, USA) contained in the BLAST + package and a custom trematode sequence database containing 174 sequences (82 sequences from the NCBI database and 92 sequences from a custom internal trematode sequence database including the sequences obtained from the amplicons generated by the SANGER sequencing on cercariae). The second BLAST analysis was performed using the online MEGABLAST tool and based on the non-redundant database without restricting parameters to achieve affiliation of OTUs that could not be assigned in the first BLAST analysis. OTUs were assigned to a species if the sequences presented a minimal blast coverage of 97% and a pairwise identity above 99% with the affiliated sequence. OTUs were assigned to a genus if the sequences presented a minimal blast coverage of 97% and a pairwise identity above 96% with the affiliated sequence. OTUs that could not be assigned to a species or genus were assigned to a higher taxonomic rank using a clustering method based on the pairwise genetic distances between the OTUs and the same set of 174 sequences as above. We then aligned these sequences with T-Coffee on EMBL-EBI [29] and built a neighbor joining phenetic tree based on the percentage of nucleotide differences from the obtained alignment using Jalview version 2.11.1.4 [30] for visualization.

Subsequent analyses were then performed on R version 4.3.1 (Lucent Technologies, Jasmine Mountain, USA). The significance levels of statistical tests were set at 0.05. Each library was normalized by rarefaction at 25,000 reads using the package Vegan version 2.6–4. The trematode composition, species number and relative abundance in each B. truncatus individual or water sample was next assessed using the package phyloseq [31]. We assessed the co-infection rates among infected B. truncatus and the number of trematode species involved in each co-infection. We also tested whether parasite aggregation at the B. truncatus individual scale occurred by comparing the observed distribution of trematode species within hosts with a theoretical distribution following a negative binomial distribution using a Chi-square test.

Abundance of B. truncatus at each sampling site

The B. truncatus eDNA copy number per liter of filtered water (\(_\)) estimated from the ddPCR results was used as a proxy of the abundance of B. truncatus at each sampling site. \(_\) was calculated from the following equation [25]:

$$C_ = \frac *V_ }} }}}} }}$$

With \(_\): the number of eDNA c/L for the amplified sample; \(_\): the copy number per reaction volume; \(_\): the total volume of eluted DNA after extraction; \(_\): the volume of extracted DNA used for ddPCR reaction; and \(_\): the total volume of filtered water.

Assessing the link between the relative abundance of Schistosoma quantified in co-infected B. truncatus and the emission status of Schistosoma cercariae by co-infected B. truncatus

To determine whether the dominance status of Schistosoma species in terms of percentage of reads of these species within B. truncatus individuals co-infected with other trematode species can be used as a proxy of Schistosoma cercarial release, we performed a generalised linear model (GLM) using the package stats version 4.3.1 implemented in R. A subset from our dataset that contained 14 B. truncatus individuals from the field work and from the experimental infection were used for this analysis. These individuals were all emitters, co-infected and harboured at least S. haematobium or S. bovis. We set the emission status (i.e. 0 and 1: non-emitting and emitting Schistosoma sp.) as the dependent variable and the relative abundance of the respective Schistosoma species in terms of percentage of sequences compared to the other co-infecting trematode species as explanatory variable. The model was built assuming a quasibinomial distribution of the data.

Effect of trematode species richness, co-infection rates and host abundance on the average dominance of schistosomes

To test our initial hypotheses, we ran three independent generalized linear mixed models (GLMM) using the package lme4 version 1.1–33 implemented in R to infer the effect of (i) the overall trematode species richness using B. truncatus as an intermediate host at the population level, (ii) the co-infection rate among the infected B. truncatus snails, and (iii) the abundance of B. truncatus; on the average dominance of schistosomes within host populations (i.e., the relative abundance of S. haematobium or S. bovis in terms of percentage of sequences compared to the other co-infecting trematode species). Only naturally infected B. truncatus from the field were considered in this analysis, excluding B. truncatus from the experimental infection. To account for pseudo-replication, we set the sampling site as a random factor. For these analyses, the co-infection rate variable was transformed into two categories (i.e. < 50% and > 50% co-infection among infected individuals). The B. truncatus abundance variable was also transformed into two categories (i.e. < 200 and > 200 \(_\)) to compare sites at low and medium to high density according to Mulero et al. 2020. The variable “trematode species richness” (that uses B. truncatus) was coded as the number of trematode species other than Schistosoma which use local B. truncatus population at the site level. We assumed a binomial error term in the three models tested.