Introduction
Cercarial dermatitis, also known as ‘swimmer's itch’ (SI), is a rash caused by avian schistosomes, snail-borne parasites related to the causative agents of human schistosomiasis (Schistosoma spp.; Brant and Loker, Reference Brant and Loker2009; Colley et al., Reference Colley, Bustinduy, Secor and King2014; Horák et al., Reference Horák, Mikeš, Lichtenbergová, Skála, Soldánová and Brant2015). Avian schistosomes are a diverse group of trematode (flatworm) parasites, mostly in the genus Trichobilharzia, that normally use water birds as their definitive hosts (Brant and Loker, Reference Brant and Loker2009). Infected snails release infectious free-living larvae, known as cercariae, which swim through the water in search of a bird to infect. Avian schistosomes cannot complete their life cycles in humans, but their cercariae sometimes mistake humans for birds and penetrate our skin resulting in the distinctive rash (Verbrugge et al., Reference Verbrugge, Rainey, Reimink and Blankespoor2004a; Horák et al., Reference Horák, Mikeš, Lichtenbergová, Skála, Soldánová and Brant2015). Anecdotal evidence suggests that SI is increasing in Michigan inland lakes, where it discourages swimming and other recreational activities impacting the local economy during the summer (Muzzall et al., Reference Muzzall, Burton, Snider and Coady2003; Verbrugge et al., Reference Verbrugge, Rainey, Reimink and Blankespoor2004b; McPhail et al., Reference McPhail, Froelich, Reimink and Hanington2022; Soper et al., Reference Soper, Raffel, Sckrabulis, Froelich, McPhail, Ostrowski, Reimink, Romano, Rudko and Hanington2023). The drivers for distribution of human schistosomes have been widely studied at multiple spatial scales (Yang et al., Reference Yang, Gemperli, Vounatsou, Tanner, Zhou and Utzinger2006; Brooker, Reference Brooker2007; Zhou et al., Reference Zhou, Yang, Yihuo, Liu, Wang, Wei and Jiang2011); however, ecological drivers of among-lake variation in avian schistosomes have been relatively understudied (McMullen and Brackett, Reference McMullen and Brackett1948; Rudko et al., Reference Rudko, Reimink, Froelich, Gordy, Blankespoor and Hanington2018; Soper et al., Reference Soper, Raffel, Sckrabulis, Froelich, McPhail, Ostrowski, Reimink, Romano, Rudko and Hanington2023).
Prior studies on schistosomes and other trematode parasites show that their distributions closely match the distributions of their molluscan intermediate hosts (Skírnisson et al., Reference Skírnisson, Galaktionov and Kozminsky2004; Dida et al., Reference Dida, Gelder, Anyona, Matano, Abuom, Adoka, Ouma, Kanangire, Owuor and Ofulla2014; Marszewska et al., Reference Marszewska, Cichy, Heese and Żbikowska2016; Gordy et al., Reference Gordy, Cobb and Hanington2018; Soper et al., Reference Soper, Raffel, Sckrabulis, Froelich, McPhail, Ostrowski, Reimink, Romano, Rudko and Hanington2023). Therefore, environmental factors that influence the abundance of intermediate host snails are likely to be major predictors of cercaria abundance (see Fig. 1. for a priori hypothesized relationships; Rohr et al., Reference Rohr, Raffel, Sessions and Hudson2008a; Paull et al., Reference Paull, Raffel, LaFonte and Johnson2015). Studies on other trematode parasites have shown that pollution from agricultural chemicals including fertilizer, herbicides and insecticides can all increase the densities of snail intermediate hosts, leading to higher infection rates (Rohr et al., Reference Rohr, Raffel, Sessions and Hudson2008a, Reference Rohr, Schotthoefer, Raffel, Carrick, Halstead, Hoverman, Johnson, Johnson, Lieske, Piwoni, Schoff and Beasley2008b; Halstead et al., Reference Halstead, McMahon, Johnson, Raffel, Romansic, Crumrine and Rohr2014). These chemicals typically act by either increasing growth of attached algae (periphyton) eaten by snails by clearing the water of phytoplankton, increasing water clarity and light penetration (Johnson and Chase, Reference Johnson and Chase2004; Brown and Lydeard, Reference Brown, Lydeard, Thorp and Covich2009) or killing off snail predators (e.g. crayfish; Rohr et al., Reference Rohr, Raffel, Sessions and Hudson2008a; Halstead et al., Reference Halstead, McMahon, Johnson, Raffel, Romansic, Crumrine and Rohr2014; Halstead et al., Reference Halstead, Hoover, Arakala, Civitello, De Leo, Gambhir, Johnson, Jouanard, Loerns, McMahon, Ndione, Nguyen, Raffel, Remais, Riveau, Sokolow and Rohr2018). However, water clarity might also be altered by invasive zebra or quagga mussels (Dreissena spp.) that filter algae out of the water column (Geisler et al., Reference Geisler, Rennie, Gillis and Higgins2016). Other potentially important factors for snail abundance include habitat characteristics such as wave action, water depth, lake size, abundance of submerged vegetation, abundance of snail predators (crayfish) or the availability of solid substrates (Laman et al., Reference Laman, Boss and Blankespoor1984; Dida et al., Reference Dida, Gelder, Anyona, Matano, Abuom, Adoka, Ouma, Kanangire, Owuor and Ofulla2014; Rohr et al., Reference Rohr, Sack, Bakhoum, Barrett, Lopez-Carr, Chamberlin, Civitello, Diatta, Doruska, De Leo, Haggerty, Jones, Jouanard, Lund, Ly, Ndione, Remais, Riveau, Schacht, Seck, Senghor, Sokolow and Wolfe2023). For example, one of the strongest positive predictors of host snails for human schistosomes is the abundance of aquatic vegetation, likely due to the increased surface area available for periphyton growth (Underwood et al., Reference Underwood, Thomas and Baker1992; Wood et al., Reference Wood, Sokolow, Jones, Chamberlin, Lafferty, Kuris, Jocque, Hopkins, Adams, Buck, Lund, Garcia-Vedrenne, Fiorenza, Rohr, Allan, Webster, Rabone, Webster, Bandagny, Ndione, Senghor, Schacht, Jouanard, Riveau and De Leo2019; Rohr et al., Reference Rohr, Sack, Bakhoum, Barrett, Lopez-Carr, Chamberlin, Civitello, Diatta, Doruska, De Leo, Haggerty, Jones, Jouanard, Lund, Ly, Ndione, Remais, Riveau, Schacht, Seck, Senghor, Sokolow and Wolfe2023). Based on the results of these past studies, we predicted greater abundance of host snails (namely Lymnaea, Planorbella and Physa spp. snails; Blankespoor and Reimink, Reference Blankespoor and Reimink1991; Muzzall et al., Reference Muzzall, Burton, Snider and Coady2003; McPhail et al., Reference McPhail, Rudko, Turnbull, Gordy, Reimink, Clyde, Froelich, Brant and Hanington2021), and thus avian schistosome cercariae (namely Trichobilharzia spp.; McPhail et al., Reference McPhail, Rudko, Turnbull, Gordy, Reimink, Clyde, Froelich, Brant and Hanington2021; Rudko et al., Reference Rudko, McPhail, Reimink, Froelich, Turnbull and Hanington2022; Soper et al., Reference Soper, Raffel, Sckrabulis, Froelich, McPhail, Ostrowski, Reimink, Romano, Rudko and Hanington2023), in sites with higher water clarity (the inverse of turbidity), faster growth of attached algae, more submerged vegetation, higher nutrient loading in water or sediment and agricultural or urban land use at the local or watershed scales (Fig. 1).
A priori hypothesized drivers of snail and avian schistosome abundance in northern Michigan lakes. Solid lines represent positive effects, and dashed lines represent negative effects.
However, cercaria abundance at any given site is not solely driven by host snail abundance. Temperature has long been known to have direct positive effects on rates of parasite development and cercaria production from snails for a variety of trematode species (McCreesh and Booth, Reference McCreesh and Booth2014; Paull et al., Reference Paull, Raffel, LaFonte and Johnson2015; Nguyen et al., Reference Nguyen, Boersch-Supan, Hartman, Mendiola, Harwood, Civitello and Rohr2021). Furthermore, snails typically have higher rates of development, reproduction and population growth at warmer temperatures (Nguyen et al., Reference Nguyen, Boersch-Supan, Hartman, Mendiola, Harwood, Civitello and Rohr2021). However, it is important to recognize that temperature can have complex and non-linear effects on snail and trematode biology depending on the temporal scale investigated (e.g. thermal stress; Paull et al., Reference Paull, Raffel, LaFonte and Johnson2015). Trematode-infected snails also have higher rates of cercaria production when they have access to higher quality food sources (Civitello et al., Reference Civitello, Fatima, Johnson, Nisbet and Rohr2018). Cercaria abundance may also be driven by definitive bird host visitation to the lake, which should increase infection prevalence in populations of intermediate host snails (Muzzall et al., Reference Muzzall, Burton, Snider and Coady2003; Byers et al., Reference Byers, Blakeslee, Linder, Cooper and Maguire2008). All of these factors may lead to more cercaria being produced and released into the water (Fig. 1). Lastly, cercaria abundance at a particular location could be impacted by factors that influence cercaria distributions within a lake, such as water currents and related variables like wind, wave action, effective fetch (the distance over which wind can travel across open water as a measure of potential wave action and surface currents), wind speed and wind direction (Fig. 1). For example, onshore wind and water currents have been found to increase cercaria abundance at local beaches, presumably by bringing in cercariae that were released off-site into the local environment (Upatham, Reference Upatham1974; Muzzall et al., Reference Muzzall, Burton, Snider and Coady2003; Rudko et al., Reference Rudko, Reimink, Froelich, Gordy, Blankespoor and Hanington2018; Sckrabulis et al., Reference Sckrabulis, Flory and Raffel2020).
We conduct an exploratory analysis investigating potential drivers of avian schistosome abundance in northern Michigan inland lakes by collaborating with a large network of local volunteers to survey 16 lakes throughout northern MI in the summer of 2016 (Fig. 2; Table S1). We focused our study on summer months, at the time of peak human recreational water use (i.e. swimming) and reported SI incidence in this region (Sckrabulis et al., Reference Sckrabulis, Flory and Raffel2020). We also quantified several dozen environmental predictors hypothesized a priori to influence cercaria and host snail abundance, based on known relationships between various abiotic and biotic factors in wetland ecosystems (Fig. 1; Table 1). Here, we take an analytical approach similar to Rohr et al. (Reference Rohr, Schotthoefer, Raffel, Carrick, Halstead, Hoverman, Johnson, Johnson, Lieske, Piwoni, Schoff and Beasley2008b), where we measured a large number of environmental factors that could potentially influence our system but restricted models to only include plausible predictors of individual response variables. We used stepwise model selection with multiple linear regression and generalized linear mixed effects models to conduct an exploratory analysis of hypothesized predictors of cercaria abundance and host snail abundance.
Map of survey sites in northern Michigan. Panel A depicts all sites at the state level, and panel B depicts all sites as an inset.
Variables included in between-site analyses of avian schistosome abundance and the best predictors of that abundance
a Predictor variable that was quantified at the lake level instead of separately for each site.
Materials and methods
Sampling sites and times
We sampled at 38 sites on 16 inland lakes in northwestern Michigan (Fig. 1, Table S1). Site locations were largely determined by stakeholder interest due to each lake individually funding the study on their lake. Despite this, these sites represented a range of lake sizes, shoreline types (e.g. beach vs marsh) and levels of human activity. Each sampling site was restricted to a 15 m stretch of shoreline and extended at least 15 m into the littoral zone (out to approximately 0.8 m depth) and 15 m into the riparian zone. We conducted all aquatic sampling within the limits of these zones. Before daily sampling began, 2 HOBO pendant temperature loggers (Onset Computer Corp., Bourne, MA, USA) and a marker buoy were installed at each site at 60 cm depth. Both loggers recorded hourly temperature readings, but 1 logger at each site also recorded hourly light intensity and was anchored to ensure the light sensor remained horizontal and facing upwards. We also placed periphyton samplers at this depth (see Supplemental Materials ‘Field survey methods – Periphyton growth’). One exception was Deer Lake, where the slope was too shallow that we decided to place equipment at only 45 cm depth to avoid going too far outside the 15 × 15 m sampling area. Daily community science sampling occurred during a 28-day period in July and August 2016, starting on July 5 (detailed below). Daily measurements started 2–3 weeks earlier at selected sites on Crystal Lake and Intermediate Lake based on individual lake support. Site visits and snail surveys started the week prior to this sampling period (the week of June 26) and occurred every 2 weeks at most sites (3 total surveys) and weekly at a few selected sites (at least 1 site at each lake; 5 total surveys).
Community science sampling – daily cercaria samples, weather and bird observations
Cercaria release into the environment is highly variable at a daily time scale for avian schistosomes (Soldánová et al., Reference Soldánová, Born-Torrijos, Kristoffersen, Knudsen, Amundsen and Scholz2022). Therefore, to capture a clear pattern of cercaria abundance at each site, we employed daily sampling utilizing a community science approach. Daily cercaria samples were collected by trained community science volunteers, who were provided with pre-labelled 15 mL centrifuge tubes, a 1 L pitcher, a custom-made 35 μm Nitex mesh filter (Fig. S1), a squirt bottle filled with tap water, another squirt bottle filled with 70% ethanol, a small funnel, and a stand to hold the filter during filtration (if necessary). Cercariae are likely to have aggregated distributions in the water, so to reduce the chances of missing small clumps of cercariae we sought to collect a distributed sample from each sampling site. We therefore asked volunteers to filter 24 separate 1-L scoops of water for each daily sample at each sampling site. Each morning, when Trichobilharzia spp. cercariae are released into the environment (~8:00–10:00 AM; Soldánová et al., Reference Soldánová, Selbach and Sures2016), volunteers entered the water near the left boundary of the site while holding the 1-liter pitcher and the Nitex mesh filter, haphazardly collecting water from the surface using the pitcher and pouring it through the filter as they moved in a zig-zag pattern throughout the littoral zone of the sampling site. Volunteers were allowed to collect their water samples in a bucket prior to filtration, allowing sediment to settle to the bottom of the bucket prior to decanting the water through the mesh filter. Avian schistosome cercariae typically move toward light and stay near the water surface (Cort and Talbot, Reference Cort and Talbot1936; Brachs and Haas, Reference Brachs and Haas2008), so they are unlikely to have been lost through this procedure. Once the lake water had been filtered, the squirt bottle containing tap water was used to rinse any residual material into 1 corner of the filter, followed by a final rinse with 70% ethanol to wash the sample into a 15 mL centrifuge tube using a small funnel. The use of 70% ethanol killed and preserved the DNA of any cercariae captured in the sample. Volunteers rinsed and back-rinsed their filters daily with tap water to keep the mesh clean. Samples were stored at room temperature until they were transported back to Oakland University at the end of the survey for further processing. Community science volunteers at each site were also encouraged to collect daily data describing weather conditions, air temperature, onshore wind speed and direction with the provided equipment. They were also asked to tally any sightings of birds within visual range of the study site.
qPCR quantification of avian schistosomes
After the survey was completed, we used qPCR to quantify the number of gene copies of avian schistosome DNA. Due to budgetary limitations, we developed a low-cost protocol for this process. Our protocol is like that of Rudko et al. (Reference Rudko, Reimink, Froelich, Gordy, Blankespoor and Hanington2018); however, we aimed to develop this protocol to be economical and scalable to a larger number of volunteer-collected samples. Briefly, we used a recently detailed protocol and custom primers that target a highly conserved region of 18S ribosomal RNA gene sequences specific to schistosome parasites (Jothikumar et al., Reference Jothikumar, Mull, Brant, Loker, Collinson, Secor and Hill2015). This assay only differentiates schistosomes vs non-schistosomes; we assumed that any amplified DNA largely belonged to the genus most associated with swimmer's itch in this geographic region, Trichobilharzia, though other SI-causing schistosomes from the genera Schistosomatium and Gigantobilharzia have also been observed in Michigan (Muzzall et al., Reference Muzzall, Burton, Snider and Coady2003; McPhail et al., Reference McPhail, Froelich, Reimink and Hanington2022; Rudko et al., Reference Rudko, McPhail, Reimink, Froelich, Turnbull and Hanington2022, Soper et al., Reference Soper, Raffel, Sckrabulis, Froelich, McPhail, Ostrowski, Reimink, Romano, Rudko and Hanington2023). We decided to divide each field sample in half prior to extraction, which we hoped would allow us to assess potential sources of error in the assay (i.e. extraction error vs qPCR error). For each field sample, we added together the qPCR estimates of cercaria abundance for the 2 extractions. qPCR data is typically lognormally distributed, so we log10-transformed the daily cercaria estimates and averaged them to obtain log10 cercaria abundance for each site. Log10 cercaria abundance was used as our proxy for avian schistosome risk for all our among-site analyses. For a full description of this protocol, see the Supplementary Materials (‘Low-cost protocol for avian schistosome DNA detection’).
Habitat assessment
At the beginning of the sampling period, a team of Raffel Lab student researchers conducted a habitat assessment in the littoral zone (15 × 15 m in water) and the riparian zone (15 × 15 m on land) using a standardized checklist modelled after the Environmental Protection Agency's 2012 National Lakes Assessment protocol (EPA 841-B-11-003). For all components of the site assessment, a numerical score was used to indicate abundance of landscape or substrate types based on the following numeric index: 0 = Absent, 1 = Sparse (<10% coverage), 2 = Moderate (10 − 40% coverage), 3 = Heavy (40 − 75% coverage) and 4 = Very Heavy (>75% coverage). The riparian ground cover was observed, checking for the abundance of vegetation (woody shrubs, saplings, herbs and grasses) and barren areas (dirt/sand/rock). Riparian canopy and human influence were also classified and recorded. Habitat in the littoral zone was assessed by noting the presence/absence and abundance of various substrates (boulder, cobble, gravel, sand, muck) and types of aquatic vegetation (submergent, emergent, floating and total plant cover).
Visual quadrat surveys
Each site was visited 3–5 times to measure the abundance and diversity of aquatic snails within the littoral zone, which we divided into 3 depth areas (0–40 cm, 40–80 cm and >80 cm). PVC quadrat sampling frames (0.09 m2) were haphazardly tossed within each depth zone until a total of 4 frames were counted in each. A clear-bottom viewing bucket was used to observe and count the number of snails on the substrate within the boundaries of each frame. Densities were recorded for snails, visually identified to the genus level in the field. Densities of any other organisms (e.g. crayfish and mussels) present in a given quadrat were also recorded as they were encountered. Snails were also collected from throughout the littoral area and preserved in 70% ethanol to verify the proportion of snails from each species during 15-minute intensive targeted searches after quadrat sampling was completed during each visit. We computed several indices of snail abundance, including total abundance, and both individual densities of each snail type, and combined densities of Lymnaea, Planorbella and Physa snails as potential predictors of cercaria abundance.
Environmental variables
Alongside the visual quadrat surveys and timed snail collections conducted at the time of each visit, we also collected samples to assess other biotic and abiotic variables that might impact snail or cercaria presence (Fig. 1). We set out several sampling apparatuses to collect mussel settling rates and periphyton/biofilm formation. We collected water and sediment samples to assess water and sediment chemistry. We conducted multiple instances of crayfish trapping and zooplankton sampling. Lastly, we also assessed lake watershed land-use information using GIS. See the Supplementary Materials (‘Field survey methods’) for a full description of methods used for each of these variables.
Statistical analyses
To investigate possible relationships between environmental factors (abiotic and biotic) and host snail cercaria abundance, we took an analytical approach similar to that of Rohr et al. (Reference Rohr, Schotthoefer, Raffel, Carrick, Halstead, Hoverman, Johnson, Johnson, Lieske, Piwoni, Schoff and Beasley2008b), where we measured a large number of variables that could plausibly influence our system. There are potential risks of obtaining spurious low P values due to model over-fitting (overly complex models) or the presence of influential outliers in analyses that contain many potential predictors relative to sample size (n = 38 sites). We sought to reduce these risks by (1) sampling a large number of sites relative to prior studies in this region, (2) limiting the analysis of each response variable to seemingly plausible explanatory variables (Fig. 1), (3) log10-transforming count variables and those with skewed distributions, (4) examining each statistical relationship graphically to ensure that it had a plausible directionality and was not driven by influential outliers and (5) limiting the size of individual models to a maximum of 5 predictors during model selection.
All data were aggregated to the site level (n = 38), using the average of each measured variable over the entire study period (e.g. average daily log10 cercaria abundance over at least 1 month, or average log10 snail density across at least 3 quadrat surveys). This resulted in all but a few predictor variables having complete data for all 38 sites (Table S3). We generated a correlation matrix to help identify pairwise relationships between variables, where relationships were considered strong candidates when the correlational statistic r was > 0.4 or < −0.4. A summary of all correlates with response variables of interest is presented in Table S2. We then used the program R (v4.3.1; R Core Team, Reference R Core Team2023) to conduct stepwise model selection to identify a set of likely predictor variables for each response variable (linear regression using the ‘regsubsets’ function in the package ‘leaps’; Lumley and Miller, Reference Lumley and Miller2020). A list of potential predictors included in the model for each response variable is provided in Table S3. To reduce chances of model over-fitting, we limited the maximum number of variables per ‘regsubsets’ model to 5 (nvmax = 5); because this function only returns models of the size specified, we sorted the list of models returned by the function by adjusted R 2 to assess model fit and predictor inclusion in the model. This generated a set of up to 5 potentially important predictor variables for each response variable. We conducted 2 versions of this analysis for each response variable that included temperature variables as hypothesized predictors, 1 with and 1 without including these predictors, due to missing temperature data from 1 site.
We then combined the top predictor variables identified as potentially important from either ‘regsubsets’ output (between 4 and 9 predictors) into a single model and used manual backward selection to remove non-significant predictors via F-tests (P < 0.05). We also considered predictors for possible removal if the effect direction was opposite of the hypothesized effect (see Fig. 2); these were examined on a case by case basis (see Discussion). At this stage, we added 1 additional hypothesized predictor to each relevant model that was not included in the ‘regsubsets’ output due to missing data from 1 site (sediment phosphorus; Table S3), assessed its significance, and again simplified the model via backward selection.
Although most predictor variables in this analysis were measured at the level of individual sites, some predictors had only a single value for each lake (e.g. maximum lake depth). We therefore sought to account for potential random effects of Lake in our analysis. The ‘regsubsets’ function only supports linear regression models (i.e. no random effects), and methods allowing for automated stepwise or all-subsets selection of mixed effects models are computationally intensive (e.g. the ‘dredge’ function in the package ‘MuMIn’; Bartón, Reference Bartón2019). We therefore added a random effect of ‘Lake’ to each final model, to guard against variable inclusion due to pseudoreplication (generalized linear mixed effects models using the ‘lmer’ function in the package ‘lme4’; Bates et al., Reference Bates, Maechler, Bolker and Walker2015). Once we had a final mixed effect model for each response variable (presented in Tables 2 & S4), we tested for spatial confoundment by adding main effects of latitude and longitude to each model, followed by another round of backward selection. If a predictor was kicked out by adding either latitude or longitude to the model, we considered it a spatially confounded predictor.
We verified the robustness of our core model outputs by also utilizing the ‘exhaustive’ method within the ‘regsubsets’ function to conduct an all-subsets analysis for each of the response variables presented in Table 2, again with nvmax = 5. We examined the top 10 models for each response variable in order of adjusted R 2. All-subsets model outputs can be difficult to interpret in analyses containing clusters of highly correlated predictor variables, but this method is less prone to variable ‘trapping’ than stepwise selection procedures. Here we used all-subsets to confirm that each of the primary predictors from our core models, or a closely related predictor, appeared in all or most of the top 10 models for each response variable.
Final models for each response variable of interest, following stepwise model selection
All final models included ‘Lake’ as a random effect. Note that the ‘Anova’ function from the ‘car’ package uses the Kenwood-Roger approximation for estimating degrees of freedom for F-tests, which can result in non-integer values (Fox and Weisberg, Reference Fox and Weisberg2019).
a Variable with missing data (one missing datapoint for sediment phosphorus)
b Lake-level variable
c Buildings became non-significant when Longitude was added to the final model (spatial autocorrelation)
In general, we are more confident of results where (1) the final model was relatively simple (i.e. 3 or fewer significant predictors), (2) predictors present in the top model were also among the top single correlates of a given response variable, (3) predictors selected by stepwise selection also appeared in the top 10 models in all-subsets analysis, (4) there was no evidence that a particular relationship was driven by spatial confoundment, and (5) there was an obvious mechanistic explanation for strength and direction of an observed statistical relationship.
Note that we also conducted an analysis of potential predictors of zebra mussel abundance (see Sckrabulis, Reference Sckrabulis2020), but this is outside the scope of the current manuscript.
Results
General patterns of snail and parasite distribution
Out of the 38 sites sampled in this study, 31 generated positive qPCR results for the presence of avian schistosome cercariae. There was a mean density of 0.05 log10 cercariae per liter across all sites. At 30 sites, we detected known host snails of avian schistosomes (15 with Lymnaea, 18 with Physa, 14 with Planorbella). At 4 sites, we detected avian schistosome DNA but did not detect known host snails of any genus.
Cercaria abundance
The strongest significant predictor of cercaria abundance was the density of Lymnaea spp. snails at each site (Log10 Lymnaea; Table 2; Figure 3A & B). Accounting for the presence of Physa or Planorbella spp. snails in any way did not increase predictive power of the model. We also found a significant negative effect of submerged aquatic vegetation on cercaria abundance (Table 2; Figure 3C & D). This effect was not driven by influential outliers and was evident with or without including Lymnaea spp. abundance in the model (Fig. 3C & D). Furthermore, we were unable to find any alternative variables we measured that correlated with vegetation and might help account for this pattern. After accounting for Lymnaea spp. abundance and submerged vegetation, we also detected a significant positive effect of sediment phosphorus (Table 2; Figure 3E & F). This pattern was completely dependent on the effects of the other 2 variables on the model, however, as there was no evidence of a direct correlation between cercaria abundance and sediment phosphorus prior to stepwise model selection (Table S2; Fig. S2). In all-subsets analysis, the top 10 models all contained Log10 Lymnaea density and sediment phosphorus (Table S6). All 10 top models also contained either submerged vegetation or the highly correlated (and biologically similar) ‘total vegetation’ (Table S6). Neither latitude nor longitude contributed significantly to the final model.
Residual plots for the main predictors of cercaria abundance (Panels A, C, E), Lymnaea density (Panels B & D), and submerged vegetation (Panel F). The top predictors of cercaria abundance were (A) Lymnaea spp. density, (C) submerged vegetation, and (E) sediment phosphorus. The top predictors of Lymnaea spp. density were (B) maximum lake depth and (D) deciduous tree cover. The top predictors for submerged vegetation were (F) sediment phosphorus and presence of buildings (not shown). All models depict each relationship after accounting for the effects of the other predictors in that model. The raw data relationships can be found in Fig. S2.
Snail abundance
The strongest predictors of total snail abundance were turbidity and effective fetch, with snails being more abundant at sites with lower turbidity (high water clarity) and high effective fetch (Table S4). However, the effect of turbidity became non-significant when longitude was added to the model (Table S5). We also found a significant positive effect of local conifer tree abundance on total snail density (Table S4). Total snails at all sites were dominated by Pleurocera spp. snails; when we analysed Pleurocera and Lymnaea spp. independently, it was clear that these 2 genera responded to very different variables, and the statistical effects of turbidity and effective fetch were driven by Pleurocera spp. (Table S4). We also detected a positive effect of mean day temperature on Pleurocera spp. abundance, but this was only evident in the full model (Table S4). In contrast, the best predictors of Lymnaea spp. snail abundance were maximum lake depth and the local abundance of deciduous trees (Table 2; Figure 3B & D). Deciduous tree cover appeared in all 10 top models in the all-subsets analysis (Table S7). Maximum lake depth appeared in 8 of the top 10 models; the remaining 2 models contained ‘lake surface area’ which was highly correlated with maximum lake depth (r = 0.88; Table S7). The Lymnaea and Pleurocera models were not affected by adding latitude or longitude, which were nonsignificant in both models.
Submerged vegetation
We examined the possible environmental predictors of submerged aquatic vegetation, which was the second-best predictor of cercaria abundance. The best single predictor of submerged vegetation was sediment phosphorus at the site (Table 2; Figure 3F). After accounting for sediment phosphorus, we also found significantly higher vegetation levels at sites with a building present in the riparian zone (Table 2). Per cent forest and per cent cropland in the watershed were also significant predictors at the end of stepwise model selection, but they became non-significant after adding ‘Lake’ as a random effect (Table 2). Sediment phosphorus appeared in all 10 top models in all subsets regression and ‘buildings’ appeared in the top 9 top models (Table S8); however, ‘buildings’ became nonsignificant when longitude was added to the model (Tables 2 & S5).
Temperature
Temperature correlated strongly with maximum lake depth, which was a strong predictor of Lymnaea spp. abundance (Table S2). We therefore examined possible environmental predictors of 3 measures of water temperature: mean (overall average) temperature, mean day temperature (average of daily maxima), and mean night temperature (average of daily minima). We limited these analyses to only include predictors related to lake hydrology (e.g. maximum lake depth) or potential sources of shade (e.g. abundance of trees or aquatic vegetation; Table S3). For all 3 measures, temperature was lower in larger lakes. The best predictors of all 3 temperature indices were either maximum lake depth or lake surface area (Table S4), which were highly correlated with each other. We also detected a significant negative effect of submerged vegetation on mean daytime temperature following stepwise model selection, which might indicate a shading effect; however, this effect became non-significant when a random effect of ‘Lake’ was added to the model (Table S4). None of the temperature models were affected by addition of latitude or longitude, neither of which was a significant predictor in any of these models.
Discussion
Cercaria abundance
Consistent with previous studies in northern Michigan, avian schistosome abundance throughout our survey was driven by the density of Lymnaea spp. snails, suggesting that these are the most important snail hosts driving variation in Trichobilharzia spp. cercaria abundance in these northern Michigan lakes (Blankespoor and Reimink, Reference Blankespoor and Reimink1991; Muzzall et al., Reference Muzzall, Burton, Snider and Coady2003; Rudko et al., Reference Rudko, Reimink, Froelich, Gordy, Blankespoor and Hanington2018). Physa and Planorbella spp. snails are also known hosts of Trichobilharzia spp. parasites, but their inclusion in our analysis led to no increase in predictive power. Physa and other genera of snails might also harbour Trichobilharzia spp. parasites that are implicated in SI in other parts of Michigan (McMullen and Brackett, Reference McMullen and Brackett1948; McPhail et al., Reference McPhail, Rudko, Turnbull, Gordy, Reimink, Clyde, Froelich, Brant and Hanington2021), though evidence suggests that Trichobilharizia spp. from Planorbella snails may contribute to SI risk throughout the state of Michigan (Soper et al., Reference Soper, Raffel, Sckrabulis, Froelich, McPhail, Ostrowski, Reimink, Romano, Rudko and Hanington2023). Additionally, the observed abundance of Physa and Planorbella snails at our local sampling sites might not have always reflected their off-site abundances in the lakes we sampled. Many of our sampling sites were (by necessity) publicly accessible sandy or rocky beach habitats, which tended to be dominated by Lymnaea and (non-host) Pleurocera snails. However, many were adjacent to more-vegetated sites where Physa and Planorbella might have been more abundant. Interestingly, we detected relatively high levels of cercariae at some sites where we observed few or no Lymnaea spp. snails (Fig. S2). The most likely source of cercariae at these sites is from snails outside of the study site, either in deeper water than we surveyed or along the shore in either direction (Laman et al., Reference Laman, Boss and Blankespoor1984; Muzzall et al., Reference Muzzall, Burton, Snider and Coady2003). Past studies have suggested that cercaria are transported from off-site via waves or water currents (Upatham, Reference Upatham1974; Leighton et al., Reference Leighton, Zervos and Webster2000; Muzzall et al., Reference Muzzall, Burton, Snider and Coady2003; Rudko et al., Reference Rudko, Reimink, Froelich, Gordy, Blankespoor and Hanington2018; Sckrabulis et al., Reference Sckrabulis, Flory and Raffel2020). We found no effect of bird abundance, at least not based on sightings by our community scientists during the study period. This result might suggest that local abundance of definitive host birds might not be a limiting factor for sustaining parasite populations in these lakes, consistent with a past study (Rudko et al., Reference Rudko, McPhail, Reimink, Froelich, Turnbull and Hanington2022). However, bird visitation likely varies seasonally, and bird visitation data from earlier in the season might have been more predictive of late-summer cercaria production by infected snails due to the time it takes for snails to develop patent infections.
We also found a negative effect of submergent vegetation on cercaria abundance in our final model (Table 2; Figure 3C). This result is opposite a pattern widely reported for human schistosomes, where aquatic vegetation is sometimes the strongest environmental predictor of human schistosomiasis risk due to providing habitat for Biomphalaria and Bulinus host snails (Klumpp and Chu, Reference Klumpp and Chu1980; Boelee and Laamrani, Reference Boelee and Laamrani2004; Rohr et al., Reference Rohr, Sack, Bakhoum, Barrett, Lopez-Carr, Chamberlin, Civitello, Diatta, Doruska, De Leo, Haggerty, Jones, Jouanard, Lund, Ly, Ndione, Remais, Riveau, Schacht, Seck, Senghor, Sokolow and Wolfe2023). The observed negative relationship in the current study was consistent across the range of vegetation levels we quantified, and it was robust to whether Lymnaea spp. abundance was included as a covariate (F 1,36 = 5.81, P(F) = 0.021). The consistency of this result suggests that it was not simply a spurious correlation. Furthermore, the Lymnaea host snails most associated with cercaria abundance in the current study have very different habitat preferences than host snails for human schistosomes, preferring silty, sandy, or rocky substrates and cool deep-water habitats (Laman et al., Reference Laman, Boss and Blankespoor1984). Thus, it is plausible that vegetation has a different effect in this system. We explored possible a posteriori hypotheses to explain a negative relationship between avian schistosome cercariae and submerged vegetation in this system. We found literature support for 3 alternative mechanisms (Warren and Peters, Reference Warren and Peters1968; Gibson and Warren, Reference Gibson and Warren1970; Christensen, Reference Christensen1979). First, Warren and Peters (Reference Warren and Peters1968) found that some species of plants can act as accidental hosts for schistosome cercariae, which might remove cercariae from the water column. However, only a small proportion of the species tested were penetrated by cercariae, and none of those species are native to Michigan lakes (Warren and Peters, Reference Warren and Peters1968). Second, Gibson and Warren (Reference Gibson and Warren1970) showed that bladderworts (genus Utricularia), a native carnivorous plant common in MI lakes, is an effective predator of schistosome cercariae and can rapidly reduce cercaria abundance in water. However, we did not collect data on the species identities of aquatic plants in our study, making it difficult to evaluate the plausibility of this hypothesis. Third, Christensen (Reference Christensen1979) found that multiple types of floating vegetation can effectively block cercaria movement (Christensen, Reference Christensen1979). Such an effect might plausibly prevent cercariae movement into a beach area from off site, especially for avian schistosome cercariae that tend to swim toward the water surface in search of floating (duck) hosts (Feiler and Haas, Reference Feiler and Haas1988). This might also explain why the aquatic vegetation pattern was especially strong at sites with few Lymnaea spp. host snails, meaning most of the cercariae detected at these sites probably came from outside the study area. This idea is also consistent with studies showing that onshore wind is an important predictor of daily variation in swimmer's itch risk, presumably due to off-site cercariae being transported towards shore in water currents generated by onshore wind (Rudko et al., Reference Rudko, Reimink, Froelich, Gordy, Blankespoor and Hanington2018; Sckrabulis et al., Reference Sckrabulis, Flory and Raffel2020).
We found that the best predictor of submerged vegetation was levels of phosphorus in the sediment (Fig. 3F), which is plausible as phosphorus is believed to be a limiting nutrient for aquatic plant growth in northern MI lakes (Bole and Allan, Reference Bole and Allan1978). The presence of an artificial structure (‘buildings’) was a positive predictor of vegetation in the final model (Table 2), but this pattern was spatially confounded with a significant longitude effect. We also discounted an apparent positive effect of 2,4-D herbicide during stepwise model selection, since it seemed unlikely that higher herbicide concentrations led to greater vegetation growth. The reverse causal direction seemed more plausible, i.e. that locals apply 2,4-D in response to increased vegetation.
To our knowledge, ours is the first field study to detect a negative effect of aquatic vegetation on SI-causing cercaria abundance. However, relatively few field studies have quantified cercaria abundance directly from water samples, so we cannot rule out the possibility that this effect might be common. We did not distinguish specific types of vegetation, so we cannot draw conclusions about the effects of specific aquatic plant species. Nevertheless, our results suggest that managing aquatic plant life near beaches might provide promising new ways to help control swimmer's itch in MI lakes. At minimum, we advise caution for lake managers involved in controlling aquatic weeds to discern target weeds from potentially beneficial native plants like bladderwort, which can be easily confused with invasive plants like Eurasian milfoil. Lake managers might also consider retaining existing beds of floating aquatic vegetation located adjacent to swim sites in temperate inland lakes, especially in locations where SI risk is known to be associated with onsite wind (Sckrabulis et al., Reference Sckrabulis, Flory and Raffel2020), since floating vegetation may help to reduce the influx of cercariae produced off site. This would be similar to other proposed management strategies aimed at blocking the influx of SI-causing cercariae, such as the floating barriers suggested by Muzzall et al. (Reference Muzzall, Burton, Snider and Coady2003) and R. L. Reimink (personal communication). However, given the positive association between aquatic vegetation and human schistosomiasis risk (Klumpp and Chu, Reference Klumpp and Chu1980; Boelee and Laamrani, Reference Boelee and Laamrani2004; Rohr et al., Reference Rohr, Sack, Bakhoum, Barrett, Lopez-Carr, Chamberlin, Civitello, Diatta, Doruska, De Leo, Haggerty, Jones, Jouanard, Lund, Ly, Ndione, Remais, Riveau, Schacht, Seck, Senghor, Sokolow and Wolfe2023), more work is needed to verify the extent to which submerged vegetation is associated with reduced SI risk in temperate lakes outside of the current study area.
We also found a significant direct positive effect of sediment phosphorus on cercaria abundance, though this was only a significant predictor after accounting for effects of host snail abundance and submerged vegetation (Fig. 3E & F). This pattern is consistent with prior studies that have found positive bottom-up effects of phosphorus on snail populations, mediated by increased availability of periphyton (Johnson et al., Reference Johnson, Chase, Dosch, Hartson, Gross, Larson, Sutherland and Carpenter2007; Rohr et al., Reference Rohr, Schotthoefer, Raffel, Carrick, Halstead, Hoverman, Johnson, Johnson, Lieske, Piwoni, Schoff and Beasley2008b). Since the observed effect of phosphorus was only evident after controlling for local Lymnaea abundance, it may be driven by increased rates of cercaria production by individual snails with better nutrition (Civitello et al., Reference Civitello, Fatima, Johnson, Nisbet and Rohr2018).
Snail abundance
We found significant effects of turbidity and effective fetch on the total snail abundance at each site, as well as an unexpected negative effect of local conifer abundance (Table 2). The observed positive relationship between water clarity (the inverse of turbidity) and snail abundance is consistent with our a priori hypothesis that increased water clarity would increase growth rates of the periphyton eaten by snails, thereby supporting greater snail population sizes. However, this turbidity effect on total snail abundance was spatially confounded and became nonsignificant when longitude was added to the model. Furthermore, most of the snail communities sampled in this survey were dominated by Pleurocera spp. snails, which are not known to harbour SI-causing avian schistosomes (Blankespoor and Reimink, Reference Blankespoor and Reimink1991). The best predictors of Pleurocera spp. abundance were water clarity and fetch, plus a positive effect of mean day temperature (Table S4), and the numerical dominance of this snail genus was the primary reason for significant effects of turbidity and fetch on total snail abundance in our survey. These patterns are consistent with Pleurocera spp. life history, which live at lower (i.e. warmer) latitudes, have thick shells to resist damage from wave action, and are commonly found in the shallow-water zone of larger rivers and lakes (Clarke, Reference Clarke1981; Tiemann et al., Reference Tiemann, Cummings and Mayer2011; Dillon et al., Reference Dillon, Jacquemin and Pyron2013).
Lymnaea spp. snail abundance was primarily associated with deeper lakes and sites with higher deciduous tree abundance (Fig. 3B & D). All Lymnaea spp. snails observed or collected in this survey appeared to fall within the Lymnaea catescopium species complex ( = Stagnicola emarginata; Correa et al., Reference Correa, Escobar, Durand, Renaud, David, Jarne, Pointier and Hurtrez-Bousses2010), which is thought to prefer hard substrate (rock/cobble) habitats in colder, deeper-water lakes (Clarke, Reference Clarke1981; Laman et al., Reference Laman, Boss and Blankespoor1984). Compared to other species, L. catescopium is found in deeper water, often at depths greater than 8 m (Clarke, Reference Clarke1981). Our quadrat surveys only extended to approx. 80 cm depth, so Lymnaea spp. abundance at our sites is likely underestimated, especially in locations where these snails might be restricted to deeper water. Additionally, L. catescopium is typically found at higher latitudes (>40°N) with Michigan being at the southern end of their range (Clarke, Reference Clarke1981), so it is thought to be more adapted to colder temperatures. In our survey, all 3 temperature measures (daytime, nighttime and daily mean temperature) were lower in lakes with greater surface area or maximum depth. It seems plausible that the observed positive relationship between Lymnaea abundance and lake depth was driven in part by the availability of deep-water habitats with cooler water temperatures, despite not detecting significant direct effects of local (site-level) water temperature measurements in our analysis.
To our knowledge, this is the first study to detect a relationship between deciduous shade trees and Lymnaea spp. abundance, and it would be necessary to conduct further study to determine whether this relationship might be causal. Nevertheless, this positive association might also be plausibly driven by providing definitive bird host habitat or a preference for cool temperatures by Lymnaea spp. snails (Fig. S2C & D). During the day when we conducted our surveys, Lymnaea spp. snails might only remain in the shallows if there are enough shade trees available to keep the water cool. We did not detect a clear effect of deciduous trees on local water temperature in our analysis, but our temperature measurements were taken at approximately 60 cm depth. Any effects of tree shade on water temperature are more likely to occur in even shallower near-shore water, where we sampled snails but did not collect temperature data.
Conclusions
It is important to emphasize that all patterns reported in this study are correlational, and follow-up experimental work would be needed to demonstrate causality. This is particularly true for novel findings, such as the apparently negative effect of vegetation on cercaria abundance, or the apparently positive effect of deciduous tree cover on local Lymnaea spp. abundance. It is also likely that the importance of some variables might depend on temporal or spatial scale. It remains possible that long-term changes in water clarity might have influenced Lymnaea populations in MI lakes over the last few decades, but if this effect was uniform across the study area, then our spatial survey would not have detected it.
Another common problem in ecological systems is that effects can be highly context dependent. For example, the most important drivers of snail abundance in northern MI lakes might differ from the most important drivers in southern MI lakes where the dominant parasite species appears to be a recently described avian schistosome that uses Planorbella snails as an intermediate host (McPhail et al., Reference McPhail, Rudko, Turnbull, Gordy, Reimink, Clyde, Froelich, Brant and Hanington2021; Soper et al., Reference Soper, Raffel, Sckrabulis, Froelich, McPhail, Ostrowski, Reimink, Romano, Rudko and Hanington2023). Prior studies found evidence of resource limitation in some northern Lymnaea populations (Cuker, Reference Cuker1983; O'Brien et al., Reference O'Brien, Barfield, Bettez, Hershey, Hobbie, Kipphut, Kling and Miller2005), suggesting that growth of attached algae (and by extension water clarity and/or nutrients) might sometimes be important drivers of these snail populations. Other studies have found an important role for top-down regulation of Lymnaea populations by fish and invertebrate predators (Hershey, Reference Hershey1990; Merrick et al., Reference Merrick, Hershey and McDonald1991), which we were unable to assess in this study. Future studies are needed to expand the geographical scope of our knowledge about SI drivers in temperate inland lakes and use manipulative experiments to test whether the relationships reported here are likely to be causal.
While this study didn't directly test methods for reducing avian schistosome abundance, identifying environmental predictors of risk does provide information that people could use to inform or develop risk management strategies. In particular, our results indicate that local vegetation (i.e. tree cover and submerged vegetation) might cercaria abundance at the site level. These are factors that could be easily modified by local landowners. If future studies confirm that these relationships are causal, this could open up new possibilities for managing swimmer's itch in inland lakes.
Supplementary material
The supplementary material for this article can be found at https://doi.org/10.1017/S0031182024000337.
Data availability statement
All code and data will be made available on GitHub upon acceptance: https://www.github.com/jasonsckrabulis/sckrabulis_etal_si_survey_2016.
Acknowledgements
We thank K.A. Altman and K.R. Sckrabulis for their assistance with field sampling. We thank M.D. Ostrowski for their assistance with zooplankton counting. We thank D.C. Szlag for their input on water chemistry methodology. We thank our community scientists for collecting daily filtered water samples and other data, without which this study would not be possible; please see a full list in the Supplementary Materials (Table S9). We thank C.L. Blankespoor and R.L. Reimink for allowing us to store field samples at their facility during the survey period. We thank an anonymous reviewer for helpful comments regarding the presentation of alternate analytical methods.
Authors’ contributions
J.P.S., M.L.M. and T.R.R. conceived and designed the study; J.P.S., M.L.M., J.S., R.B.M., A.M.H. and A.B. conducted the field survey; M.L.M. developed the qPCR protocol; M.L.M. and J.S. conducted qPCR; H.D.A. conducted all nutrient analyses; A.M.H. quantified all mussel and crayfish data. A.B. quantified agrochemical concentrations. J.P.S. and T.R.R analysed the data; J.P.S. wrote the manuscript; all authors reviewed and edited the manuscript.
Financial support
J.P.S. and M.L.M. were supported by Oakland University Teaching Assistantships. A.M.H. and A.B. were supported by Oakland University Honor's College research grants. Partial support was provided by NSF awards to T.R.R (CAREER: IOS-1651888 and IOS-1121529). Additional support for J.P.S. was provided by an NSF award to J.R. Rohr (DEB-2017785). Partial support was provided by a grant awarded by Michigan Swimmer's Itch Partnership and managed by the Tip of the Mitt Watershed Council to T.R.R. and D.M. Soper. External funding was provided by the following lake associations and organizations: Crystal Lake (Crystal Lake & Watershed Association), East & West Twin Lakes (Twin Lakes Property Owners Association), Elk & Skegemog Lakes (Elk-Skegemog Lakes Association), Glen Lake (Glen Lake Association Inc.), Hamlin Lake (Hamlin Lake Preservation Society), Higgins Lake (Higgins Lake Swimmer's Itch Organization), Intermediate Lake (Intermediate Lake Association), Lake Margrethe (Lake Margrethe Property Owners Association and Community Foundation of Northeast Michigan), Lime and Little Traverse Lakes (Lime Lake Association, Little Traverse Lake Association and the County of Leelanau), Platte Lake (Platte Lake Improvement Association), Portage Lake (Portage Lake Association and Alliance for Economic Success).
Competing interests
None.
Ethical standards
Not applicable.
Open access
