Limited impact of schistosome infection on Biomphalaria glabrata snail microbiomes

TL;DR

This study finds that infection with the blood fluke Schistosoma mansoni has limited effects on the microbiome of Biomphalaria glabrata snails, which are important in the spread of schistosomiasis.

Contribution

The novel finding is that schistosome infection does not significantly alter snail microbiome composition or bacterial density, contrary to expectations.

Findings

Schistosome infection had no effect on bacterial density in snail hemolymph and hepatopancreas.

Microbiome composition was only slightly affected during the pre-patent period of infection.

Sample type (hemolymph vs. hepatopancreas) was the main factor influencing microbiome differences.

Abstract

The microbiome of disease vectors can be a key determinant of their ability to transmit parasites. Conversely, parasite infection may modify vector microbiomes. We explore the interactions between the Biomphalaria glabrata snail microbiome and the blood fluke Schistosoma mansoni, responsible for an estimated 200,000 human deaths each year. We have previously shown that the snail hemolymph (i.e. blood) and organs harbor a diverse microbiome. Here, we investigate the impact of schistosome infection on snail microbiomes, hypothesizing that invading schistosomes can alter the snail microbiomes in both composition and abundance over the course of infection, as developing schistosome parasites are in close contact with the host tissues.

We generated cohorts of uninfected and S. mansoni-infected snails. We collected snail hemolymph and hepatopancreas (i.e. liver) at eight timepoints during the pre-patent and patent periods of schistosome infection. We quantified bacterial density using qPCR and profiled the microbiome composition of all samples by sequencing the V4 region of the 16S rRNA.

Surprisingly, schistosome infection had no effect on bacterial density and limited effect on the microbiome composition, affecting mainly the hemolymph during the pre-patent period (at days 7 and 21). Organ and hemolymph microbiomes were relatively stable over time for both infected and uninfected snail cohorts. The sample type (hemolymph, hepatopancreas) was the major driver of the differences observed in microbiome composition.

The limited impact of schistosome infection on the host snail microbiomes might be explained by the long-term interaction of the two partners. Further investigations into the interactions among snails, their microbiomes and schistosome parasites may suggest strategies to disrupt the parasite lifecycle and, consequently, schistosomiasis transmission.

The online version contains supplementary material available at 10.1186/s13071-026-07299-z.

Figures

Experimental design. *Biomphalaria glabrata* Bg26 snails were randomly assigned to the exposed (*n* = 55) or control (*n* = 55) cohort. Six snails were sampled on D0 prior to miracidial or mock exposure in both the exposed and the control cohorts. Snails were then sampled weekly, through the pre-patent (*n* = 8 snails per cohort—D7 to D21) and patent (*n* = 6 snails per cohort—D28 to D39) periods. At these timepoints, the hemolymph and hepatopancreas of each snail were collected under sterile conditions, along with an environmental water sample from each tray

Hemolymph bacterial densities between control and schistosome-exposed snails. The longitudinal evolution of the number of 16S rRNA gene copies per microliter of snail hemolymph (D0–D39, measured by qPCR) showed no difference between cohorts at any of the timepoints. The dashed line marks the period before snails were exposed to miracidia (D0) and after snails had been exposed (D7 onwards). Statistical differences within cohorts at each of the timepoints were assessed using Kruskal-Wallis and Dunn’s post hoc tests for non-parametric data. Letters generated from Dunn’s post hoc tests are shown a

Longitudinal comparison of α-diversity metrics between control and schistosome-exposed snails. Boxplots showing values for A observed ASVs, B Simpson evenness and C Faith’s phylogenetic diversity between control (green) and exposed (red) snail cohorts at each of the eight timepoints, and for each sample type (Hm: hemolymph, Hp: hepatopancreas, W: water). The longitudinal analysis of the observed ASVs and the Faith’s PD metrics reveals specific timepoints where schistosome infection significantly impacts the microbial composition of hemolymph (D21 and D39) or hepatopancreas (D28 and D39). The d

Longitudinal taxonomic diversity between control and schistosome-exposed snails. The relative abundance of the top 20 most dominant phyla is shown for both control (left) and exposed (right) cohorts at each of the eight timepoints, including snail tissues (hemolymph and hepatopancreas) and environmental water samples. Samples were merged by sample type and day: each bar represents two water samples and 6–8 snail samples for hepatopancreas or hemolymph. The taxonomic compositions of both control and exposed snails are relatively stable over time, dominated by phyla Proteobacteria, Bacteroidetes

Full text

Backgrounds

The microbiome plays numerous roles in host fitness and functions, including in metabolism [1–3], immune system development and homeostasis [4, 5] and protection against infections or diseases [5, 6]. While host microbiomes have been shown to impact invading pathogens [7, 8], pathogen infection can also influence the composition of the host microbiome [9, 10]. However, the impact of pathogens on the host microbiome may vary significantly depending on whether the pathogens are transient or engage in long-term interactions with the host (i.e. parasitic relationships), especially when pathogen fitness is closely tied to host survival [9, 11–15]. The effect of parasites on the host microbiomes has mainly been investigated in vertebrates, with a strong emphasis on gut microbiomes [16]. While parasites often use invertebrate vectors for transmission, their impact on invertebrate vector microbiomes has received comparatively less attention. Investigating the impact of parasites on invertebrate microbiomes across a wide variety of host-parasite systems, which include a diverse range of vector species (insects, crustaceans, mollusks, etc.) and varying interaction durations, could provide valuable insights into the dynamics of host microbiomes during long-term infections. A better understanding of these interactions could open new avenues for controlling parasite transmission, as demonstrated by the use of Wolbachia bacteria in mosquitoes to manage dengue virus [17, 18].

Biomphalaria snails are an excellent model to understand the impact of parasites on host microbiomes. These snails have diverse microbiomes [19, 20] and are vectors of a wide range of parasites [21, 22], which colonize their host for weeks [23]. Among these parasites, the schistosome blood flukes are of significant biomedical importance, infecting 230 million people in tropical and subtropical areas and responsible for approximately 200,000 deaths per year [24]. Infected people release parasite eggs from which miracidia larvae hatch [25]. These miracidia then actively search for a snail host and, once it has been found, they penetrate its head foot. The parasite then develops for several weeks within its host. This development is divided into a pre-patent and patent periods [26]. The pre-patent period corresponds to the first 4 weeks of infection, where the miracidium transforms into primary sporocyst at the site of penetration, and starts producing second-generation sporocysts, which will be released after 10–11 days [27]. The secondary sporocysts then migrate through snail organs to the hepatopancreas (liver-like organ) and ovotestis [28]. These sporocysts reside within these organs and, when mature, produce either additional sporocysts or cercariae, the free-swimming larvae infecting humans. The cercariae are shed from the snail host by breaking through the snail tissues [28]. The start of cercarial shedding marks the beginning of the patent period, which occurs around the 4th–5th week of infection.

During their development within its snail host, parasites are exposed to the snail microbiomes, which are not just diverse but also vary across organs and tissues [19, 20]. For instance, the hepatopancreas, where parasites grow, and the hemolymph (i.e. blood), which bathes organs and parasites, have dramatically different microbiomes [19]. Most microbiome studies of Biomphalaria and other parasite-transmitting snails have focused on whole snails [29–34]. However, our recent characterization of hemolymph, organ and whole snail microbiomes revealed that whole snail microbiomes are composite and do not accurately represent the microbiome compositions of the different parts of the snail [19, 20]. Whole snail microbiomes could, therefore, mask subtle changes or generate spurious shifts in microbiome diversity and composition.

Only a few studies have investigated the impact of schistosome and other trematode infections on the snail microbiome. Portet et al*.* characterized the microbiomes of whole snails exposed to schistosomes at three timepoints during the pre-patent period (day 1, 4 and 25 after exposure) and did not reveal a significant disturbance in α-diversity during this period [32]. McCann et al. examined the whole snail microbiomes of field-collected Galba truncatula and found differences between uninfected snails and those infected with liver fluke larvae [31]. However, the impact of parasites on specific tissues and organs during the pre-patent and patent periods is unknown.

To better understand the dynamics within relevant snail organ compartments, we characterized the impacts of Schistosoma mansoni infection on the microbiome of Biomphalaria glabrata hemolymph (to which parasites are exposed) and hepatopancreas (where parasites reside) over the course of infection. We hypothesize that parasite infection progressively and significantly affects the microbiomes of its snail host over the course of the infection. To test this, we performed weekly collections of infected and non-infected snails over a 6-week period, quantifying the bacterial burden and characterizing the microbiome composition of sampled tissues. Our results showed a limited and specific impact of schistosome infection on the composition of the hemolymph microbiome.

Methods

Snail infection

We used Biomphalaria glabrata snails from the inbred line 26 (Bg26, derived from 13 to 16-R1 line [35]) to limit the impact of the host genetic variation on our study. Snails were reared in 10-gallon aquariums with well water at 26–28 °C on a 12-h light cycle and fed ad libitum on green lettuce [26]. All snails used in the study had a shell diameter between 5 and 7 mm and were randomly sampled from two aquariums and acclimated in trays (55 snails per tray) 7 days before exposure to parasites. We prepared miracidia by emulsifying livers of two Syrian golden hamsters as previously described [36]. We exposed Bg26 snails (n = 110) to 5 SmLE miracidia in 1 ml freshwater in a 24-well plate under artificial light on D0 overnight [36]. We prepared a separate cohort of 110 control snails into plates to undergo the same procedure of exposure without miracidia (Fig. 1). The next morning, we moved snails to trays with well water and green lettuce. Each cohort had 2 trays, and each tray had 55 snails (control trays: CA, CB; exposed trays: EA, EB). We kept trays in a temperature-controlled room at 26–28 °C on a 12-h light cycle. During normal maintenance of laboratory-infected snails, we covered the trays with a black plexiglass lid starting on week 3 post-exposure to parasite to prevent cercarial shedding from infected snails. To better mimic natural field conditions, the experimental snails were covered only with a clear plexiglass lid for the duration of the experiment. Well water and lettuce were replaced daily. Lettuce was extensively washed with tap water and spun dry. To minimize the impact of water or lettuce-associated microbes, the same batches of water and lettuce were used for all control and exposed trays on a given day. Snails were visually assessed for infection during dissection, starting at day 21, by looking at the presence of schistosome sporocysts or cercariae in ovotestis and hepatopancreas. The infection rates for D7 and D14 were estimated based on the overall infection rate observed from D21 to D39.Fig. 1. Experimental design. Biomphalaria glabrata Bg26 snails were randomly assigned to the exposed (n = 55) or control (n = 55) cohort. Six snails were sampled on D0 prior to miracidial or mock exposure in both the exposed and the control cohorts. Snails were then sampled weekly, through the pre-patent (n = 8 snails per cohort—D7 to D21) and patent (n = 6 snails per cohort—D28 to D39) periods. At these timepoints, the hemolymph and hepatopancreas of each snail were collected under sterile conditions, along with an environmental water sample from each tray

Hemolymph and hepatopancreas collection

Snails were collected at eight timepoints during the course of infection, on D0 (prior to exposure), D7, D14, D21, D28, D35, D37 and D39. For each collection, 3–4 snails were sampled from each tray (6–8 snails per cohort; 12–16 snails per timepoint). We sampled 500 µl of water from each tray at these timepoints (environmental samples). Each snail was then processed individually. We wiped the shell three times with 70% ethanol for disinfection. The snail was killed by complete withdrawal of the hemolymph via heart puncture using a sterile 1 ml syringe and 23.5-gauge needle [20]. We then gently crushed the snail shell between two sterile microscope slides and removed shell pieces with sterilized tweezers. We collected the hepatopancreas (liver) and rinsed it three times with 1 ml of sterile water. Each collected sample was immediately placed in a sterile 1.5 ml pre-cooled microtube stored on dry ice.

16S rRNA gene library preparation and sequencing

We randomized all samples prior to both DNA extraction and triplicate PCR assays for library preparation to minimize any potential batch effects. Samples were randomized based on factors including tray, snail number, sample type (water, hemolymph, hepatopancreas) and day of collection. We extracted DNA from the snail and tray water samples using the DNeasy Blood and Tissue Kit (Qiagen) with unmodified buffers from the kit and following the manufacturer’s protocol with minor adjustments [20]. “Kitome” controls were included to check for any potential contamination in the DNA extraction kit [37]. We used 40 µl of hemolymph and environmental water sample for each DNA extraction. We manually homogenized the hepatopancreas organs using sterile micropestles prior to extraction. Samples were incubated for 1 h at 56 °C in a water bath. We recovered the final gDNA in 50 µl of elution buffer.

The 16S rRNA gene libraries were prepared in a biosafety cabinet using sterile equipment and materials to avoid contamination. We performed triplicate PCR for each sample to limit potential amplification biases introduced during the PCR reaction. The 16S rDNA V4 region was amplified using primers (515f and 806rB) from the Earth Microbiome Project [38]. Each 515f primer was barcoded for sample identification. We followed the methods described in Chevalier et al*.* [20] for generating the 16S V4 libraries, with some modifications: For each PCR reaction mix, we used 1 µl sterile water, 5 µl AccuStart II PCR Supermix, 1 µl of each primer at 2 µM and 2 µl of DNA template. Libraries were purified using KAPA Pure Beads and quantified using Picogreen assay (Invitrogen) following the manufacturer’s protocol. Two equimass library pools (124 samples each) were made using 30 ng of each library. These pooled libraries were quantified by qPCR with the KAPA library quantification kit, following the manufacturer’s protocol, and sequenced by Admera Health on an Illumina MiSeq platform (250-bp paired-end sequencing). Raw sequencing data are accessible from the NCBI Sequence Read Archive under BioProject PRJNA1171869.

Measure of bacterial density

We measured bacterial density within the hemolymph by amplifying the V4 region of the 16S rRNA gene using primers 515f and 806rB (without adaptors and barcodes) by qPCR as previously described [20]. We determined absolute 16S rRNA gene quantities by using standard curves of six tenfold dilutions (6 × 10^1^ to 6 × 10^7^) of purified 16S rRNA gene PCR product from E. coli. Results were normalized by the amount of hemolymph used in the DNA extraction (40 µl) and the qPCR dilution factor (1/10th dilution). Samples with a normalized quantity > 500,000 copies (outliers in the data distribution) were excluded from the analysis.

Bioinformatic and statistical analysis

We performed sequence processing and analysis using QIIME2 (v2021.4) [39] and R (v4.3.1) [40]. Scripts used for sequence processing and analysis are available in a Jupyter notebook (10.5281/zenodo.18262726). We processed the sequencing data as previously described [20]. Briefly, we denoised and clustered the sequences into amplicon sequence variants (ASVs) using the dada2 module with a maximum expected error of 5. We denoised only the forward reads because of the low quality of the reverse reads. We determined the taxonomy of the ASVs using the SILVA database (release 132). We blasted ASVs with unassigned taxonomy against the NCBI nt database using megablast from BLAST+ to identify eukaryotic contaminants for removal in downstream analysis. We aligned the ASVs and masked the highly variable positions using the mafft and mask commands from the alignment module to build a phylogenetic tree. Finally, we built and rooted a tree from the alignment using the fasttree and midpoint-root commands.

We imported the QIIME2 files in R using the qiime2R package (v0.99.6) and converted them into phyloseq objects using the phyloseq package (v1.46.0). We assessed alpha diversity using the number of observed ASVs and Simpson evenness using the package microbiome (v1.24.0). None of the α-diversity data followed a normal distribution (Shapiro-Wilk test, P < 0.05). We therefore used a non-parametric (Wilcoxon or Kruskal-Wallis and Dunn’s post hoc) test for the statistical comparison of the data. We assessed β-diversity on rarefied data, as best practices suggest [41], using R packages vegan (v2.6-4) and pairwiseAdonis (v0.4.1) to perform β-dispersion and PERMANOVA tests including factors cohort and day with 1000 permutations. We performed rarefaction and statistical tests for each sample type (hemolymph, hepatopancreas, water). We assessed the data for homogeneity of variances by computing the centroid distances and assessing differences using the Tukey HSD post hoc test. Additive and interaction linear mixed models including fixed factors (cohort, sample type, day) and random factors (snail_ID, tray) were assessed for each α- and β-diversity index using the packages stats (v4.3.1), lme4 (v1.1-34) and lsmeans (v2.30-0). We used the Bayesian information criterion (BIC) scores for each statistical model tested to determine the best fitting models to report.

We also assessed bacterial density data using the non-parametric Wilcoxon or Kruskal-Wallis and Dunn’s post hoc tests for statistical comparisons.

Results

Experimental design and library statistics

We sampled the hemolymph and hepatopancreas of B. glabrata snails uninfected and infected with S. mansoni parasites over the course of 39 days (Fig. 1). For each of the eight timepoints, we collected 6–8 uninfected (control, unexposed cohort) and infected (exposed cohort) snails. The snail infection rate was 100% after visual assessment of the hepatopancreas during dissection starting D21 (presence of daughter sporocysts) and onwards (presence of cercariae). We assumed snails from prior D21 were also 100% infected in the S. mansoni exposed cohort as it is unlikely that only uninfected snails were sampled during the first 2 weeks. Such a high infection rate was expected given the high susceptibility of Bg26 snails to our SmLE population [36]. The survival rates of snails not sampled for microbiome analysis were similar between replicates within each treatment group (control or S. mansoni-exposed) [(CA: 11/22, 50% survival; CB: 15/22, 68.18% survival; log rank test = 2.2, P = 0.1) (EA: 3/24, 12.5% survival: EB: 2/24, 8.33% survival; log rank test = 0, P = 1)] (SuppFigure1). However, S. mansoni infection had a significant impact on survival, with control snails having a higher survival rate than the exposed snails, as expected (log rank test = 13.8, P = 2e-04) [42]. We observed a marked increase in mortality in the exposed snails around D30.

A total of 216 snail samples (108 hemolymph and 108 hepatopancreas) and 32 environmental samples (tray water) were collected and processed for microbiome library preparation. We used high-throughput sequencing of the V4 regions of the 16S rRNA gene with the Illumina MiSeq platform and obtained a total of 14,428,659 reads from 248 samples with an average of 58,180 ± 1406 (mean ± SE) input reads per sample. After filtering, we had 9,836,712 reads with an average of 39,664 ± 918 filtered reads per sample (SuppFigure2, SuppTable1). We identified 4923 ASVs and removed ASVs matching mitochondria and chloroplasts (82 ASVs, 1.7% of ASVs) or eukaryotes (149 ASVs, 3.1% of ASVs). Among the 4692 ASVs retained, 43.27% (2030) had unassigned taxonomy. Rarefaction curves showed that all samples reached a plateau, confirming that our sequencing effort was sufficient to capture all the bacterial diversity present in our samples (SuppFigure3).

Bacterial density is stable in uninfected and infected snail hemolymph over time

We measured the impact of S. mansoni infection on the bacterial density of snail hemolymph by comparing bacterial loads between infected and uninfected snails at each timepoint using qPCR. Infected snails showed stable bacterial density in the hemolymph over time (Kruskal-Wallis test: X^2^ = 13.19, df = 7, P = 0.0676; Fig. 2). Bacterial density measured in the hemolymph of the uninfected snail cohort also remained relatively stable, with the only statistical difference observed between days 28 and 37 (Kruskal-Wallis test: X^2^ = 17.171, df = 7, P = 0.0163, comparison between D28 vs. D37: P = 0.0482; Fig. 2), likely because of lower variation in bacterial density at day 28. In addition, there were no significant differences at any given timepoint between the infected and uninfected cohorts (Table 1; Wilcoxon test: W = 10–27, P = 0.27–1). Overall, schistosome infection did not impact the bacterial density in the hemolymph. In addition, the hemolymph microbiome showed a strong temporal stability, suggesting tight control of the microbiome by the snail.Fig. 2. Hemolymph bacterial densities between control and schistosome-exposed snails. The longitudinal evolution of the number of 16S rRNA gene copies per microliter of snail hemolymph (D0–D39, measured by qPCR) showed no difference between cohorts at any of the timepoints. The dashed line marks the period before snails were exposed to miracidia (D0) and after snails had been exposed (D7 onwards). Statistical differences within cohorts at each of the timepoints were assessed using Kruskal-Wallis and Dunn’s post hoc tests for non-parametric data. Letters generated from Dunn’s post hoc tests are shown above the boxplots. Groups under the same letter are not statistically different: Groups marked “ab” are not statistically different from “a” or “b”, while groups marked “a” are statistically different from groups marked “b”. Control: green and simple letters. Exposed: red and letters with prime symbolTable 1Bacterial density statisticsDayWP0121.0007210.27014270.79821211.00028160.92735170.93637140.57539100.749

Similarly, we attempted to assess the impact of schistosome infection on the bacterial density of the hepatopancreas organ, where schistosome parasites reside. However, we were unable to amplify the snail reference gene (piwi4) in most of the samples for undetermined reasons, preventing normalization of the data.

Schistosome parasite infection has minimal impact on the microbiome diversity of its snail host

α-Diversity

We assessed α-diversity by measuring the total number of observed ASVs (species richness), Simpson evenness (species evenness) and Faith’s phylogenetic diversity (phylogenetic richness) of each snail hemolymph and hepatopancreas sample.

Species richness

We first assessed how much variation in species richness was explained by each factor using a fixed effect model (observed ASVs ~ cohort + day + sample type). Sample type (snail tissues [hemolymph, hepatopancreas], water) explained most of the variation in microbiome richness, followed by day (Table 2). Next, we explored which factors specifically shaped species richness in each snail tissue sample type using a second fixed effect model (observed ASVs ~ cohort + day). Day was the only significant factor for both hemolymph and hepatopancreas, while the cohort factor (exposed to S. mansoni parasites or control) had no significant effect (Table 2). Table 2. Factors influencing the α-diversityObserved ASVsSimpson evennessPhylogenetic diversityF valueP valueF valueP valueF valueP valueAll sample types Cohort1.420.231.010.322.010.16 Sample type78.42 < 2.2e-1610.534.2E-05110.79 < 2.2e-16 Day8.354.1E-090.740.647.061.1E-07Hemolymph Cohort2.800.102.240.143.110.08 Day3.661.5E-031.010.432.849.7E-03Hepatopancreas Cohort1.910.173.040.081.620.21 Day2.630.011.370.232.370.03We tested the impact of several factors (cohort, sample type and day) using linear mixed-effects models followed by an ANOVA on the different α-diversity metrics: observed ASVs, Simpson evenness and Faith’s phylogenetic diversity. The additive models included all three sample types (hemolymph, hepatopancreas, environmental water) or were run individually for hemolymph and hepatopancreas samples. Factors with P < 0.05 show a significant influence on the metrics tested. The degree of influence is highlighted by the F value, a high F value corresponding to a greater effect of the factor

We only observed a single timepoint with a significant difference in the number of ASVs between the hemolymph of control and exposed snails, on D21 (Wilcoxon test: w = 55, P = 0.018), with control snails having a higher number of ASVs (Fig. 3A). Water showed extensive variation in observed ASVs over time (range: 92–489, mean: 233.406); however, this variation did not appear to affect the number of ASVs in the snail hemolymph or hepatopancreas. The species richness index showed a significantly higher number of ASVs in the hemolymph than in the hepatopancreas and environmental water (Fig. 3A). We found approximately two times more ASVs in the hemolymph (262.407 ± 8.422) than in the hepatopancreas (142.564 ± 5.341) samples. The average number of observed ASVs for water samples (233.406 ± 19.283) was comparable to that of hemolymph; however, the number of observed ASVs in water showed dramatic variations over time compared to the hemolymph samples (SuppTable2).Fig. 3. Longitudinal comparison of α-diversity metrics between control and schistosome-exposed snails. Boxplots showing values for A observed ASVs, B Simpson evenness and C Faith’s phylogenetic diversity between control (green) and exposed (red) snail cohorts at each of the eight timepoints, and for each sample type (Hm: hemolymph, Hp: hepatopancreas, W: water). The longitudinal analysis of the observed ASVs and the Faith’s PD metrics reveals specific timepoints where schistosome infection significantly impacts the microbial composition of hemolymph (D21 and D39) or hepatopancreas (D28 and D39). The dashed line marks the period before snails were exposed to miracidia (D0) and after snails have been exposed (D7 onwards)

Species evenness

We measured the distribution of ASVs using the Simpson evenness index (Fig. 3B). The sample type factor explained all the observed variation in species evenness, based on the best fitting model for this index (evenness ~ cohort + day + sample type) (Table 2).

This community evenness was low across all samples (Fig. 3B), exhibiting a wider range in the hemolymph (range: 0.014–0.228) than in the hepatopancreas (range: 0.023–0.19) (Fligner-Killeen test: X^2^ = 12.914, df = 1, P = 0.0003), as previously observed [19]. In the hemolymph, we observed a significant decrease in species evenness only in exposed snails compared to control on day 39 (Wilcoxon test: W = 33, P.adj = 0.02) (SuppTable3). In the hepatopancreas, significant differences in evenness between cohorts were observed on day 28 and 39 (D28: Wilcoxon test: W = 2, P.adj = 0.013; D39: Wilcoxon test: W = 0, P.adj = 0.005).

Phylogenetic diversity

We estimated the phylogenetic differences between microbial communities using Faith’s phylogenetic diversity index. Using the same fixed model as for species richness and evenness, we found that the sample type factor explained most of the variation observed in the phylogenetic diversity (Table 2). When we applied the model to each snail tissue separately, the day factor significantly explained the variation observed in both the hemolymph and hepatopancreas samples.

Similar to species richness, only one timepoint—day 21 in the hemolymph—showed significant differences between cohorts (Wilcoxon test: W = 57, P.adj = 0.01) (SuppTable3). At day 21, control snails had significantly higher microbial phylogenetic diversity (32.381 ± 1.154) in their hemolymph than exposed snails (29.562 ± 1.238) (Fig. 3C). The phylogenetic diversity between tissues (Fig. 3C) mirrored the observed ASVs diversity, with water samples showing similar diversity levels to hemolymph and hepatopancreas samples exhibiting the lowest diversity (SuppTable3).

β-Diversity

We conducted the β-diversity analysis on hemolymph and hepatopancreas separately because of the strong effect of the sample type factor observed in the α-diversity analysis. To ensure that differences observed by PERMANOVA were due to group distance rather than group dispersal, we tested for the homogeneity of variances using a β-dispersion analysis. Nearly all of the hemolymph and hepatopancreas samples showed similar group dispersion between cohorts, except for three hemolymph comparisons and a maximum of four hepatopancreas comparisons (SuppFigure4).

We then explored differences between sample groups using Bray-Curtis and weighted UniFrac distances. Both metrics quantify dissimilarity between samples, with Bray-Curtis accounting for species abundance and weighted UniFrac also considering the phylogenetic relatedness of ASVs. We tested the effect of the factors cohort and day on the hemolymph groups using PERMANOVA analysis (distance ~ cohort + day) (Table 3). Both cohort and day significantly explained the variation observed with the Bray-Curtis distance (cohort: PERMANOVA: F = 3.9705, P = 0.0009; day: PERMANOVA: F = 4.6403, P = 0.0009) (Fig. 4A, SuppFigure5A). In contrast, neither factor cohort nor day explained the overall variation observed with weighted UniFrac distance (cohort: PERMANOVA: F = 0.5674, P = 0.5215; day: PERMANOVA: F = 1.6243, P = 0.0959) (Fig. 4B, SuppFigure5B). These results suggest that parasite infection may alter the relative abundance of certain taxa without significantly affecting the broader evolutionary composition of the microbiome. Additionally, we found that the hemolymph microbiome composition was significantly different between infected and uninfected snails, and within the infected cohort at days 7, 21, 37 and 39 with Bray-Curtis and days 7 and 21 with weighted UniFrac (Fig. 4). This suggests specific microbiome changes occurring at defined developmental stages of the schistosome parasite within its snail host. Table 3. Factors influencing the β-diversitydfSum of squaresR^2^F**P valueBray distance ~ cohort + day Hemolymph  Cohort10.880.033.979.99E-04  Day77.190.244.649.99E-04  Residual9921.900.73 Hepatopancreas  Cohort10.310.011.950.04  Day74.420.213.939.99E-04  Residual9915.890.77Weighted UniFrac distance ~ cohort + day Hemolymph  Cohort14.28E-035.10 E-030.570.52  Day70.090.101.620.10  Residual990.750.89 Hepatopancreas  Cohort13.64E-030.011.300.26  Day70.040.121.940.01  Residual990.280.87We tested the impact of cohort and day factors using PERMANOVA on additive models for the Bray-Curtis and weighted UniFrac β-diversity metrics. Models were run separately for each snail sample type (hemolymph and hepatopancreas). Factors with P < 0.05 had a significant influence on the metrics tested. The degree of influence is highlighted by the F value, a high F value corresponding to a greater effect of the factor. df degree of freedomFig. 4Longitudinal comparison of β-diversity metrics between control and schistosome-exposed snails for hemolymph and hepatopancreas. Heatmaps depicting the statistical significance between microbiome sampled on different days determined by using a pairwise post hoc test on the PERMANOVA models for A, B hemolymph (Hm) and C, D hepatopancreas (Hp), using the A, C Bray-Curtis and B, D weighted UniFrac metrics. The lower triangle shows pairwise F values, the upper triangle shows P values, and the days compared are marked on the axes. The hemolymph microbiome varied significantly between infected and uninfected snails, especially at D7 and D21, showing stage-specific changes during schistosome development. The hepatopancreas microbiome showed minimal impact from infection

We performed the same analysis with the hepatopancreas samples. Similarly, both cohort and day factors explained the variation observed with the Bray-Curtis distance (cohort: PERMANOVA: F = 1.95, P = 0.0389; day: PERMANOVA: F = 3.9336, P = 0.0009) (Fig. 4C, SuppFigure5C). However, only the day factor explained the observed dissimilarities with weighted UniFrac distance (cohort: PERMANOVA: F = 1.3002, P = 0.2587; day: PERMANOVA: F = 1.9374, P = 0.01099) (Fig. 4D, SuppFigure5D). While we observed some significant differences between days in the hepatopancreas microbiome, there were no consistent patterns at these timepoints. Additionally, we found no significant effect of the cohort in the microbiome composition in this tissue. These results suggest that schistosome infection does not significantly impact the microbiome composition of the hepatopancreas.

Taxonomy

We investigated the impact of S. mansoni infection on the microbiome community structure of B. glabrata snail hemolymph and hepatopancreas over time in both uninfected and infected snails. The microbiomes in snail samples (hemolymph and hepatopancreas) and environmental (water) samples were dominated by Proteobacteria, Bacteroidetes, Tenericutes and unassigned phyla (Fig. 5). While the microbiome composition of snail hemolymph and hepatopancreas was consistent with previous findings [19], the taxonomic composition for the water samples was unexpected. In our previous studies [19, 20], water samples were dominated by Actinobacteria, Bacteroidetes and Proteobacteria, whereas Actinobacteria was almost absent in the current study. This difference may be attributed to the water source and/or the frequency of the water change. In the present study, snails were kept in well water changed daily, whereas snails from previous studies were sampled from tanks containing several-month-old well water [19, 20].Fig. 5. Longitudinal taxonomic diversity between control and schistosome-exposed snails. The relative abundance of the top 20 most dominant phyla is shown for both control (left) and exposed (right) cohorts at each of the eight timepoints, including snail tissues (hemolymph and hepatopancreas) and environmental water samples. Samples were merged by sample type and day: each bar represents two water samples and 6–8 snail samples for hepatopancreas or hemolymph. The taxonomic compositions of both control and exposed snails are relatively stable over time, dominated by phyla Proteobacteria, Bacteroidetes and Tenericutes. We observed significant differences in the microbiome composition between snail tissues (hepatopancreas and hemolymph) and the water environment

As observed with Bray-Curtis distances, significant differences in microbiome compositions were evident between each sample type. While microbiome compositions remained relatively consistent over time, the proportion of the different phyla exhibited variations (SuppFigure 6). For instance, the proportion of Fusobacteria dramatically decreases in both hemolymph and water after the start of the experiment, for both controls and exposed cohorts. We observed a similar decrease for Actinobacteria in the hepatopancreas. In contrast, Acidobacteria proportions increased in the hemolymph, particularly in the control cohorts starting at day 21. A similar increase was observed for Nitrospirae and Cyanobacteria at day 35. Infected snails also showed a dramatic increase in Firmicute on day 39, in both the hemolymph and hepatopancreas.

Discussion

Limited and specific effect of schistosome on snail hemolymph microbiome during infection

We investigated the impact of schistosome parasite infection on the microbiomes of both the hemolymph and hepatopancreas in Biomphalaria snails (Fig. 1). We selected these two sample types for three reasons: (i) both harbor diverse and specific microbiomes [19], (ii) the hemolymph bathes the snail’s organs and parasites and has the greatest microbial diversity [19, 20] and (iii) the hepatopancreas is the primary organ where schistosome parasite sporocysts reside and develop [43]. Comparisons of microbiomes between control and exposed snails during the pre-patent (before cercariae larvae are produced and shed) and patent periods (when cercariae are produced and released) revealed no significant overall impact of parasitism. Taken together, the α- and β-diversity results indicate that the primary difference observed was between sample type (hemolymph vs. hepatopancreas), consistent with our previous observations in uninfected snails [19]. This difference in microbial composition of the snail tissues persisted over time, regardless of infection status. We also observed an effect of time (day factor), but to a lesser extent.

While infection status had no overall impact on the host microbiome, we observed differences at very specific timepoints. During the pre-patent period, the hemolymph microbiome of infected snails showed differences in composition at day 7 and 21. Day 21 is a critical timepoint in schistosome development: by this stage, schistosome sporocysts have migrated through host tissues and are actively colonizing and multiplying within the snail’s hepatopancreas [43]. During the patent period, we observed differences in community evenness at day 28 for the hepatopancreas and at day 39 for both the hemolymph and hepatopancreas. Heavy colonization of the snail hepatopancreas and ovotestis [43, 44] and tissue damage due to cercarial shedding could explain these differences. Previous studies using schistosome parasites have shown that the microbiomes of the mammalian [45, 46] and the invertebrate hosts [32] can be marginally impacted by infection. Portet et al. reported no significant changes in α-diversity of the whole snail microbiome following infections of B. glabrata snails (BgBRE) with S. mansoni (SmBRE and SmVEN) on days 1, 4 and 25 after primary infection, except for day 4 in BgBRE infected with SmBRE [32]. Our results, focusing on organ and hemolymph microbiomes, indicate similar findings at the beginning of the pre-patent period but suggest specific changes toward the end of this period. This difference could be attributed to the use of SmLE parasites, which are known to have a stronger impact on the snail host physiology [36], or it may have been revealed by the analysis of tissue-specific microbiomes that were otherwise masked by the use of whole snails [32].

We observed no significant differences in bacterial density between control and exposed snails at any timepoint, or within cohorts over time, suggesting that bacterial density within the host hemolymph remains relatively stable. The minimal change in bacterial density may reflect the constant environment conditions within snails, or tight control bacterial density by the host snail, and is unaffected by schistosome infection.

Our results contrast with those obtained from bivalves infected with pathogens [9], where mollusk microbiomes are significantly disturbed. However, those studies involved transient pathogens that do not establish long-term relationships with their hosts and do not rely on host fitness for transmission. In cases of co-dependence, the impact on microbiomes can differ, although examples of this “strategy” in invertebrates are scarce. The impact of the parasitic snail Coralliophila violacea on the microbiome of its coral host, Porites cylindrica, offers an interesting parallel [15]. This sedentary snail employs a feeding strategy that minimizes tissue damage for its coral host, disturbing the coral microbiome only at the feeding site. In contrast, other mobile predatory snails tend to disturb the coral microbiome on a larger scale. The limited effect of parasites on host microbiomes could be under strong natural selection, as parasite fitness is closely dependent on host fitness.

Stability of the microbiomes over time

The microbiomes of the snail hemolymph and hepatopancreas remained stable across most timepoints. The weighted UniFrac distances showed only minor inconsistencies between hemolymph or hepatopancreas over time, with no specific patterns observed. In contrast, differences became apparent when using Bray-Curtis distances. Unlike Bray-Curtis, the UniFrac distance weights microbial abundance by the phylogenetic distance between microbiome community members (i.e. the greater the phylogenetic distance, the greater the weight). Therefore, the differences detected with Bray-Curtis distances but not with weighted UniFrac strongly suggest that changes occurred primarily among closely related taxa. Additionally, microbiome diversity over time did not appear to be affected by the environmental variations: the high variability in microbial diversity in the water did not impact the diversity of the snail hemolymph or hepatopancreas microbiomes.

The apparent stability of host microbiomes may be attributed to survivorship bias, but it could also have a biological origin. Parasite infection negatively impacts the survival of snails (SuppFigure1), meaning that only the surviving snails are sampled. However, the survival of the infected snails during the pre-patent period was similar to that of the control snails and decreased sharply only during the patent period. Therefore, any observation made during the pre-patent period should only be minimally impacted by this bias. Alternatively, the minimal disturbance in the host microbiome could result from parasite action or microbiome resilience. Schistosomes and their snail hosts have co-evolved for over 200 million years [23], and their fates are intimately linked during a relatively long developmental period of several weeks. Parasites have likely been under strong selective pressure to reduce their impact on host microbiomes to ensure their fitness and transmission. They may actively control their host microbiome using mechanisms such as extracellular vesicles [47, 48] or may passively avoid disturbing the host microbiome by evading or controlling the immune system [49, 50]. Interestingly, the immune status of infected snails during schistosome development, especially during the patent period, has received little attention. For example, phenoloxidase activity (a component of the humoral immune response in invertebrates) does not differ between infected and non-infected snails before week 7 (D49) post-infection, suggesting minimal impact of the parasite during the early stages of the patent period [26]. The minimal disturbance we observed could also be independent of the interaction between the host and parasite, instead explained by the resilience of the microbial community, which contributes to maintaining homeostasis within the host [51]. The snail tissue microbiomes are diverse [19], which can enhance the stability of the microbiome when challenged by stressful conditions [52]. Similarly, in oysters, which also exhibit diverse and distinct tissue microbiomes, changes in the hemolymph microbiome due to environmental stressors are limited to specific taxa [53].

What would be the consequences of manipulating the snail host microbiome for schistosome parasites?

If parasites do not disturb their host microbiomes to ensure their fitness and transmission, what would be the consequence of manipulating the host microbiome on the parasites? A limited number of studies have explored this topic. For instance, snails first treated with antibiotics and then exposed to schistosomes showed a shorter pre-patent period but produced fewer cercariae [54]. Antibiotic treatment of schistosome-infected snails led to the suppression of cercarial production [55]. While these results may suggest a complex interaction between the snail microbiome and the schistosome parasite, it is challenging to discern whether the impact on schistosome development was due to changes in the microbiome or the effects of the antibiotics themselves. The use of axenic (germ-free) snails, generated without the administration of antibiotics [3], could help elucidate the role of the host microbiome in the interactions between snails and schistosomes. Absence of snail microbiome, which does not prevent snail infection [56], may affect schistosome fitness. The addition of a synthetic microbial community to axenic snails prior to or after schistosome infection may also reveal specific interactions.

Limitations

Our study has several limitations that should be addressed in future studies to better understand the impact of schistosome infection on snail host microbiomes. First, our sample sizes per group, with 6–8 snails per cohort at each timepoint, were relatively small and may have limited our ability to detect subtle effects of parasitism on the host microbiomes. Increasing the group sample sizes will improve the statistical power to detect differences in the microbiome. Second, we sampled different individuals from a snail population over time rather than resampling the same individuals. While we selected individuals from an inbred snail population [35, 36], we cannot exclude the potential effect of snail genotypes on our measurements, which could have masked small effects of the parasites on the host microbiomes as well. Development of a technique that allows resampling of tissues (e.g. hemolymph) from individual snails would be critical. Third, we used a single host-parasite combination in these experiments: further replication with different snail-schistosome combinations is needed to test the robustness of our findings and their possible generalization. Fourth, a recent study has shown that miracidial stage may harbor a microbiome [57]. The parasite microbiome was not investigated in our study. However, it is unlikely that the presence of this parasite microbiome could significantly impact the microbiomes of the host, as we have shown no large impact of schistosome infection on the microbiome of the snail tissues in the pre-patent period. Fifth, the snail environment was quite controlled: water and lettuce were changed daily, and lettuce was the only source of food. Clean water may also have limited the exchange of microbes between individuals. Different rearing conditions may reveal varying impacts of parasites on the microbiomes of their hosts. Finally, sequencing the V4 of the 16S gene may provide only a partial view of microbiome changes during infection, as it cannot resolve differences at the microbial strain level or detect alterations in microbiome metabolic pathways that may significantly impact the host. Full-length 16S sequencing and metagenomics analyses using long-read sequencing will be required to obtain a more comprehensive understanding of the microbiome dynamics.

Conclusions

Overall, schistosome infections have a limited impact on the microbiome of the hemolymph and hepatopancreas of their snail hosts. Bacterial density and microbiome diversity do not significantly differ between control and infected snails, suggesting either maintenance of microbiome homeostasis by the snail and/or the parasite during infection, or a lack of interaction between the parasite and snail microbiomes. Differences in microbiome composition between infected and uninfected snails were observed only at specific timepoints, potentially corresponding to key stages in parasite development. The stability of the microbiome during the prepatent period may reflect selective pressure on the parasite to avoid disturbing the host microbiome, thereby ensuring the survival of both host and parasite. This microbiome stability warrants further investigation into the underlying mechanisms and the consequences of microbiome perturbations on host-parasite interactions.

Supplementary Information

Supplementary material 1: Table 1. Average number of reads.Supplementary material 2: Table 2. Mean α-diversity over time.Supplementary material 3: Table 3. α-diversity statistics.Supplementary material 4: Figure 1. Survival curve.Supplementary material 5: Figure 2. Average reads.Supplementary material 6: Figure 3. Rarefaction curves.Supplementary material 7: Figure 4. β-diversity: homogeneity of variances.Supplementary material 8: Figure 5. β-diversity: additive models day comparisons.Supplementary material 9: Figure 6. Abundance of the phyla showing the greatest variation during the course of infection.