Trematode infection affects shell shape and size in Bulinus tropicus.

Trematodes can increase intraspecific variation in the phenotype of their intermediate snail host. However, the extent of such phenotypic changes remains unclear. We investigated the influence of trematode infection on the shell morphology of Bulinus tropicus, a common host of medically important trematodes. We focused on a snail population from crater lake Kasenda (Uganda). We sampled a single homogeneous littoral habitat to minimize the influence of environmental variation on shell phenotype, and barcoded snails to document snail genotypic variation. Among the 257 adult snails analysed, 99 tested positive for trematode infection using rapid-diagnostic PCRs. Subsequently we used high-throughput amplicon sequencing to identify the trematode (co-)infections. For 86 out of the 99 positive samples trematode species delineation could discriminate among combinations of (co-)infection by 11 trematode species. To avoid confounding effects, we focused on four prevalent trematode species. We performed landmark-based geometric morphometrics to characterize shell phenotype and used regressions to examine whether shell size and shape were affected by trematode infection and the developmental stage of infection (as inferred from read counts). Snails infected by Petasiger sp. 5, Echinoparyphium sp. or Austrodiplostomum sp. 2 had larger shells than uninfected snails or than those infected by Plagiorchiida sp. Moreover, the shell shape of snails infected solely by Petasiger sp. 5 differed significantly from that of uninfected snails and snails infected with other trematodes, except from Austrodiplostomum sp. 2. Shape changes included a more protuberant apex, an inward-folded outer apertural lip and a more adapically positioned umbilicus. Size differences were more pronounced in snails with 'late' infections (>25 days) compared to earlier-stage infections. No phenotypic differences were found between snails infected by a single trematode species and those harbouring co-infections. Further work is required to assess the complex causal links between trematode infections and shell morphological alterations of snail hosts.

Introduction

Trematodes are parasitic flatworms with complex life cycles that include a snail as first intermediate host and a variety of vertebrates as final hosts (Esch et al., 2001). The heavy medical and veterinary burden caused by schistosomiasis (bilharzia) and fasciolosis after infection with certain Schistosoma and Fasciola species, respectively, has generated ample scientific interest (Hotez et al., 2014). Trematodes are intimately linked to their snail hosts through coevolutionary dynamics (Lockyer et al., 2004) and because several trematode species can infect the same snail species, or even individual, trematode communities can develop within a snail population (Lafferty et al., 1994). Inside the snail, trematodes occupy the mantle, the digestive and/or reproductive organs and reproduce asexually. They can disrupt the endocrine system of the host and may cause serious damage either by diverting nutrients or by feeding directly on host tissue (Lafferty and Kuris, 2009). Depending on the specific interaction between host and parasite, infections can increase snail mortality and/or decrease fecundity, up to castration of the host (Sorensen and Minchella, 2001). Upon infection, parasite fitness partly depends on the phenotype of the snail host, which determines immunological responses and the amount of nutrients and space available for parasite survival and reproduction. Notably, shell phenotype may be a key determinant of trematode fitness due to its protective role for the snail and its size and shape limits the space and, therewith, nutrient availability for the parasite. Therefore, some trematode species have been suspected to manipulate the shell phenotype of their host to improve their reproductive potential (e.g., McCarthy et al., 2004; Levri et al., 2005). However, changes in snail phenotype upon trematode infection could also result from histological damage or endocrine disruption linked to the pathology caused by the parasite, rather than from parasite manipulation (Hay et al., 2005). The range of effects may even be more complex under co-infection by multiple parasite species, as competition among parasites and various strategies to exploit host resources could alter snail phenotype in various directions. Given this complexity, the impact of (co-)infection on snail phenotype is currently understudied (but see Miura and Chiba, 2007 and Lagrue et al., 2007). Observed phenotypic changes linked to trematode infections include decreases or increases in shell size (Raymond and Probert, 1993; Krist and Lively, 1998; Zakikhani and Rau, 1999; Probst and Kube, 1999; McCarthy et al., 2004), alterations of shell shape or ornamentation (Krist, 2000; Hay et al., 2005; Żbikowska and Żbikowski, 2005; Lagrue et al., 2007) and behavioural changes (Miura et al., 2006). Independent of parasite infection, intraspecific phenotypic variation is omnipresent in natural populations due to genetic differences among individuals and diverse environmental influences (Bolnick et al., 2011). In snails, intraspecific variation in shell size, shape and ornamentation are linked to a variety of biotic and abiotic environmental factors, e.g., presence of predators, wave velocity and water chemistry (Krist, 2002; Trussel, 1997). Overall, understanding the causal relationship between snail shell phenotype and trematode infection and physiological feedback mechanisms is a challenging task that requires experimental studies. However, detailed characterization of the intraspecific variability of snail shell morphology in a documented natural setting with limited environmental variation is also required to understand the complexity of influencing factors, and how they interact. Characterizing the interactions between trematode parasites and their snail host is important to understand the ecology and co-evolutionary dynamics of parasite transmission.

Despite the widespread distribution of Bulinus throughout Africa, southern Europe, and the Middle East (Brown, 1994) and the diverse trematode communities that are hosted by Bulinus species (Chingwena et al., 2002; Loker et al., 1981; Mohammed et al., 2016), the impact of trematode infections on Bulinus shell phenotype remains poorly understood. Indeed, only laboratory analyses investigating morphological changes induced by Schistosoma margrebowei on Bulinus natalensis (Raymond and Probert, 1993) and by Schistosoma haematobium on Bulinus truncatus, Bulinus senegalensis and Bulinus globosus (Fryer et al., 1990) have been published. Those studies showed that shell size may increase or decrease under Schistosoma infections, depending on the host-parasite system analysed (e.g., infection correlated with larger shells in the B. senegalensis - S. haemotobium association but with smaller shells for B. globosus - S. haemotobium). However, laboratory experiments have limitations as they usually focus on infections by a single trematode species in a highly simplified experimental environment. Perhaps such experiments do not reflect natural conditions and dynamics well, as greater ecological variability and co-infections are common in nature. Thus, detailed studies of trematode infections on snail phenotypic variation in natural snail populations are needed.

Against the abovementioned complexity, we here study the role of natural trematode infections on the phenotype of the intermediate snail host, for which three main elements are required. First, we need to be able to accurately characterize snails and their infecting trematodes. Second, we require an accurate documentation of shell shape and size. Finally, we need to effectively characterize both of the above in a natural system with sufficient environmental control for hypothesis testing. Newly developed techniques combining the screening of snail populations for infections using diagnostic PCRs (Schols et al., 2019) with high-throughput amplicon sequencing (HTAS, Hammoud et al., 2022) are promising in that they allow sensitive detection of trematode (co-)infections and accurate identification provided closely matching reference sequences are available in public DNA databases (Schols et al., 2020). Studies of the effect of trematode infections on snail shell morphology currently rely mainly on traditional measurements of shell morphological features (e.g., aperture length or width) or categorical indices of ornamentation (Levri et al., 2005; Lagrue et al., 2007). Methods that preserve geometric relationships, e.g., landmark-based morphometrics, allow a more holistic evaluation of shell shape variation, however (Van Bocxlaer and Schultheiß, 2010; Gustafson and Bolek, 2016). Finally, studies addressing the impact of trematode infections on snail shell morphology in natural environments often assume homogeneous snail populations and environments (e.g., Hay et al., 2005; Levri et al., 2005; Lagrue et al., 2007), but variation in environment and host genotype also influence snail morphology (Boulding and Hay, 1993; Gustafson and Bolek, 2016; Westram et al., 2018). In light of these confounding factors, we propose that the ideal site to study how parasite infections affect snail morphology is a geographically-restricted and isolated system inhabited by a single population of Bulinus with a high prevalence of trematode infections. After surveying multiple crater lakes located in western Uganda in which Bulinus populations had been reported (Tumwebaze et al., 2019), we retained Lake Kasenda as a promising study system.

In this study, we examine the impact of trematode infections on shell morphology in a wild population of Bulinus tropicus by combining recent methodological advances that allow simultaneous genotyping of snails and their infecting trematodes with detailed phenotypic characterization of snail hosts using landmark-based geometric morphometrics. More specifically, we test: 1) whether trematode infections affect shell size and/or shape in B. tropicus; 2) whether the observed effects are trematode species-specific; and 3) whether the effect of co-infection differs from that of infection by a single trematode species?

Material and methods

Study area and specimen collection

We collected the Bulinus snails for this study on February 20th, 2019, in Lake Kasenda of the Ndali-Kasenda cluster of volcanic crater lakes located in western Uganda between the Rwenzori mountains and Kibale forest. Characterization of the Bulinus fauna in western Ugandan crater lakes by Tumwebaze et al. (2019) suggested that Bulinus tropicus is the only Bulinus species that occurs in Lake Kasenda. Our sampling was restricted to a 20-m stretch of forested shoreline on the north-eastern side of the lake (N 0° 25′ 58.8612″, E 30° 17′ 29.1474″), with an environmentally uniform aquatic habitat of gently sloping muddy lake bottom covered with live and dead aquatic macrophytes. Bulinus specimens were collected in 45 min by two people (CH, JT), through visual inspection of the submerged vegetation and scooping of vegetation and sediments using a metal sieve mounted on a handle (down to 0.5 m water depth). Snails were sacrificed by heat shock, immediately preserved in 80% ethanol, and later stored in the Royal Museum for Central Africa (RMCA, Tervuren, Belgium) until further processing.

Genetic characterisation of Bulinus tropicus

To examine snail population structure, we genetically characterized all snails in our study. Snails were first visually inspected and measured using a ruler, to exclude specimens that are potentially juvenile (<7 mm) or with damaged shells from further analyses. All soft tissue was removed from each snail and the shell was kept for subsequent phenotypic characterization. The tissues were placed on absorbing paper to remove traces of ethanol and subsequently they were homogenized using a sterilized scalpel. DNA was extracted using the E. Z.N.A. Mollusc DNA Kit (OMEGA Biotek, Norcross, GA, USA) and final extracts were diluted in ultrapure water in a 1:10 concentration. Shells and DNA extracts are deposited in the RMCA (specimen codes BE_RMCA_MOL_DNA.000548 to BE_RMCA_MOL_DNA.000966). We barcoded each snail with a fragment of the gene cytochrome c oxidase subunit 1 (COX1) using primers of Folmer et al. (1994). We performed Polymerase Chain Reaction (PCR) using the Qiagen™ Taq DNA polymerase kit containing 1.5 mM PCR buffer (Qiagen™), 0.6 mM dNTP mix (Qiagen™), 1.5 mM MgCl2, 0.45 units of Taq Polymerase (Qiagen™), 0.8 mM μM HCO primer and 0.8 μM LCO primer. PCR was performed in a Tprofessional Thermocycler (Biometra™) with initial denaturation at 95 °C for 5 min, followed by 39 cycles of 95 °C for 45 s, 50 °C for 45 s and 72 °C for 45 s and a final elongation step at 72 °C for 10 min. PCR products were subjected to gel electrophoresis on a 1.2% agarose gel with Midori Green Direct® staining at 120 V for 30 min, followed by visualization under UV light. All PCR products were purified using the ExoSAP (Fermentas™) kit and sequenced by Macrogen™. For each sample, the HCO and LCO sequences were assembled into a consensus sequence in Geneious v.2020.2.4. The resulting chromatograms were checked for ambiguities to optimize sequence quality. All COX1 sequences were then aligned in Geneious using MUSCLE (Edgar, 2004). Aligned DNA sequences were translated to amino acids to avoid inclusion of nuclear mitochondrial DNA (NUMT, Lopez et al., 1994). We assessed genetic diversity by computing the genetic distances among all specimens and illustrated population genetic structure with a TCS haplotype network (Clement et al., 2002) using PopART v. 1.7 (Leigh and Bryant, 2015).

Characterization of trematode infections

We characterized trematode infections in B. tropicus specimens by performing a Rapid Diagnostic PCR (RD-PCR; Schols et al., 2019) on the whole-body DNA extracts of all snail specimens. Then, we subjected all infected snail specimens to the HTAS workflow to identify their infecting trematode species (Hammoud et al., 2022). The procedure and conditions of RD-PCR reactions are provided in supplementary Text 1. As the studied snails were barcoded (see above), we focused our HTAS characterization on trematodes. First, we amplified five trematode markers (beyond four snail markers) in a single multiplex PCR for each snail extract: two fragments of COX1, one fragment of NADH dehydrogenase subunit 1 (NAD1), one fragment of cytochrome b (cytb) and one fragment of the internal transcribed spacer 2 of the rRNA cluster (ITS2). The PCR products were tagged in a second PCR reaction using sample-specific identifiers, pooled, and sequenced using an Illumina MiSeq v3 sequencer (300 bp pair-ended). Subsequently, the sequenced reads were trimmed, filtered and denoised followed by inference of trematode Amplicon Sequence Variants (ASVs) and chimera removal with the dada2 pipeline (Callahan et al., 2016). The trematode ASVs profiles of each snail specimen were inspected to identify cross-contaminating sequences and to assess the overall level of cross-contamination in the sequencing run after which they were filtered based on read support to remove any potentially unidentified cross-contamination. In addition to the infected Bulinus tropicus specimens, our sequencing run also included multiple PCR1 duplicates of 8 laboratory-reared Biomphalaria glabrata, either uninfected or experimentally infected with Schistosoma mansoni (2 different strains) or Schistosoma rodhaini. Those B. glabrata specimens were used as negative and positive controls in lab work and sequencing and to quantify cross-contamination (Supplementary Table 1). The quality-controlled ASV profile obtained for each sample was then used to delineate phylogenetic species of trematodes per marker using the Assemble Species by Automatic Partitioning (ASAP) method (Puillandre et al., 2021). For each marker, our unique trematode ASVs were aligned with overlapping, identified GenBank sequences using MUSCLE in Geneious. The ASAP procedure requires a range of value of genetic distances within which the threshold for species delimitation is expected. We decided to provide a broad range (Supplementary Table 2) as the primers used for HTAS were designed to amplify highly variable regions of the targeted markers (Hammoud et al., 2022). Then, we ran the ASAP algorithm on the alignments of each marker and picked the partitioning scheme with optimal ASAP score within the range set (Puillandre et al., 2021). To evaluate the performance of the ASAP procedure, we assessed, for each marker, whether the GenBank sequences included in the alignments were correctly partitioned into their recognized species entities. More specifically, we calculated the occurrence of incorrect species delimitation for GenBank sequences (i.e., whether conspecific specimens were split across partitions or specimens belonging to different species were pooled together). Finally, to characterize trematode species across markers, we linked ASAP-inferred ASV clusters for our five markers based on their co-occurrence in snail specimens. To identify the infecting trematodes to the lowest possible taxonomic level, we analysed the resulting ASVs clusters per marker and inferred trematode species by BLASTing against NCBI GenBank (Camacho et al., 2009). If trematode infection causes changes in snail shell morphology, we would expect differences to accumulate over time. Just after infection morphological changes would be minimal, as hardly any shell growth would have taken place, whereas effects would accumulate for long-lasting infections. Trematode read proportions (i.e., the number of reads allocated to ASVs of trematode markers divided by the total number of quality-controlled reads for a sample) semi-quantitively reflect the increase of parasite biomass during the development of infections in their hosts (Hammoud et al., 2022). Therefore, we used the proportion of reads supporting trematode amplicon sequence variants (ASVs) in a sample to infer the developmental stage of infections, which we characterized only for single-species infections by our prevalent trematode species. Based on the characterization performed by Hammoud et al. (2022), we classified samples with a trematode read proportion ≥45% as ‘late’ stage (corresponding to >25 days after infection), and the others as ‘early’ stage.

Characterization of shell phenotype

Shells were cleaned with a needle in diluted bleach (5%) and photographed on a macrophotography stand in apertural view with the columellar axis horizontal and the aperture as such that a vertical tangent could be fit to its most peripheral margin (Van Bocxlaer and Schultheiß, 2010). Several high-resolution pictures per specimen taken with a Canon EOS 1200D camera were stacked with Zerene® software (T2019-10-07-1410). The resulting pictures were post-processed (cropped and rotated) using the GNU Image Manipulation Program v. 2.10.22 (GIMP Development The GIMP Development Team, 2019) and grouped into a single TPS-file using tpsUtil v. 1.78 (Rohlf, 2019). We scaled images and digitized 10 landmarks and 4 curves of semi-landmarks, which contained 15, 25, 15, and 7 equally spaced semi-landmarks, respectively (Fig. 1), on each individual in tpsDig v. 2.31 (Rohlf, 2017). Subsequently, all specimens were aligned by generalized Procrustes superimposition in CoordGen8 of the Integrated Morphometrics Package (Sheets, 2014), after which semi-landmarks were slid along curves in SemiLand8 (Sheets, 2014) to minimize Procrustes distance. During this procedure 8, 13, 8 and 4 semi-landmarks from the curves I, II, III and IV, respectively, were defined as helper points, resulting in a final dataset of 10 landmarks, 29 semi-landmarks and centroid size (CS).

Fig. 1

Illustration of a Bulinus shell with indication of our ten landmark points and four semi-landmark curves, numbered with roman numbers I to IV and indicating the original number of equidistant semi-landmarks between brackets. Shell height = 8 mm.

Analyses of the impact of trematode (co-)infections on shell morphology

Exploratory analyses

First, we examined whether the probability of trematode infection differs among snail haplotypes with a G-test of independence (DescTools package v. 0.99.44; Signorell, 2021), i.e., does the proportion of infected snails within each haplotype follow random expectations. Whereas this test considers general trematode infections (i.e., all infections collectively), we likewise tested, with a separate G-test, whether trematode-specific infections (i.e., per trematode species) differ among snail haplotypes. For these tests and in subsequent analyses we removed all snail specimens that did not have barcoding or HTAS results. We also removed snail haplotypes that were represented by less than five specimens. Second, we used a Kruskal–Wallis test to examine whether shell size, i.e., CS, differs among haplotypes. The distribution of CS by haplotypes was visualized using boxplots. All statistical analyses were performed in R v. 4.1.1 (R Core Team, 2019). Finally, we tested whether shell shape differed among haplotypes with a permutational multivariate analysis of variance (MANOVA) of the Procrustes shape coordinates with haplotypes as grouping factor, using the adonis function of the vegan package v. 2.5–7 (Oksanen et al., 2017) and 1000 iterations of residual randomization. For all statistical analyses the significance level was fixed at α = 0.05.

Influence of general trematode infections on shell phenotype

We investigated whether general trematode infections influence shell size of B. tropicus while controlling for the potential influence of snail genotype by building a linear regression model including both general infection status and haplotype as predictors of CS. Then, we used a Procrustes regression as implemented in the R package geomorph v. 4.0.1 (Adams et al., 2021) to investigate the effect of general trematode infections, allometry (i.e., phenotypic change with size) and snail genotype on shell shape. The regression model was built using the Procrustes shape coordinates as response variable, and snail haplotype, general infection status, CS and interaction between CS and infection status as explanatory variables. The interaction term was included to test whether allometric shape variation differs between infected and uninfected specimens. The Procrustes regression was fitted using 1000 iterations of residual randomization.

Differential influence of infecting trematode species on shell phenotype

To further investigate whether the effect of infection on shell size differed among trematode species, we built a linear regression model using CS as response variable and snail haplotype and trematode-specific infection status as predictors (Table 2 model 3). To avoid potential confounding effects related to co-infection, we only included cases where infection occurred by a single trematode species. Furthermore, we only focused on trematode species that infected at least ten snail specimens to ensure sufficient statistical power. We investigated whether the effect of infections on shell shape was trematode species-dependent using a Procrustes regression with the Procrustes shape coordinates as response variable, and haplotype, species-specific infection status, CS and interaction between infection status and CS as response variables (Table 2: model 4). As the interaction term indicated that allometric shape variation differed among the 5 levels of species-specific infection status (i.e., uninfected, or infected by either one of the 4 prevalent trematode species), we also performed pairwise comparisons of the allometric slopes for each of those infections (Table 2: test 5) by using the pairwise function of RRPP v. 1.1.2 (Collyer and Adams, 2021). The Procrustes regression and pairwise comparisons of allometric slopes were fitted using 1000 iterations of residual randomization. Finally, to document differences in shell shape among infected vs. uninfected snails, or even among snails infected with various trematode species, we reconstructed the morphospace occupation of the studied B. tropicus population. We performed non-metric multidimensional scaling (NMDS) in two dimensions on the landmark and semi-landmark dataset using functions of MASS v. 7.3–54 (Venables and Ripley, 2002) and vegan v. 2.5–7. As the stress value of this analysis was high, we investigated NMDS robustness by adding a third dimension. Finally, we analysed whether shell size differed among uninfected snails versus those with ‘early’ and ‘late’ infections for each of the four prevalent trematode species using Wilcoxon tests (Table 2: test 6) and boxplots. Likewise, we used Procrustes regressions and NMDS biplots to examine differences in shell shape among uninfected snails and those with ‘early’ and ‘late’ infections (Table 2: models 5 to 10).

Table 2

Summary of the statistical analyses performed and their output. Predictors with p-values below the significance level of 0.05 are marked in bold. N = sample size for the analysis, Df = degree of freedom. Regression coefficients are provided with standard error.

Analysis	Type	Response	Predictor	Test statistic	N	Df	Coefficient	p-value	R2
Test 1	G-test	Infection (general)	Haplotype	0.48 (W)	198	1	-	0.486	-
Test 2	G-test	Infection (species-specific)	Haplotype	22.10 (W)	198	17	-	0.181	-
Test 3	Kruskal–Wallis analysis of variance	CS	Haplotype	1.51 (K–W chi-squared)	198	1	-	0.219	-
Test 4	Permutational MANOVA	Procrustes shape coordinates	Haplotype	0.98 (F)	198	1	-	0.404	0.005
Model 1	Linear regression	CS	Intercept		198	1	2.86 ± 0.04	<0.001	0.049
			Infection (general)	10.95 (F)		1	0.14 ± 0.05	0.001
			Haplotype: H2	1.28 (F)		1	−0.05 ± 0.04	0.259
Model 2	Procrustes regression	Procrustes shape coordinates	CS	22.12 (F)	198	1	-	0.001	0.100
			Infection (general)	2.29 (F)		1	-	0.033	0.010
			Haplotype	1.10 (F)		1	-	0.324	0.005
			CS x Infection (general)	1.64 (F)		1	-	0.117	0.007
Model 3	Linear regression	CS	Intercept		182	1	2.86 ± 0.04	<0.001	0.162
			Haplotype: H2	2.20 (F)		1	−0.07 ± 0.04	0.140
			Infection: Plagiorchiida sp.	9.49 (F)		4	−0.07 ± 0.06	0.288
			Infection: Austrodiplostomum sp. 2				0.17 ± 0.08	0.046
			Infection: Petasiger sp. 5				0.42 ± 0.09	<0.001
			Infection: Echinoparyphium sp.				0.26 ± 0.09	0.003
Model 4	Procrustes regression	Procrustes shape coordinates	CS	19.44 (F)	182	1	-	0.001	0.094
			Infection (species-specific)	1.93 (F)		4	-	0.015	0.037
			Haplotype	0.92 (F)		1		0.435	0.004
			CS x Infection (species-specific)	2.02 (F)		4	-	0.006	0.039
Test 5	Pairwise comparison of Procrustes regression slopes (model 4)	Uninfected	Infection: Austrodiplostomum sp. 2	1.21 (Z)	182		-	0.112	-
		Uninfected	Infection: Petasiger sp. 5	2.23 (Z)			-	0.016	-
		Uninfected	Infection: Plagiorchiida sp.	1.24 (Z)			-	0.115	-
		Uninfected	Infection: Echinoparyphium sp.	−0.07 (Z)			-	0.542	-
		Infection: Austrodiplostomum sp. 2	Infection: Petasiger sp. 5	1.39 (Z)			-	0.080	-
		Infection: Austrodiplostomum sp. 2	Infection: Plagiorchiida sp.	1.89 (Z)			-	0.038	-
		Infection: Austrodiplostomum sp. 2	Infection: Echinoparyphium sp.	0.71 (Z)			-	0.250	-
		Infection: Petasiger sp. 5	Infection: Plagiorchiida sp.	2.08 (Z)			-	0.019	-
		Infection: Petasiger sp. 5	Infection: Echinoparyphium sp.	1.82 (Z)			-	0.041	-
		Infection: Plagiorchiida sp.	Infection: Echinoparyphium sp.	0.83 (Z)			-	0.194	-
Test 6	Wilcoxon test	CS	Petasiger sp. 5 ‘early’ vs. ‘late’	27 (W)	13		-	0.354	-
			Petasiger sp. 5 ‘early’ vs. uninfected	634 (W)	126		-	0.106	-
			Petasiger sp. 5 ‘late’ vs. uninfected	516 (W)	123		-	0.005	-
			Plagiorchiida sp. ‘early’ vs. ‘late’	56 (W)	25		-	0.267	-
			Plagiorchiida sp. ‘early’ vs. uninfected	643 (W)	129		-	0.963	-
			Plagiorchiida sp. ‘late’ vs. uninfected	582 (W)	132		-	0.072	-
Model 5	Procrustes regression	Procrustes shape coordinates	CS	4.45 (F)	13	1	-	0.005	0.272
			Petasiger sp. 5 ‘early’ vs. ‘late’	1.88 (F)		1	-	0.071	0.115
Model 6	Procrustes regression	Procrustes shape coordinates	CS	11.82 (F)	126	1	-	0.001	0.085
			Petasiger sp. 5 ‘early’ vs. uninfected	4.14 (F)		1	-	0.003	0.029
Model 7	Procrustes regression	Procrustes shape coordinates	CS	9.64 (F)	123	1	-	0.001	0.073
			Petasiger sp. 5 ‘late’ vs. uninfected	2.21 (F)		1	-	0.046	0.017
Model 8	Procrustes regression	Procrustes shape coordinates	CS	4.13 (F)	25	1	-	0.020	0.153
			Plagiorchiida sp. ‘early’ vs. ‘late’	0.81 (F)		1	-	0.497	0.029
Model 9	Procrustes regression	Procrustes shape coordinates	CS	11.96 (F)	129	1	-	0.001	0.086
			Plagiorchiida sp. ‘early’ vs. uninfected	0.58 (F)		1	-	0.751	0.004
Model 10	Procrustes regression	Procrustes shape coordinates	CS	7.33 (F)	132	1	-	0.001	3.856
			Plagiorchiida sp. ‘late’ vs. uninfected	1.35 (F)		1	-	0.234	0.730
Test 7	Wilcoxon test	CS	Co-infection (single vs. co-infection)	418 (W)	80	1	-	0.596	-
Model 11	Procrustes regression	Procrustes shape coordinates	CS	12.73 (F)	80	1	-	0.001	0.141
			Co-infection (single vs. co-infection)	0.61 (F)		1	-	0.77	0.007
			Haplotype	0.76 (F)		1	-	0.60	0.008
			CS x Co-infection (single vs. co-infection)	0.78 (F)		1	-	0.58	0.009

Impact of co-infection on shell phenotype

Due to the relatively low occurrence of co-infection and the variability of associations of trematode species in those co-infections, our current dataset does not present sufficient observations to allow us to resolve questions related to species-specific patterns of co-infections. Therefore, we decided to investigate whether co-infections in general impacted shell phenotype differently compared to cases of infection with a single trematode species. We compared shell size between both categories with a Wilcoxon test (Table 2: test 7). For shape, we used Procrustes regression with Procrustes shape coordinates as response variable and haplotype, infection status (single infection versus co-infection) and CS as response variables (Table 2: model 11).

Results

Characterization of the snail population

From a total of 419 Bulinus tropicus specimens we excluded 134 potential juveniles (Supplementary Fig. 1) and 28 severely damaged specimens. Of the remaining 257 specimens, 227 were successfully barcoded with our 492 bp fragment of COX1 (88.3%, GenBank accession numbers ON955535-ON955760). We identified 12 unique snail haplotypes differing from one another by one to four single nucleotide polymorphisms (SNPs; see Fig. 2 for haplotype network and Supplementary Table 3 for genetic distances). Two haplotypes (H1 and H2) that differ by one SNP accounted each for 104 specimens, representing 91.6% of the dataset; the others each occurred <5 × and were therefore excluded from the statistical analyses (19 specimens in total).

Fig. 2

Haplotype network of the 227 Bulinus tropicus specimens from Lake Kasenda that were barcoded for a fragment of COX1. Circle area is proportional to the number of specimens of each haplotype (smallest circle = 1 specimen, largest = 104). Single nucleotide polymorphisms between the haplotypes are represented by dashes on the connecting branches. For each haplotype the proportion of snails that are uninfected or infected (any trematode infection), as derived from RD-PCR results, is indicated in grey and white, respectively. Labels refer to the haplotypes used in Supplementary Table 3.

Characterization of the trematode species assemblage

RD-PCRs indicated that 99 of 257 snails were infected (38.5%). The number of quality-controlled reads obtained with HTAS for each of these 99 infected snails ranged from 0 to 53,294 (mean and standard deviation: 19,743 ± 7763). For 86 of these infected snails (86.9%), the HTAS workflow allowed to successfully characterize the infecting trematode(s). For the remaining 13 specimens, trematode infections could not be characterized, either due to a lack of reads (<1000) from Illumina sequencing (4 specimens) or because of failure to amplify or sequence the trematode markers. Those 13 snail specimens were therefore excluded from the statistical analyses (for 3 of these snails COX1 barcoding had failed additionally). Examination of control B. glabrata snails with known infection status revealed low levels of background cross-contamination, which was effectively removed during subsequent filtering (Supplementary Table 4). We deposited the quality-controlled sequences obtained in GenBank under the following accession numbers: ON939714-ON939728 for COIX1 I, ON939657-ON939697 for COX1 II, ON970175 – ON970210 for cytb, ON970211 – ON970251 for NAD1 and ON964584-ON964736 for ITS2.

Delimitation of species with ASAP resulted in correct classification for almost all COX1 I, COX1 II, cytb and NAD1 sequences (98–100%) from GenBank, but the performance for ITS2 was less good (92%, see Supplementary Table 2). Nevertheless, our results suggest reliable delimitation and 11 trematode taxa were molecularly delineated in our study (Table 1). Seven of these trematode taxa were not identified beyond their order or family as close matches were lacking from GenBank. The other four taxa were identified more accurately and linked to species found in other studies. Based on 99.8% sequence similarity in NAD1 we identified Petasiger sp. 5 as described by Laidemitt et al. (2019) from cercariae shedding from Bulinus ugandae, Bulinus globosus and Bulinus truncatus collected in Kenya. The COX1 sequences from Petasiger sp. 5 closely match (99.7% similarity) those of an unidentified echinostomate by Schols et al. (2020) in Bulinus globosus from Zimbabwe. Laidemitt et al. (2019) tentatively identified Petasiger sp. 5 morphologically as Petasiger variospinosus, a species infecting an amphibian as second intermediate and a bird as final host (King and Van As, 2000). Our second trematode species shares 98.3% NAD1 similarity with Echinoparyphium sp. described from cercaria shedding from a Kenyan B. tropicus by Laidemitt et al. (2019). Third, we report a taxon closely related to the diplostomid cercaria type 2 (99.5% identity in ITS2) found in Biomphalaria obstructa snails from catfish aquaculture in the U.S.A. (Rosser et al., 2016). Rosser et al. (2016) related this parasite to Austrodiplostomum sp. 2 from Locke et al. (2015) based on high similarity in COX1, which we also adopt here. Finally, we identified Euclinostomum heterostomum based on 99.1% similarity in ITS2 with metacercariae collected from singhi fishes in India (Athokpam et al., 2014) and 97.0% similarity in COX1 to metacercariae collected from cichlid fishes in Israel (Caffara et al., 2016). The life cycle of E. heterostomum involves a fish as second intermediate host and a piscivorous bird as final host (Dönges, 1974; Jhansilakshmibai and Madhavi, 1997). Hitherto, only Indoplanorbis exustus in India (Jhansilakshmibai and Madhavi, 1997) and Bulinus globosus in Nigeria (Dönges, 1974) were known as first intermediate host for E. heterostomum, but our results indicate that B. tropicus also hosts this trematode. Four trematode species (Plagiorchiida sp., Austrodiplostomum sp. 2, Petasiger sp. 5 and Echinoparyphium sp.) were each retrieved in >10 snail specimens, and jointly account for >80% of all detected infections. The remaining seven species were rarer (one to three occurrences), hampering statistical assessment of their impact on shell phenotype (Table 1). We detected 11 instances of interspecific co-infections: nine snails were infected by two trematode species, one by three species, and one by four species.

Table 1

Composition of the trematode community (11 taxa) in the population of Bulinus tropicus from Lake Kasenda. For each trematode species the overall prevalence of infection and the proportion of co-infections are indicated. Finally, we report which markers were successfully sequenced from our HTAS protocol for each trematode species (v = success; x = failure).

Trematode species	Prevalence (n = 257)	Proportion in co-infections	COX1 I	COX1 II	ITS2	NAD1	cytb
Plagiorchiida sp.	12.8%	24.2%	x	v	v	x	v
Austrodiplostomum sp. 2	8.2%	33.3%	x	x	v	v	x
Petasiger sp. 5	5.8%	13.3%	v	v	v	v	v
Echinoparyphium sp.	4.7%	0.0%	x	v	v	v	v
Euclinostomum heterostomum	1.2%	66.7%	v	x	v	x	v
Heterophyidae sp.	1.2%	66.7%	x	x	v	x	x
Allocreadiidae sp.	0.8%	100.0%	x	x	v	x	x
Strigeidae sp.	0.4%	0.0%	v	v	v	v	v
Crassiphialinae sp.	0.4%	0.0%	x	v	v	v	x
Allocreadiidae sp.	0.4%	100.0%	x	x	v	x	x
Digenea sp.	0.4%	100.0%	x	x	v	x	x

Analysis of the influence of trematode infections on shell phenotype

Our core dataset included 198 adult snails that were successfully barcoded, did not have rare COX1 haplotypes and, if infected, for which trematode characterization was successful. Snail haplotype did not influence general infection status, shell size or shape neither when considering general and species-specific trematode infections (Table 2: tests 1 and 2). Our Kruskal–Wallis test did not indicate significant differences in CS between the two snail haplotypes (Table 2: test 3; Supplementary Fig. 2). Finally, no differences in shell shape were found between haplotypes with permutational MANOVA (Table 2: test 4).

Our linear regression of infection status and snail haplotypes on CS showed a significant positive effect of trematode infection, but there was no significant effect of snail haplotype (Table 2: model 1, Fig. 3). The goodness of the fit of the model was poor, however, suggesting that the model explained only a low proportion of the variance of CS. Procrustes regression of infection status, haplotype and size on Procrustes shape coordinates indicated that shell size and infections both significantly affect shell shape, but snail haplotype and the interaction between infection status and CS had no significant effect (Table 2: model 2). Allometry explained a larger share of the variance in shell shape than did trematode infections (R2 of 0.100 and 0.010 respectively, Table 2: model 2), and the differences between uninfected and infected snails in shell shape are illustrated in the NMDS morphospace (Fig. 4). Compared to uninfected snails, infected snails have a more protuberant apex, and an inward-folded outer lip of the aperture (Fig. 5B).

Fig. 3

Variation in centroid size (CS) of Bulinus tropicus shells depending on infection status. Median CS values are illustrated as thick lines within the boxes, first and third quartiles as the upper and lower hinges. The whiskers extend from the hinges to the highest or lowest values (for upper and lower whisker, respectively) within 1.5 x the inter-quartile range of the corresponding hinge (Wickham 2016). Each dot represents an individual snail.

Fig. 4

Morphospace occupation plot reconstructed by non-metric multidimensional scaling (NMDS) on the Procrustes shape coordinates of 198 individuals of Bulinus tropicus from Lake Kasenda. Blue circles indicate uninfected specimens, orange symbols infected specimens, with ∇ = single infection by Echinoparyphium sp.; Δ = single infection by Austrodiplostomum sp. 2; □ = single infection by Plagiorchiida sp. I; ◊ = single infection by Petasiger sp. 5; ○ = any other infection. Filled symbols indicate the mean shape for each infection group, and the filled orange circle represents the mean of all trematode-infected snails. Vectors indicate the shape changes by any of these infection types compared to uninfected specimens. These differences are further illustrated in Fig. 5. (For interpretation of the references to colour in this figure legend, the reader is referred to the Web version of this article.)

Fig. 5

Difference in mean shell shape between an uninfected Bulinus tropicus (A), an infected B. tropicus (all types of infections collectively; B), and a B. tropicus infected by Petasiger sp. 5 (C). Each of the displayed specimens is the specimen nearest to the group mean in morphospace (Fig. 4). The differences in Procrustes shape coordinates among these specimens are visualized with vectors (vector length magnification = ×2.5).

To investigate whether the influence of infections on shell phenotype varied among trematode species, we excluded the 16 snails that are infected by rare trematodes or that have co-infections (see Material and Methods). Thus, the subsequent analyses are based on 182 snails, each either uninfected or infected by one of our four common trematode species: Echinoparyphium sp., Austrodiplostomum sp. 2, Petasiger sp. 5 and Plagiorchiida sp. Linear regression of trematode species-specific infections and haplotype on CS indicated that snail specimens infected by Austrodiplostomum sp. 2, Echinoparyphium sp. and Petasiger sp. 5 were significantly larger than uninfected snails (on average 6, 9 and 15% larger, respectively), whereas Plagiorchiida sp. did not differ significantly in CS from uninfected specimens (Table 2: model 3, Fig. 3). Snail haplotype had no significant effect. Assessing the effect of trematode infections on CS by considering individual trematode species improved the goodness-of-fit of the model compared to when all infections are treated collectively (R2 = 0.162 and 0.049 respectively, Table 2: models 2 and 4). Procrustes regression of trematode-specific infections, CS, and haplotype on Procrustes shape coordinates indicated that CS, trematode-specific infections, and their interaction all had a significant impact on shell shape, whereas snail haplotype did not (Table 2: model 4). CS explained a larger share of the variance of Procrustes coordinates (R2 = 0.094) compared to species-specific trematode infections (R2 = 0.037) and the interaction between CS and infections (R2 = 0.039). Pairwise comparisons of the Procrustes regression of shell shape revealed significant differences among the allometric slopes of uninfected snails and snails infected by Petasiger sp. 5; among snails infected by Petasiger sp. 5 versus those by Plagiorchiida sp.; among those by Petasiger sp. 5 versus those by Echinoparyphium sp.; and among those by Austrodiplostomum sp. 2 versus those by Plagiorchiida sp. (Table 2: test 5). Differences in allometric trajectories for the first and second NMDS shape axes for the various infecting species are illustrated in Supplementary Fig. 3. Finally, significant differences in shell shape are illustrated in our morphospace occupation plot by proportionally larger distances in the morphospace among the mean shape of snails infected with Petasiger sp. 5 versus that of uninfected snails, or among snails infected by Plagiorchiida sp. or Echinoparyphium sp. (Fig. 4). Compared to uninfected snails, those infected by Petasiger sp. 5 had a more protuberant apex, the inward-folded outer lip of the aperture, and an umbilicus displaced further away from the basal margin of the shell (Fig. 5C). Based on trematode read proportions, we classified 11 Plagiorchiida sp. infections as ‘late’ and 14 as ‘early’, and 5 Petasiger sp. 5 infections as ‘late’ and 8 as ‘early’. For Echinoparyphium sp. and Austrodiplostomum sp. 2, we observed dramatically lower read proportions (means and standard deviations of 3.1 ± 8.3% and 0.3 ± 0.5% respectively) and could not determine the infection stage. For Plagiorchiida sp. and Petasiger sp. 5, we detected no significant differences in shell size and shape between ‘early’ and ‘late’ infections by either species, but snails with early infections have intermediate sizes and shapes compared to uninfected snails and snails with ‘late’ infections (Supplementary Figs. 4 and 5). Only snails with ‘late’ infections by Petasiger sp. 5 were significantly larger than uninfected snails (Table 2: test 6)). Snails with ‘late’ infections by Plagiorchiida sp. had shells that trend towards significantly smaller size compared to uninfected snails (Table 2: test 6, Supplementary Fig. 4). Finally, shell shape did not significantly differ between ‘early’ and ‘late’ infections of either species (Table 2: models 5 and 8, Supplementary Fig. 5). Snails with ‘early’ or ‘late’ Pagiorchiida sp. infections did not differ in shell shape from uninfected snails (models 9 and 10) whereas the shell of snails with both ‘early’ and ‘late’ Petasiger sp. 5 infections differed significantly from those with no infection (models 6 and 7). Differences in the displacement along the second axis of the NMDS morphospace (which shows independent allometric change in shape for Petasiger sp. 5, Supplementary Fig. 3) were more important for ‘late’ versus ‘early’ infections (Supplementary Fig. 5).

Impact of co-infection on shell phenotype

Our dataset to evaluate the impact of co-infections on shell phenotype included 69 specimens that were infected by a single trematode species and 11 with co-infections. The Wilcoxon test as well as linear regression indicated that CS did not differ significantly between snails infected by one trematode species and those co-infected by multiple species (Table 2: test 7, Fig. 3). Procrustes regression showed that shell shape was significantly influenced by CS, but that there was no significant difference in shell shape between snails infected with a single versus multiple trematode species (Table 2: model 11). Snail haplotypes again had no significant effect.

Discussion

Trematode infections in Bulinus tropicus from Lake Kasenda

The overall trematode infection prevalence (38.5%) and trematode species richness (11 taxa) observed in B. tropicus from Lake Kasenda were both high. However, care should be taken when comparing our results with values from other studies because various methods were used to diagnose and identify infections, which have different levels of sensitivity. More specifically, cercarial shedding may underestimate trematode prevalence because it does not detect immature infections and it is less sensitive than RD-PCR (Born-Torrijos et al., 2014; Carolus et al., 2019) and high-throughput amplicon sequencing (HTAS, Hammoud et al., 2022). Furthermore, studies that rely on morphological identification of cercariae may underestimate species richness since diagnostic morphological characters are scarce in larval trematodes (Nolan and Cribb, 2005). Using cercarial shedding and morphological identification, Loker et al. (1981) reported infections by 2–7 parasite morphotypes with a prevalence ranging from 2.5 to 9.7% for several species of Bulinus from Tanzanian streams (sample size ranging from 186 to 1503 snails). Chingwena et al. (2002) analysed 4080 B. tropicus collected from streams and dams in Zimbabwe, 13.1% of which were infected by combinations of 5 morphotypes of cercariae, and Mohammed et al. (2016) reported that 46.2% of 1403 B. truncatus collected form irrigation canals in Sudan were infected by an assemblage of 8 morphotypes of cercariae. Beyond sampling size, various factors likely drive variation in observed trematode species richness among snail populations, including differences in susceptibility to infection among snail species and final host species richness (Hechinger and Lafferty, 2005). Discussing these extensively is beyond our scope, but the infection dynamics reported here combined with the properties of the study system mentioned above render Bulinus tropicus from Lake Kasenda promising to study trematode-intermediate host dynamics in detail. Although our molecular approach towards linking parasites and their intermediate hosts has high sensitivity and allows accurate characterization of trematode taxa, we were still limited in producing specific trematode identifications, because trematodes are underrepresented in public DNA databases (see also Schols et al., 2020). We are, however, optimistic that integrative molecular and morphological characterizations of trematodes (e.g., Laidemitt et al., 2019) will improve future species-level identification, and trematode representation in existing databases. Despite this limitation, our study indicates that relying on multiple molecular markers increased our chances to find close matches in DNA databases and it allowed us to link sequences of different markers belonging to the same trematode species, which can subsequently help attributing GenBank data of unverified/uncertain origin to trematode genera or even species a posteriori.

Factors influencing Bulinus tropicus shell phenotype

We assessed the effect of trematode infections on B. tropicus shell phenotype while controlling for the influence of snail genotype and environmental factors by sampling a single homogeneous lake habitat. Barcoding of the B. tropicus population from Lake Kasenda revealed overall low genotypic diversity (p-distances ranging from 0.002 to 0.008, Supplementary Table 3), with most of the 227 analysed snails belonging to two genetically similar haplotypes. One of these haplotypes (H2) had been detected in Lake Kasenda by Tumwebaze et al. (2019). However, Nalugwa et al. (2010) earlier reported three B. tropicus haplotypes from Lake Kasenda which all differ substantially from those we found (p-distance ranging from 0.012 to 0.018, Supplementary Table 3). Those differences suggest that Nalugwa and colleagues have sampled a different neighbouring system, as confusion exists in lake nomenclature in the Ndali-Kasenda crater lake cluster (obs. pers.), and the provided geographic coordinates map to the shore of the Kazinga Channel. Alternatively, and much less likely, the B. tropicus population in Lake Kasenda may have been completely replaced during the ∼10 years separating both studies. In any case, our genetic data indicate a homogenous population of B. tropicus in Lake Kasenda with limited genotypic diversity. Furthermore, snail phenotype and the probability of trematode infection do not differ significantly among snail haplotypes.

Shell size and shape in B. tropicus differed significantly depending on trematode infections. More specifically, compared to uninfected specimens, infected snails (by all trematode taxa collectively) had larger shells, slightly more protuberant apices, and an inward-folded outer apertural lip. In pulmonate snails, the association between increased body size and trematode infection had previously mostly been documented in laboratory experiments (e.g., B. truncatus and B. senegalensis artificially infected with S. haematobium; Fryer et al., 1990), and it has only been reported in wild populations of a few species such as Lymnaea stagnalis (Żbikowska and Żbikowski, 2005) and Galba truncatula (Chapuis, 2009). Multiple mechanisms have been proposed to explain the higher incidence of trematode infections in larger snails. The most frequent explanations (see Sorensen and Minchella, 2001 for a comprehensive review) are 1) castration of the snail due to destruction of gonadal tissue by asexual stages of trematodes may reallocate metabolic resources from reproduction to growth (gigantism, Sousa, 1983) or 2) larger snails may be older and, therewith, would have been exposed longer to parasites, increasing the likelihood of infection (Sousa, 1983; Sorensen and Minchella, 2001). Examining the effect of specific trematode infections on body size in B. tropicus we found that snails infected by Petasiger sp. 5, Echinoparyphium sp. or Austrodiplostomum sp. 2 had larger shells than uninfected snails, whereas those infected by Plagiorchiida sp. did not. Secondly, comparing ‘early’ (<25 days, <45% of reads attributed to the infection) and ‘late’ infections (>25 days, >45% of reads) with Petasiger sp. 5 and Plagiorchiida sp., we found that snails with early infections have sizes intermediate to uninfected snails and snails with ‘late’ infections. As our results indicate specific relationships between shell size and species-specific trematode infections, as well as to the developmental stage of infection, the hypothesis that size changes are caused by trematode infections is more likely in our study system. Sorensen and Minchella (2001) argued that infections by trematode species that produce rediae (an asexual stage with a mouth, feeding directly on host tissue) cause gigantism more often than species producing sporocysts (an asexual stage without mouth that feeds by absorbing nutrients through the body wall), whereas Probst and Kube (1999) proposed the opposite. Petasiger sp. 5 and Echinoparyphium sp. likely produce rediae (King and Van As, 2000; Huffman and Fried, 2012) whereas Austrodiplostomum sp. 2 probably produces sporocysts (Cribb et al., 2003). The identification of Plagiorchiida sp. is taxonomically too coarse to deduce the type of asexual stage. The seven other encountered trematodes were too rare to assess their influence on shell size. In any case, the type of asexual stage does not appear to solely explain why infection by the three above-mentioned trematode species may have led to larger shell size whereas infection by Plagiorchiida sp. did not.

Our significant correlation between shell shape and centroid size indicates that shape changes with growth, i.e., we observe substantial allometric changes, as is common in many organisms (Zelditch et al., 2004). Therefore, we included allometry into our Procrustes models to assess the relative importance of shape changes explained by other variables than shell size (Outomuro and Johansson, 2017). As verification, we have recuperated the residuals from a regression of shell shape to size to retest all Procrustes models constructed with the residuals (instead of raw shape variables) as response variable, which produced highly similar results (not shown). Similar to the abovementioned differences in shell size we found significant differences in shell shape between uninfected and infected snails and the effect of infection on shell shape was found to be parasite-specific. We also observed that in snails infected by Petasiger sp. 5 and Plagiorchiida sp., snails with ‘early’ infections are on average of intermediate shape to those of uninfected specimens and those with ‘late’ infections. Upon examining the effects of ‘early’ and ‘late’ infections on shell size and shape, a first important constraint is that this distinction decreased sampling size and therewith the power of statistical test. Nevertheless, ‘early’ infections were always found to be intermediate to uninfected snails and snails with ‘late’ infections, which is consistent with the hypothesis that prolonged exposure to parasite infections results in accumulated changes in shell size and shape. Whereas our study is mainly correlative in nature, the differences between various infection stages are consistent with expectations under a causal link. Our results on the species-specific effects of infections on shell morphology are congruent with the findings of Żbikowska and Żbikowski (2005) who showed that infections by different trematode species correlated with differences in shell height/width ratio in Lymnaea stagnalis populations from Poland. The shape changes observed under Petasiger sp. 5 infections were substantial and of similar magnitude to those documented for Littorina saxatilis infected by Microphallus piriformes (McCarthy et al., 2004), Lymnaea stagnalis infected by various trematode species (Żbikowska and Żbikowski, 2005) or Cominella glandiformis infected by Curtuteria australis (Thieltges et al., 2009). We report (to our knowledge) the first changes in the shape of Bulinus shells related to trematode infection. The effect of infection on shell phenotype differed among trematode species, which strongly suggests that these trematodes exploit their host differently. These differences may either be the by-products of the pathology caused by infection or manipulative changes allowing the parasite to increase its reproductive potential (Hay et al., 2005). For example, McCarthy et al. (2004) showed that morphological changes in the shell of the marine snail Littorina saxatilis caused by infection with Microphallus piriformes increased the space within the shell, without increasing the volume of snail tissues and organs, i.e., the extra space is used for asexual reproduction of the parasite. Thus, these authors suggested that this phenotypic alteration may be triggered by the parasite to improve its reproductive potential without affecting the viability of its host. We suspect that the elongation of the apex of B. tropicus snails infected by Petasiger sp. 5 may play a similar role. Our results corroborate the hypothesis that trematode infections alter shell size and shape in their intermediate host and that the patterns are highly specific to trematode species. Further work, including experimental infections, is required to assess the physiological link between trematode infection and intermediate host morphological alterations in more detail.

In conclusion, we used state-of-the-art methods to reliably characterize trematode (co-)infections and to document shell phenotype, which allowed us to accurately document how trematode infections affect B. tropicus shell size and shape under natural conditions. Interestingly, effects are highly trematode-specific and accumulate with the time since infection. This level of complexity may explain the contrasting results in the variety of studies on the impact of trematode-snail interactions on snail phenotype (see above). As new molecular methods enable detailed documentation of the link between parasites and intermediate hosts, it becomes possible to decipher the variety of interactions and mechanisms at play. Given the high complexity of multiple interactions, our results incline us to anticipate large variability in shell phenotypic changes linked to co-infections, which can only be disentangled with larger sampling sizes.

Declaration of competing interest

The authors of this manuscript declare that they have no competing interests.