Introduction
Interactions between free-living species are the hallmarks of animal and plant communities (Gotelli and Graves, Reference Gotelli and Graves1996; Morin, Reference Morin2011). Positive interactions reflect symbiotic and mutualistic relationships, as for example among nectar-feeding insects pollinating flowers, while negative interactions in which the abundance of 1 species is reduced in the presence of another, may arise through predation as part of the food chain or through competition for limited resources (Morin, Reference Morin2011).
Similarly, in disease ecology interactions between different parasite species in hosts may be either positive or negative (Montgomery and Montgomery, Reference Montgomery and Montgomery1990; Haukisalmi and Henttonen, Reference Haukisalmi and Henttonen1998). The latter arise because for 1 or more reasons the host environment is made less suitable for the second species when the first is present. In consequence, the second species may be less likely to establish in the host, its growth and development may be delayed or impaired entirely, it may be excluded from the optimum site for its development in the host and/or its longevity may be curtailed in hosts concurrently infected with the first species. The underlying mechanisms, often based on cross immunity, generation of non-specific inflammatory responses, parasite-induced host mortality or competition for resources in concurrently infected hosts, have been reviewed comprehensively in earlier publications (Cox, Reference Cox2001; Behnke et al. Reference Behnke, Bajer, Siński and Wakelin2001a; Pedersen and Fenton, Reference Pederson and Fenton2007; Cattadori et al. Reference Cattadori, Boag and Hudson2008).
Positive interactions among infectious agents, resulting in co-occurrence in hosts, may arise because the first species modifies the host’s internal environment in a way that improves the establishment, development and survival of the second species. The most likely such mechanism is through depression of the host’s defences by the first species (Pritchard and Behnke, Reference Pritchard and Behnke1985; Hewitson et al., Reference Hewitson, Harcus, Murray, van Agtmaal, Filbey, Grainger, Bridgett, Blaxter, Ashton, Ashford, Curwen, Wilson, Dowle and Maizels2011; Maizels et al., Reference Maizels, Hewitson, Murray, Harcus, Dayer, Filbey, Grainger, McSorley, Reynolds and Smith2012), as part of its survival strategy, with the consequent knock-on effect of also benefiting the survival of other co-infecting species (Behnke et al., Reference Behnke, Wakelin and Wilson1978, Reference Behnke, Ali and Jenkins1984, Reference Behnke, Eira, Rogan, Gilbert, Torres, Miquel and Lewis2009).
Interactions among parasites are an important interest because they affect the way parasite communities assemble and in turn, the progression of disease within individual hosts and the onward transmission of infection between hosts (Graham et al., Reference Graham, Cattadori, Lloyd-Smith, Ferrari and Bjørnstad2007; Bradley and Jackson, Reference Bradley and Jackson2008). Through acting as regulatory feedbacks, interactions could even affect the long-term stability of host–parasite assemblages (Jackson et al., Reference Jackson, Bajer, Behnke-Borowczyk, Gilbert, Grzybek, Alsarraf and Behnke2023). Positive and negative interactions between parasites might sometimes be detectable as statistical associations in count or prevalence data from surveys of free-running systems. Unfortunately, however, a number of coincidental processes may produce such statistical associations, confounding the ability to detect true interactions that could act as genuine regulatory feedbacks (Bradley and Jackson, Reference Bradley and Jackson2008; Barger, Reference Barger2023). Such coincidental processes result from correlated susceptibility or correlated exposure (Bradley and Jackson, Reference Bradley and Jackson2008). Correlated exposure may arise because the infective stages of multiple species co-locate in proximity in the same external environment as that resided in by the host and concurrent transmission is therefore likely. For example, the definitive host of multiple parasite species may feed on an animal utilized by all of these parasites as an intermediate host for their larval stages or on one that acts as a paratenic host for the parasites concerned (Behnke et al., Reference Behnke, Jackson, Gilbert, Mohallal and Bajer2024). Correlated susceptibility may occur because of genetic variation in the host that simultaneously affects resistance to multiple species, or because of environmental factors that affect host quality and the general ability of hosts to defend themselves against parasitic infection.
Parasite–parasite interactions have been demonstrated clearly many times in laboratory systems (Holmes, Reference Holmes1961; Cox, Reference Cox2001; Behnke et al., Reference Behnke, Bajer, Siński and Wakelin2001a) but the evidence for their importance in more complex real-world situations, where infection exposures and host phenotypes may be very different to those in the laboratory, remains somewhat contradictory (Poulin, Reference Poulin2001; Telfer et al., Reference Telfer, Birtles, Bennett, Lambin, Paterson and Begon2008, Reference Telfer, Lambin, Birtles, Beldomenico, Burthe, Paterson and Begon2010: Knowles et al., Reference Knowles, Fenton, Petchey, Jones, Barber and Pedersen2013; Keegan et al., Reference Keegan, Pedersen and Fenton2024) and an active area of research. While field-scale experiments (Knowles et al., Reference Knowles, Fenton, Petchey, Jones, Barber and Pedersen2013; Fenton et al., Reference Fenton, Knowles, Petchey and Pedersen2014) or analysis of individual time series (longitudinal study designs) (Telfer et al., Reference Telfer, Birtles, Bennett, Lambin, Paterson and Begon2008) provide a potentially powerful mode of inference in natural systems, that could partition genuine interactions from other processes and reveal their true sign, nonetheless they are also subject to very significant practical and technical difficulties. Thus, large-scale observational surveys (cross-sectional study designs), that are easier to execute and replicate, and that search for evidence of an interaction through statistical association in species counts or prevalence, are still relevant. In particular, consistent patterns of association, linked to individual taxa and predictable from the known biological propensities of those taxa, might still be taken as important corroborating evidence supporting the existence of interactions in the real world. Here, we present such an example, based on a very large survey of a community that included species from a group with known immunosuppressive properties.
Nematode parasites of the family Heligmosomidae are parasites of rodents, and 1 species in particular, Heligmosomoides bakeri has been shown through extensive laboratory based experimental infections to give rise to chronic infections (Robinson et al., Reference Robinson, Wahid, Behnke and Gilbert1989) and to be immunomodulatory in its murine host as part of its survival strategy in the face of the adaptive immune system (Behnke et al., Reference Behnke, Hannah and Pritchard1983; Hewitson et al., Reference Hewitson, Harcus, Murray, van Agtmaal, Filbey, Grainger, Bridgett, Blaxter, Ashton, Ashford, Curwen, Wilson, Dowle and Maizels2011; Maizels et al., Reference Maizels, Hewitson, Murray, Harcus, Dayer, Filbey, Grainger, McSorley, Reynolds and Smith2012; Smith et al., Reference Smith, Filbey, Reynolds, Hewitson, Harcus, Boon, Sparwasser, Hammerling and Maizels2016; Filbey et al., Reference Filbey, Mehta, Meijlink, Pellefigues, Schmidt and Le Gros2020). A possible consequence of such a strategy is enablement of longer survival of other concurrently resident species (Jenkins and Behnke, Reference Jenkins and Behnke1977; Behnke et al., Reference Behnke, Wakelin and Wilson1978, Reference Behnke, Ali and Jenkins1984). The mechanism of immunodepression is complex and has been reviewed (Maizels et al., Reference Maizels, Hewitson, Murray, Harcus, Dayer, Filbey, Grainger, McSorley, Reynolds and Smith2012, Reference Maizels, Smits and McSorley2018; Shepherd et al., Reference Shepherd, Navarro, Wangchuk, Wilson, Daly and Loukas2015; Yap and Gause, Reference Yap and Gause2018; McManus and Maizels, Reference McManus and Maizels2023). Significantly, H. bakeri’s closely related congeneric species H. polygyrus has been shown to accumulate with age in naturally occurring Apodemus sylvaticus hosts (Gregory, Reference Gregory1992; Gregory et al., Reference Gregory, Montgomery and Montgomery1992; Behnke et al., Reference Behnke, Lewis, Mohd Zain and Gilbert1999; Friberg et al., Reference Friberg, Little, Ralli, Lowe, Hall, Jackson and Bradley2013; Loxton et al., Reference Loxton, Lawton, Stafford and Holland2017). The occurrence of H. polygyrus is furthermore associated with increasing species richness of other helminth species under natural conditions (Behnke et al., Reference Behnke, Gilbert, Abu-Madi and Lewis2005, Reference Behnke, Eira, Rogan, Gilbert, Torres, Miquel and Lewis2009), although antagonistic (negative) interactions with coccidia (Protozoa) have also been documented in field experiments (Knowles et al., Reference Knowles, Fenton, Petchey, Jones, Barber and Pedersen2013). Apart from these studies on H. polygyrus, little is known about whether other naturally occurring heligmosomid nematodes employ immunodepression for survival in definitive hosts.
Two heligmosomid nematodes, Heligmosomoides glareoli and Heligmosomum mixtum are dominant helminths of bank voles (C. glareolus; Haukisalmi and Henttonen, Reference Haukisalmi and Henttonen1993a, Reference Haukisalmi and Henttonen1993b) in continental Europe. The prevalence of both species in bank voles may be high, while abundance seldom exceeds 10 worms/host (Haukisalmi, Reference Haukisalmi1986; Haukisalmi and Henttonen, Reference Haukisalmi and Henttonen1993a; Grzybek et al., Reference Grzybek, Bajer, Bednarska, Alsarraf, Behnke-Borowczyk, Harris, Price, Brown, Osborne, Siński and Behnke2015). Like H. polygyrus, both accumulate in their hosts with increasing host age (Grzybek et al., Reference Grzybek, Bajer, Bednarska, Alsarraf, Behnke-Borowczyk, Harris, Price, Brown, Osborne, Siński and Behnke2015), suggesting that they give rise to chronic infections. If these species, like their congeners modify the host environment to benefit their own long-term survival, then we may expect the prevalence and abundance of other helminth species, as well as helminth species richness to be enhanced in heligmosomid-infected bank voles. To test this prediction we analysed the database, utilized by Grzybek et al. (Reference Grzybek, Bajer, Bednarska, Alsarraf, Behnke-Borowczyk, Harris, Price, Brown, Osborne, Siński and Behnke2015), which comprises over 900 bank vole records of animals sampled at 3- to 4-year intervals over more than a decade in 3 disparate sites in NE Poland, specifically seeking evidence for a positive association between the heligmosomids and other helminth taxa.
Materials and methods
Study sites and source of data
The following analysis utilized the database of Grzybek et al. (Reference Grzybek, Bajer, Bednarska, Alsarraf, Behnke-Borowczyk, Harris, Price, Brown, Osborne, Siński and Behnke2015). These authors reported on the long-term spatiotemporal stability and changes in the helminth communities of bank voles from 3 study sites (Urwitałt, Tałty and Pilchy) in the temperate forests of NE Poland, named after nearby villages. The animals were trapped during 4 surveys conducted initially at 3, then 4-year intervals (1999, 2002, 2006 and 2010) in late summer (August–September) in each of these years. Descriptions of the sites, the methodology used in fieldwork and for subsequent laboratory analysis of worm burdens have all been described comprehensively in earlier publications (Behnke et al., Reference Behnke, Barnard, Bajer, Bray, Dinmore, Frake, Osmond, Race and Siński2001b, Reference Behnke, Bajer, Harris, Newington, Pidgeon, Rowlands, Sheriff, Kuliś-Malkowska, Siński, Gilbert and Barnard2008a, Reference Behnke, Bajer, Harris, Newington, Pidgeon, Rowlands, Sheriff, Kuliś-Malkowska, Siński, Gilbert and Barnard2008b; Grzybek et al., Reference Grzybek, Bajer, Bednarska, Alsarraf, Behnke-Borowczyk, Harris, Price, Brown, Osborne, Siński and Behnke2015). Based on morphometric data collected at autopsy and dried eye lens weight, the animals were allocated to 3 age classes (age classes 1, 2 and 3, corresponding to immature, young, adult and mature older voles, respectively; see Behnke et al., Reference Behnke, Barnard, Bajer, Bray, Dinmore, Frake, Osmond, Race and Siński2001b) and were sexed (male or female).
Definition of terms employed
In this paper, we refer to Clethrionomys glareolus (as recommended by the Mammal Diversity Database; Upham et al., Reference Upham, Burgin, Widness, Liphardt, Parker, Becker, Rochon and Huckaby2023; Tesakov et al., Reference Tesakov, Lebedev, Bannikova and Abramson2010; Kryštufek et al., Reference Kryštufek, Tesakov, Lebedev, Bannikova, Abramson and Shenbrot2020), rather than Myodes glareolus (see Carleton et al., Reference Carleton, Gardner, Pavlinov and Musser2014), following several recent revisions of the generic attribution of bank voles.
We focus on the effect of 2 species of heligmosomid nematodes (Heligmosomidae; Heligmosomoides glareoli and Heligmosomum mixtum) on other helminth taxa. These 2 species were treated as 1 taxon and are referred to hereafter as the heligmosomids. We refer also to H. bakeri, a parasite of Mus spp., previously known as H. polygyrus and Nematospiroides dubius (Behnke et al., Reference Behnke, Keymer and Lewis1991; Cable et al., Reference Cable, Harris, Lewis and Behnke2006). All other species of helminths were pooled into the following groups for analysis:
1. Non-heligmosomid helminths (all other helminth species, referred to also as other helminths),
2. Non-heligmosomid nematodes (all other nematode species, referred to also as other nematodes),
3. Oxyuroidea (Aspiculuris tianjinensis and Syphacia spp.),
4. Cestodes (all cestode species and individual worms, both immature and mature intestinal strobilate stages and tissue dwelling larvae),
5. Intestinal cestodes (adult, strobilate stages from intestines) and
6. Larval cestodes (all larval stages, found in the pleural and peritoneal cavities, and associated organs).
For further details of all the species involved, see Table 1.
Table 1. The overall prevalence and abundance of helminths recovered from Clethrionomys glareolus in 3 study sites in NE Poland during 4 surveys at 3- to 4-year intervals (n = 922)
Most of the data in this table were included in the text in Grzybek et al. (Reference Grzybek, Bajer, Bednarska, Alsarraf, Behnke-Borowczyk, Harris, Price, Brown, Osborne, Siński and Behnke2015).
Statistical analysis
For data subsets, prevalence values (% of voles infected with the helminth taxon referred to) are given with 95% exact confidence limits (CL95) in the tables and text, and 95% confidence intervals (CI95) are illustrated in the figures. These were calculated in bespoke software based on the tables of Rohlf and Sokal (Reference Rohlf and Sokal1995). For quantitative data (i.e. abundance as reflected in worm burdens), mean values ± standard error of the mean (SEM) are provided and these refer to abundance as defined by Bush et al. (Reference Bush, Lafferty, Lotz and Shostak1997) and Margolis et al. (Reference Margolis, Esch, Holmes, Kuris and Schad1982) and are derived for all voles in a data subset, including those uninfected by the worms of the relevant taxon.
The frequency distributions of worm burdens of all the taxa employed for analysis were tested by x 2 for goodness of fit to positive binomial, negative binomial (negbin), Poisson and Gaussian distributions as relevant. We provide values for the index of dispersion (mean/variance ratio, I), index of discrepancy (D; Poulin, Reference Poulin1993), and the negative binomial exponent (k) for each of the heligmosomid taxa separately and when combined (Figure S1).
For consistency with our previous work (Grzybek et al., Reference Grzybek, Bajer, Bednarska, Alsarraf, Behnke-Borowczyk, Harris, Price, Brown, Osborne, Siński and Behnke2015), prevalences were analysed using maximum likelihood techniques based on log linear analysis of contingency tables in the software package IBM SPSS, version 28 (Armonk, New York, USA). Full factorial models included SEX (2 levels, male and female voles), AGE (3 levels, juvenile, young adult and fully mature voles), geographical locality, i.e. SITE (3 levels corresponding to the 3 study sites as specified above), YEAR of survey (1999, 2002, 2006 and 2010), and these are referred to as the intrinsic (sex and age class) and extrinsic (site and year) factors. Models also included presence/absence (P/A) of infection with the heligmosomid worms. This approach is based on categorical values of factors of interest, which are used not only to fit hierarchical log-linear models to multidimensional cross-tabulations using an iterative proportional-fitting algorithm but also to detect associations between the fitted factors. Multi-factorial models were fitted as described previously (Grzybek et al., Reference Grzybek, Bajer, Bednarska, Alsarraf, Behnke-Borowczyk, Harris, Price, Brown, Osborne, Siński and Behnke2015), beginning with the full factorial model, followed by simplification of the model until only significant terms remained. Thus, for each level of analysis in turn, beginning with the most complex model, involving all possible main effects and interactions, those combinations that did not contribute significantly to explaining variation in the data were eliminated in a stepwise fashion, beginning with the highest-level interaction (backward selection procedure). A minimum sufficient model (MSM) was then obtained, for which the likelihood ratio of Chi-squared was not significant, indicating that the model was sufficient in explaining the data. The importance of each term in interactions involving P/A of a specified taxon in the final model was assessed by the probability that its exclusion would alter the model significantly and these values are given in the text, assessed by a likelihood ratio test between models with and without each term of interest.
For analysis of quantitative data, we first fitted full factorial (FFM) generalized linear (GLM) statistical models in R (v. 4.3.2, R Core Development Team) with both extrinsic (YEAR and SITE) and intrinsic (SEX and AGE) fixed factors for all dependent variables, testing models that were based on negative binomial or Poisson error structures, as relevant. The dependent variables were worm burdens of the groups listed in Table 2 and helminth species richness of all taxa other than the heligmosomids. Wherever possible, we used the ‘step’ command in R to refine these to minimum sufficient models (MSMs). However, in some cases, full factorial GLMs could not be simplified, and in others only reduced models (best-fit model, BFMs) would converge appropriately, and in this latter case, we specify in the text or in parenthesis the factors that were successfully fitted. Next, we combined the FFMs, MSMs, or BFMs with the addition of P/A of heligmosomid worms as a factor at 2 levels (present or absent). To evaluate whether worm burdens of other taxa differed significantly between voles with/without heligmosomids we compared models with and without (P/A) heligmosomids.
Table 2. Prevalence and abundance of non-heligmosomid helminth (i.e. Not Heligmosomidae) species and higher taxa in bank voles with and without heligmosomids. Data have been pooled, and include all years and sites, and both host sex and age classes
The higher prevalences and abundances at each level are highlighted in bold.
To test whether increasing worm burdens of heligmosomids were associated with increasing worm burdens of specified taxa (i.e. whether there was quantitative covariance between the worm burdens of these taxa), we combined models of those taxa, as fitted above, with the (deviance) residuals of the MSM for heligmosomid worm burdens but fitted as a covariate. This approach removes correlation between heligmosomid worm burdens and the other explanatory variables and is conservative (i.e. only assigning variance explained by the heligmosomid count variable over and above that already explained by other correlated explanatory variables). Furthermore, effects from these models could be visualized on the scale of the response worm burdens. We also examined in R, linear covariance of the residuals of FFMs, MSMs or BFMs, as relevant, of specified taxa with the residuals of the MSM for heligmosomid worms and then tested whether slopes (ß ± SEM) were positive and whether they differed significantly from zero. This approach not only deals with correlation between the explanatory variable of interest and other explanatory variables (as above) conservatively, but also preserves the maximum amount of degrees of freedom for the key test of the hypothesis. In all these cases, we predicted positive slopes, and hence we employed 1-tailed values for P. Figures of the resulting relationships were constructed in the ggplot2 and sjPlot packages of R. Values of likelihood ratios (LR) for models based on negative binomial errors and deviances (Dev) for models with Poisson errors are provided.
For analysis of ordinal data (e.g. age class) we used the Mann–Whitney U test in the software package IBM SPSS, version 28.
The covariance distribution test of Haukisalmi and Henttonen (Reference Haukisalmi and Henttonen1998), which focuses on the frequency distribution of pairwise covariances (i.e. positive and negative covariances [cov] between all pairs of species in a community), was applied to the data. We used both log10 + 1 and square root transformations of worm counts, and then obtained significance values from assessment of the number of times the observed covariance for each pair of species was less than or greater than that derived from each of 5000 cycles of randomization of the data. The test was then repeated with subsample constraints, i.e. with year, site, and host age and sex, taken into consideration.
Results
Age of voles that were infected or not infected with heligmosomid nematodes
Helminths are known to accumulate in hosts exposed to natural infections, such that worm burdens and helminth species diversity both increase with host age (Montgomery and Montgomery, Reference Montgomery and Montgomery1989; Behnke et al., Reference Behnke, Lewis, Mohd Zain and Gilbert1999; Grzybek et al., Reference Grzybek, Bajer, Bednarska, Alsarraf, Behnke-Borowczyk, Harris, Price, Brown, Osborne, Siński and Behnke2015). Therefore, it was necessary to determine first whether the 2 key groups in this study (voles infected with heligmosomid nematodes and those without heligmosomids) differed in age. Since the voles had been allocated to 3 age classes (1, 2 or 3) based on a range of morphological features, we averaged the age classes. Voles without heligmosomids had an average age class of 1.81 ± 0.039, while those infected with heligmosomids 2.22 ± 0.034. This difference is significant (Mann–Whitney U test, z = 7.754, P < 0.001).
Overall prevalence and abundance of helminths in bank voles
Table 1 gives the overall prevalence and mean worm burdens of the helminths recovered from C. glareolus, with data combined across the 3 study sites, 4 surveys, both sexes and the 3 age classes. Collectively, nematodes accounted for most of the helminths recovered (80.5% of all the helminths recovered were nematodes). For a breakdown of these data by the 4 factors and relevant analyses, see Grzybek et al. (Reference Grzybek, Bajer, Bednarska, Alsarraf, Behnke-Borowczyk, Harris, Price, Brown, Osborne, Siński and Behnke2015).
Across the whole study period the most prevalent species was A. tianjinensis (42.1%), and individually the 2 heligmosomid species were the next most prevalent, but when combined their prevalence was 55.5%. The frequency distributions of worm burdens with H. mixtum, H. glareoli and the 2 heligmosomids combined are shown in Fig. S1. For each taxon the best-fit distribution was the negative binomial: for H. mixtum (χ 211 = 14.1, P = 0.229); H. glareoli (χ 28 = 16.9, P = 0.031); heligmosomids combined (χ 213 = 16.3, P = 0.234).
Does the presence of heligmosomids affect prevalence of other helminths?
The prevalence of non-heligmosomid helminths was higher in bank voles with (69.7% [66.76–72.55%]), compared to those without heligmosomid nematodes (54.4% [48.17–60.55%]). In a model without extrinsic (YEAR and SITE) or intrinsic (AGE and SEX) factors, this difference was highly significant (χ 21 = 22.9, P < 0.001; Table 2). However, when all 4 factors were fitted, this relationship was dependent on both site and age of the voles (P/A heligmosomids × P/A other helminths × SITE × AGE CLASS, χ 24 = 11.3, P = 0.024), as illustrated in Figure 1. Thus, while prevalence of non-heligmosomid helminths was numerically higher overall in bank voles that carried heligmosomid nematodes in each major subset of data (i.e. by year, site, sex, and age; Table 3), there was a significant difference in how prevalence varied with host age between the 3 study sites. Although, it can be seen in Figure 1 that prevalence of non-heligmosomid helminths was mostly higher in voles also carrying heligmosomids in each of the 3 age classes, especially at Pilchy, the bias was in the opposite direction among age class 1 voles from Urwitałt and age class 2 voles from Tałty.
Figure 1. Prevalence of non-heligmosomid helminths in bank voles infected or not infected with heligmosomid nematodes, by year and site.
Table 3. Prevalence of non-heligmosomid helminths (i.e. not Heligmosomidae) and non-heligmosomid nematodes in bank voles with and without heligmosomids by year, site, host sex and age class
The higher prevalence at each level is highlighted in bold.
Does the presence of heligmosomids affect prevalence of other nematodes?
A similar picture emerged when the analysis was restricted to non-heligmosomid nematodes. Prevalence of non-heligmosomid nematodes was higher overall in bank voles with (59.6% [56.45–62.62%]) compared with those without heligmosomid nematodes (48.3% [42.07–54.51%]), and in each of the data subsets (Tables 2 and 3). In a model without extrinsic (YEAR and SITE) or intrinsic (AGE and SEX) factors, this difference was highly significant (χ 21 = 11.7, P < 0.001) and retained independent significance when all 4 factors and their interactions were taken into consideration (P/A heligmosomids × P/A other nematodes, χ 21 = 4.69, P = 0.030).
Restricting the analysis further to the Oxyuroidea (Table 2), yielded a weaker relationship (P/A heligmosomids × P/A Oxyuroidea, χ 21 = 4.30, P = 0.038) when the 4 factors were not included in the model, but lost significance when the 4 factors were taken into account (χ 21 = 3.45, P = 0.063). We found a similar outcome for A. annulosa, with a significant difference when no factors were included (χ 21 = 6.51, P = 0.011; Table 2) but loss of significance when we controlled for both intrinsic and extrinsic factors (χ 21 = 0.52, P = 0.47). Likewise, the overall prevalence of M. muris was also numerically higher in voles with heligmosomid nematodes compared to those without and this difference was highly significant (χ 21 = 14.4, P < 0.001; Table 2). However, in this case inclusion of the 4 factors in the model yielded 2 significant but complex terms (P/A heligmosomids × P/A M. muris × YEAR × SITE × SEX, χ 26 = 15.5, P = 0.017; P/A heligmosomids × P/A M. muris × SITE × SEX × AGE, χ 24 = 20.2, P < 0.001), which were not explored further.
Does the presence of heligmosomids affect prevalence of cestodes?
In a model without the 4 factors, the prevalence of cestodes (all stages combined) was significantly higher among voles infected with heligmosomid nematodes (Table 2; χ 21 = 31.3, P < 0.001). This retained significance when the 4 factors were taken into account but was confounded by interactions with YEAR and SITE (P/A heligmosomids x P/A cestodes × SITE × YEAR, χ 26 = 17.8, P = 0.007). This interaction arose because while the difference between prevalence in heligmosomid infected and not-infected voles was always biased in favour of the former (in all 12 data subsets illustrated in Figure 2), the extent of the difference varied between years and sites.
Figure 2. Prevalence of cestodes (all species and stages combined) in bank voles infected or not infected with heligmosomid nematodes, by age class and site.
Prevalence of intestinal, strobilate stages of cestodes was significantly higher in voles infected with heligmosomid nematodes (Table 2), and independent of all 4 factors and interactions between them (P/A heligmosomids x P/A intestinal cestodes, χ 21 = 4.22, P = 0.04). The bias in favour of higher prevalence of intestinal cestodes among heligmosomid-infected voles was evident in each year and site, and in both sexes and all 3 age classes (Table 4). A similar model with larval cestode stages revealed that there was also a bias for higher prevalence among heligmosomid-infected animals, but that it was confounded by host sex (P/A heligmosomids × P/A larval cestodes × SEX, χ 21 = 5.52, P = 0.019). While prevalence of larval cestodes was higher in heligmosomid-infected voles in each year, site, sex and age class (Table 4), in the case of sex, prevalence was 1.9 times higher in heligmosomid-infected male voles, but the disparity was much greater in female voles (7.2 higher in heligmosomid-infected voles; Table 4).
Table 4. Prevalence of intestinal strobilate and larval cestodes in bank voles with and without heligmosomids by year, site, host sex and age class
The higher prevalence at each level is highlighted in bold. Sample sizes are as in Table 3.
Does the presence of heligmosomid nematodes affect species richness with non-heligmosomid helminths?
The overall mean species richness of non-heligmosomid helminths was 0.91 ± 0.029. However, in the presence of heligmosomid nematodes this was higher (heligmosomids present, 1.07 ± 0.042 and absent, 0.70 ± 0.038; GLM with Poisson errors, Dev 1,920 = 35.3, P < 0.0001). As the data in Table 5 show, the mean species richness of non-heligmosomid helminths was numerically higher in heligmosomid carrying voles in each of the 4 years of survey, in each of the 3 study sites, in both sexes and all 3 age classes. With both intrinsic and extrinsic factors taken into account this difference was highly significant (GLM with Poisson errors, MSM for non-heligmosomid helminth species richness + P/A heligmosomids, Dev 1,906 = 6.680, P = 0.0097), but explained only 0.89% of the variance.
Table 5. Mean species richness of non-heligmosomid helminths (i.e. not heligmosomidae) in bank voles with and without heligmosomids by year, site, host sex and age class
The higher prevalence at each level is highlighted in bold.
Sample sizes are as in Table 3.
Does the presence of heligmosomid nematodes affect abundance of other helminths?
Mean abundances of all the different taxa, in voles that were infected or not infected with heligmosomid nematodes, are shown in Table 2. In each case, mean worm burdens were numerically higher in voles that were infected with heligmosomids. Analysis of all the non-heligmosomid helminths, with all significant intrinsic and extrinsic factors and their interactions taken into account, showed that this difference between voles with and without heligmosomids was significant (GLM with negbin errors, FFM for non-heligmosomid helminths + P/A of heligmosomids, LR 1,850 = 6.228, P = 0.0126), explaining 0.12% of the variance. As the summary data in Table 6 show, mean worm burdens were numerically higher in voles with heligmosomids in each of the 4 surveys, in each of the 3 study sites, in both sexes and in each of the 3 age classes.
Table 6. Abundance of non-heligmosomid helminths (i.e. not Heligmosomidae) and cestodes (all stages combined) in bank voles with and without heligmosomids by year, site, host sex and age class
The higher abundance at each level is highlighted in bold. Sample sizes are as in Table 3.
When the analysis was confined to non-heligmosomid nematodes, the difference in abundance between voles with or without heligmosomids was also highly significant (GLM with negbin errors, P/A of heligmosomids, LR 1,920 = 29.68, P < 0.0001), but lost significance when all the intrinsic and extrinsic factors and their interactions were taken into account (GLM with negbin errors, FFM for non-heligmosomid nematodes + P/A of heligmosomids, LR 1,850 = 2.191, P = 0.139; Table S1).
Additional refinement of the analysis, confined to oxyuroid nematodes, showed that there was a highly significant effect of P/A of heligmosomid nematodes on the abundance of oxyuroid nematodes (GLM with negbin errors, P/A of heligmosomids, LR 1,920 = 22.26, P < 0.0001), accounting for 0.51% of variance. However, the full factorial model did not converge. A model with all 4 factors as main effects and with 2 interactions (SITE × YEAR and SEX × AGE) plus P/A heligmosomids converged but was marginally non-significant (GLM with negbin errors, BFM for oxyuroid nematodes + P/A of heligmosomids, LR 1,904 = 3.505, P = 0.061). A model without SEX, but with all 3 remaining main effects (YEAR, SITE and AGE) and three 2-way interactions (GLM with negbin errors, BFM for oxyuroid nematodes + P/A of heligmosomids, LR 1,897 = 3.928, P = 0.047) converged successfully but accounted for only 0.1% of variance. These results suggest a weak effect, after controlling for the intrinsic and extrinsic factors and the interactions that were taken into account. The summary data in Table S1 show that with 1 exception (in the year 2002), mean oxyuroid worm burdens were numerically higher in heligmosomid-infected hosts in other years (1999, 2006 and 2010), in all 3 sites, both sexes and all age classes.
We selected 2 of the nematode species that showed a prevalence exceeding 10% for analysis. The most prevalent oxyuroid nematode was A. tianjinensis (Table 1), and as above, without controlling for confounding factors, abundance was higher in voles with heligmosomid nematodes (Table 2; GLM with negbin errors, P/A of heligmosomids on abundance of A. tianjinensis, LR 1,920 = 4.33, P = 0.0374). Again, FFM for A. tianjinensis would not converge, however a model with YEAR + SITE + AGE + all three 2-way interactions did (GLM with negbin errors, BFM for A. tianjinensis + P/A of heligmosomids, LR 1,897 = 10.377, P = 0.0013, R 2 = 0.28%; see Table S2 for breakdown by main effects and levels within).
The abundance of M. muris was also numerically higher in voles concurrently harbouring heligmosomids (Table 2), but the difference was not significant when we fitted a GLM with negbin errors without controlling for intrinsic and extrinsic factors (LR 1,920 = 1.577, P = 0.209). Full factorial models for M. muris also did not converge but a model comprising all 4 main effects and 5 of the six 2-way interactions (omitting YEAR × AGE) did, revealing that the P/A of heligmosomids was not a significant effect (Table S2; GLM with negbin errors, P/A of heligmosomids, LR 1,895 = 2.212, P = 0.137, R 2 = 0.188%). Although, in 8 of the data subsets (Table S2) abundance was numerically higher in heligmosomid carrying voles, in 4, abundance was higher in those without concurrent heligmosomids.
The mean abundance of cestodes (all species combined) was numerically 10 times higher in voles with heligmosomids (Table 2), and without controlling for the intrinsic and extrinsic factors this was a highly significant difference (GLM with negbin errors, P/A of heligmosomids on abundance of cestodes, LR 1,920 = 46.39, P < 0.00001). However, models controlling for all possible confounding factors proved highly problematic. The most complex acceptable model was a GLM with negbin errors, SITE × YEAR × AGE and P/A of heligmosomids on abundance of cestodes (LR 1,885 = 14.39, P = 0.00015, R 2 = 0.773%). As the data in Table 6 show, abundance was numerically higher in voles also carrying heligmosomids in all 4 surveys, all 3 sites, in both sexes and in all 3 classes. The abundance of intestinal stages of cestodes and that of larval (intermediate stages) cestodes was also higher overall in heligmosomid carrying voles (Table 2; twice and 12 times higher, respectively), and in both cases significantly higher without controlling for confounding factors (GLMs with negbin errors and P/A of heligmosomids on abundance, LR 1,920 = 9.95, P = 0.0016 and LR 1,920 = 14.108, P = 0.00017, respectively). In voles concurrently infected with heligmosomids, mean abundances were numerically higher in all data subsets corresponding to the main effects (Table 7), however, full factorial models were problematic in each case. For intestinal cestodes, a BFM (comprising all four main effects + five 2-way interactions [excluding YEAR × AGE] and two 3-way interactions (SITE × YEAR × SEX and SITE × SEX × AGE) + P/A of heligmosomids on abundance, LR 1,885 = 0.885, P = 0.35, R 2 = 0.084%), indicated that the difference in abundance of intestinal cestodes in voles with and without heligmosomids was not significant. For larval cestodes only a simple model (SITE × AGE + P/A heligmosomids, LR 1,912 = 6.118, P = 0.0134, R 2 = 0.634%) converged and indicated a significant effect of the P/A of heligmosomid nematodes on abundance of larval cestodes, but must be interpreted with caution because neither YEAR, nor SEX were taken into account.
Table 7. Abundance of intestinal and larval cestodes in bank voles with and without heligmosomids by year, site, host sex and age class
The higher abundance at each level is highlighted in bold.
Sample sizes are as in Table 3.
Does the increasing intensity of heligmosomid burdens covary with increasing species richness of non-heligmosomid helminths?
To control for the possible confounding effects of intrinsic and extrinsic factors known to affect heligmosomid abundance and non-heligmosomid species richness, covariance between these was tested in the MSM for non-heligmosomid helminth species richness either combined with or without the residuals of the MSM for heligmosomids as a covariate. The residuals of the MSM for the heligmosomid burden were a highly significant explanatory term in this model (Table 8, R 2 = 1.32%) with a positive gradient (Figure 3A).
Figure 3. Covariance of the residuals of the Minimum Sufficient Model (MSM) for heligmosomid worm burdens with the residuals of the MSM for non-heligmosomid helminth species richness. Panel A illustrates the predictions of the model in Table 8, whereas panel B shows the regression (R 2 = 0.0148). The shaded areas show the 95% confidence region. For additional details, see text and Table 9.
Table 8. Quantitative covariance of non-heligmosomid species richness and worm burdens of non-heligmosomid helminth taxa with worm burdens of the heligmosomid nematodes. Models comprised the residuals of the minimum sufficient model (MSM) for heligmosomid burdens as a covariate fitted as an additional explanatory factor to full factorial (FFM), minimum sufficient (MSM) or best-fit (BFM) models of the specified taxa
a In each case, we first either fitted full factorial, minimum sufficient or best-fit models, derived for the specified taxon from full factorial models (YEAR × SITE × AGE × SEX), and then added the residuals of the minimum sufficient model for heligmosomid worm burdens (WB) as a covariate.
b Test statistic (LR, likelihood ratio, distributed as chi-squared or Dev, deviance) with degrees of freedom, its value and P for comparison of models with and without heligmosomid WB (covariate). P-values are 1-tailed since we had predicted positive covariance in all cases.
c Species richness of non-heligmosomid helminths. For this dependent variable, the MSM was based on Poisson errors. For all the other dependent variables, models were based on negative binomial error structures.
Additionally, we also tested the linear relationships between the residuals of the MSMs for non-heligmosomid helminth species richness and those of the MSM for heligmosomid burdens. This relationship was positive (Table 9 and Figure 3B, r = 0.126) and highly significant (t 920 = 3.854, P = 0.00006), accounting for 1.48% of the variance.
Table 9. Quantitative covariance of non-heligmosomid species richness and worm burdens of non-heligmosomid helminth taxa with heligmosomid worm burdens, as reflected in the linear regressions between the residuals of full-factorial, best-fit or minimum sufficient models (as in Table 8). The explanatory variable in all models is the residuals of the MSM for heligmosomid worm burdens
a The sample size for all = 922.
b The R 2 is for the regression, followed by the gradient of the slope and the SEM of this estimate.
c Student’s t, degrees of freedom (DoF) and P-values (1-tailed) for test of whether the slopes differ from zero, since our prediction in all cases was for a positive slope.
d Species richness of non-heligmosomid helminths. For this dependent variable the MSM was based on Poisson errors. All the others were based on models with negative binomial error structures.
Does the increasing intensity of heligmosomid burdens covary with increasing intensities of non-heligmosomid helminth taxa?
The same approach was implemented to assess the degree of covariance of heligmosomid burdens with non-heligmosomid helminth, non-heligmosomid nematode, oxyuroid nematode, combined cestode, intestinal and larval cestode worm burdens. In all cases, these relationships were positive and significant (Tables 8; Figure 4), except in the case of oxyuroid nematodes and intestinal cestodes where the relationships were marginally non-significant.
Figure 4. Covariance of the residuals of the Minimum Sufficient Model (MSM) for the heligmosomid burdens with the residuals of Full Factorial Models (FFM) for all non-heligmosomid helminths (A and B) and non-heligmosomid nematode (C and D) worm burdens. The panels on the left (A and C) illustrate the predictions of the relevant models in Table 8, whereas those on the right (B and D) show regressions with FFMs of non-heligmosomid helminth (B; adjusted R 2 = 0.0111) and non-heligmosomid nematode (C; R 2 = 0.0090) worm burdens. The shaded areas show the 95% confidence region. For additional details, see text and Table 9.
In linear models testing covariance between the residuals of best-fit models of these taxa and the residuals of the MSM for heligmosomid nematodes, positive slopes were obtained in all cases (Table 9). However, highly significant covariance was found only for non-heligmosomid helminths (R 2 = 1.11%; Figure 4A,B) and non-heligmosomid nematodes (R 2 = 0.9%; Figure 4C,D), with weaker but significant relationships for the Oxyuroidea (R 2 = 0.27%) and combined cestode worm burden (R 2 = 0.39%) evident. As above, the relationships for intestinal cestodes (R 2 = 0.11%) was non-significant, but using this approach, that for larval cestodes (R 2 = 0.15%) also failed to achieve significance.
The database contained multiple cases of animals that were not infected with specific taxa (see Table 1), and because concentrations of double negatives can generate false positive slopes in models of covariance, we repeated the analysis of covariance of residuals from models of relevant taxa, but this time selecting only voles that were infected with at least 1 heligmosomid worm and had at least 1 worm of the specified taxon. Two examples are illustrated in Figure S2 with statistical analyses provided in the figure legend. Significant positive covariance was found for non-heligmosomid helminths and non-heligmosomid nematodes (Figure S2A,B). For intestinal cestodes (n = 107, r = 0.140, ß = 0.153 ± 0.1055, t 105 = 1.447, P = 0.0754, Adjusted R 2 = 0.0102) and for A. tianjinensis n = 231, ß = 0.090 ± 0.0609, t 229 = 1.479, P = 0.0703, Adjusted R 2 = 0.0051) covariance was close to significance, but for all cestodes and larval cestodes covariances were not significant. With the exception of the model for larval cestodes, all revealed positive slopes, but none explained more than 0.9% of variance.
In a further approach, the covariance distribution test of Haukisalmi and Henttonen (Reference Haukisalmi and Henttonen1998) was applied to the data. All the covariances of the specified taxa with heligmosomid worm burdens (Table 10) were positive and most retained significance when subset grouping (the 4 factors) was taken into account, although P values were generally lower. However, covariance of heligmosomids with intestinal cestodes, and with the most prevalent species, C. henttoneni, lost significance when we controlled for the 4 factors (subsample constraints).
Table 10. Covariance of heligmosomid worm burdens with other taxa, as detected by the covariance distribution test of Haukisalmi and Henttonen (Reference Haukisalmi and Henttonen1998)
a Significance (P values) of the test on data without subsamples being taken into account, and with subsample constraints taken into account (i.e. with YEAR, SITE, SEX and AGE accounted for). Probability was calculated by comparison of observed data to 5000 randomizations of the data.
b Covariance, as calculated in the covariance distribution test of Haukisalmi and Henttonen (Reference Haukisalmi and Henttonen1998).
c Probability (P).
Discussion
In this paper, we have focused specifically on heligmosomid nematodes and the consequences of their presence for other species of helminths concurrently residing in bank voles. Our study therefore differs from many others seeking evidence for interactions between helminth species in wild caught rodents because our study was based on a hypothesis derived from laboratory experimental studies, and on its clear predictions.
In Western Europe, bank voles are typically subject to infection by 2 closely related heligmosomid species, Heligmosomum mixtum and Heligmosomoides glareoli (Haukisalmi and Henttonen, Reference Haukisalmi and Henttonen1993a, Reference Haukisalmi and Henttonen1993b; Grzybek et al., Reference Grzybek, Bajer, Bednarska, Alsarraf, Behnke-Borowczyk, Harris, Price, Brown, Osborne, Siński and Behnke2015), and although these are classified in different genera within the family Heligmosomidae, both on morphological grounds and genetically, they are closely related (Durette-Desset, Reference Durette-Desset1985; Zaleśny et al., Reference Zaleśny, Hildebrand, Paziewska-Harris, Behnke and Harris2014; Harris et al., Reference Harris, Zaleśny, Hildebrand, Paziewska-Harris, Behnke, Tkach, Hwang and Kinsella2015). They are also closely related to H. bakeri, a species that has been used as a laboratory model of chronic nematode infections for over 60 years and to H. polygyrus, a parasite of wood mice in the W. Palearctic (Durette-Desset, Reference Durette-Desset1985; Zaleśny et al., Reference Zaleśny, Hildebrand, Paziewska-Harris, Behnke and Harris2014; Harris et al., Reference Harris, Zaleśny, Hildebrand, Paziewska-Harris, Behnke, Tkach, Hwang and Kinsella2015), that also gives rise to long-lasting chronic infections in its host (Gregory et al., Reference Gregory, Keymer and Clarke1990). Given the now well-established immunomodulatory properties of H. bakeri (Behnke et al., Reference Behnke, Hannah and Pritchard1983; Hewitson et al., Reference Hewitson, Harcus, Murray, van Agtmaal, Filbey, Grainger, Bridgett, Blaxter, Ashton, Ashford, Curwen, Wilson, Dowle and Maizels2011; Maizels et al., Reference Maizels, Hewitson, Murray, Harcus, Dayer, Filbey, Grainger, McSorley, Reynolds and Smith2012; Zhang et al., Reference Zhang, Hardy, Kapita, Hall, Arbeeva, Campbell, Urban, Belkaid, Nagler and Iweala2023; Lynch and Grencis, Reference Lynch and Grencis2025), and evidence that the presence of H. polygyrus in wood mice is associated with increased species richness of other helminths (Behnke et al., Reference Behnke, Gilbert, Abu-Madi and Lewis2005, Reference Behnke, Eira, Rogan, Gilbert, Torres, Miquel and Lewis2009), we hypothesized that bank voles carrying heligmosomid nematodes should also be subjected to depressed intestinal resistance to infections arising from the long-term survival strategy of these nematodes. On this basis, we predicted that heligmosomid-infected bank voles are more likely to also carry other species of helminths and heavier infections with the latter, compared to bank voles without current heligmosomid infections. Moreover, we expected to see positive covariance between appropriately adjusted heligmosomid worm burdens (i.e. with both intrinsic and extrinsic factors taken into account) and those of other taxa. All these predictions are clearly supported by our results and associated data analysis.
The strongest support for our hypotheses was evident when we examined the relationships of heligmososmid nematodes with all other helminths, i.e. at the highest taxonomic level for the latter, the non-heligmosomid taxon with the highest prevalence in bank voles (62.9%). In heligmosomid carrying bank voles, prevalence of non-heligmosomid helminths was significantly higher (although marginally context dependent), abundance of non-heligmosomid helminths was significantly higher and evidence for quantitative covariance was found in GLMs of the former with heligmosomid worm burdens (fitted as a covariate), in linear regressions of both taxa, after control for intrinsic and extrinsic factors, and in the covariance distribution test, in which context was controlled for (i.e. subset constraints). Moreover, adjusted worm burdens covaried positively and significantly, when the analyses were confined to voles that carried at least 1 specimen of each taxon. Much the same was found when the analyses were confined to non-heligmosomid nematodes (essentially 5 species combined in this taxon, with a collective prevalence of 54.6%; see Table 1), although the relationships were weaker, and the abundance of non-heligmosomid nematodes did not achieve significance when intrinsic and extrinsic factors were taken into account.
The prevalence and abundance of cestodes (all species combined) was also significantly higher in heligmosomid infected bank voles in all the categories of analysis that were undertaken, despite the much lower prevalence of 20.6%. However, quantitative analyses based on GLMs, all proved difficult to implement and we resorted to fitting best-fit reduced models. None of the tests yielded significance for intestinal strobilate worms, for which there was a high preponderance of animals carrying just a single worm (single worm infections were recorded in 95 of the 150 bank voles infected with intestinal cestodes).
In none of our models for non-heligmosomid taxa, even those showing the highest level of significance, did the presence of heligmosomid nematodes account for more than 1.0% of the variance in other taxa. The most marked effects were evident when we fitted species richness of non-heligmosomid helminths (i.e. by summing the number of species present in bank voles, rather than that of individual worms). For covariance of non-heligmosomid helminth species richness with heligmosomids, and after adjusting for the presence of intrinsic and extrinsic factors, our models accounted for up to 1.48% of variance.
Previous studies have generally found associations of intestinal helminths in rodents to be difficult to detect, rare and mostly context dependent. With respect to heligmosomids, Mongomery and Montgomery (Reference Montgomery and Montgomery1990) found significant positive but mostly weak co-occurrence of H. polygyrus (referred to as N. dubius in their paper) with Corrigia vitta and Syphacia stroma. Although Lewis et al. (Reference Lewis, Morley and Behnke2023) recorded an arithmetically higher prevalence of C. vitta in H. polygyrus-infected wood mice compared to those without the latter species, this was not a significant difference, and there was no quantitative interaction between these species.
Among the earlier attempts to identify whether specific helminths co-occur together more frequently than not, Kisielewska (Reference Kisielewska1970) found that H. halli (but probably more likely H. glareoli or H. mixtum) co-occurred with C. pusilla (C. henttoneni) in bank voles in Eastern Poland. Behnke et al. (Reference Behnke, Gilbert, Abu-Madi and Lewis2005) reported significantly higher prevalence of C. pusilla in wood mice carrying H. polygyrus compared to those without the latter species, and significant positive correlation of worm burdens, even after controlling for context. However, while Lewis et al. (Reference Lewis, Morley and Behnke2023) found numerically higher prevalence of C. pusilla in H. polygyrus-carrying wood mice compared to those without the latter species, this was not a significant difference when context was not accounted for. With context taken into account, Lewis et al. (Reference Lewis, Morley and Behnke2023) found that there was a significant interaction with age; prevalence of C. pusilla was similar among H. polygrus-infected juvenile (7.9%) and adult mice (9.4%), but among mice without H. polygrus prevalence of C. pusilla was higher in adults (10.7%) compared to juveniles (0%). Moreover, in analyses confined to mice carrying both H. polygyrus and C. pusilla, there was a weak positive significant correlation of worm burdens, and this remained significant when the data were controlled for both intrinsic and extrinsic factors.
Another approach used in earlier studies was to assess the effect of H. polygyrus on prevalence and abundance of other helminths (i.e. as reflected in combined worm counts and species richness of non-heligmosomid helminths). Behnke et al. (Reference Behnke, Gilbert, Abu-Madi and Lewis2005) reported significant correlation of H. polygyrus worm burdens with species richness of other helminth species both after fitting raw data and data adjusted for context, but not with the total worm burdens of the other species. In a subsequent study based on analysis of 3 independent data sets (2 localities in England and a dataset based on 6 sites in Portugal), Behnke et al. (Reference Behnke, Eira, Rogan, Gilbert, Torres, Miquel and Lewis2009) found much the same, higher prevalence of other helminths and higher non-heligmosomid species richness in mice infected with H. polygyrus. In contrast to these studies and the current work, Lewis et al. (Reference Lewis, Morley and Behnke2023) found that while prevalence of other helminths (all non-heligmosomids combined) and non-heligmosomid species richness increased with increasing worm burdens of H. polygyrus, the significance was lost when the data were adjusted for context.
Where found to be significant in previous studies, associations between helminth species and correlations in quantitative data have generally been found to be weak and often context dependent. This may be partly due to sample sizes that are too small to detect subtle effects, especially in data sub-sets, and the use of inappropriate statistical methods. Moreover, there are many reasons why a process of functional importance might nonetheless appear to be of small effect size in cross-sectional correlations. This could be due to masking by any of the confounding processes arising from correlated susceptibility or correlated exposure considered above (see Introduction). Additionally, the statistical adjustment for variables such as sex or size or season or locality, although highly desirable when attempting to increase the inference that 1 species might influence another (the key interest here), may also partition out a genuine signal that is collinear with the adjusted variable(s). For example, adjustment for size may mask facilitatory interactions, as these might increase in magnitude as hosts age, grow and are exposed to more parasites (with age, size and parasite exposure likely to be collinear to some extent).
Some have argued that conclusions of the occurrence of quantitative interactions among multiple parasite species in host populations may be unreliable when derived from cross-sectional studies at single time points and associated worm burdens at autopsy (Poulin, Reference Poulin2005) and are unlikely to detect interactions accurately (Knowles et al., Reference Knowles, Fenton, Petchey, Jones, Barber and Pedersen2013; Fenton et al., Reference Fenton, Knowles, Petchey and Pedersen2014; Sweeny et al., Reference Sweeny, Albery, Venkatesan, Fenton and Pederson2021; Barger, Reference Barger2023). Nevertheless, when similar relationships can be detected in different data subsets, collected at different sites in different years, as in the current work, some confidence in their repeatability and hence robustness is warranted (Behnke et al., Reference Behnke, Jackson, Gilbert, Mohallal and Bajer2024). Moreover, where there is consistency with known biological properties (such as the well-known immunosuppressive properties of heligmosomids leading to an expectation of positive interaction and hence positive association) then a stronger causal inference can be sustained.
Haukisalmi and Henttonen (Reference Haukisalmi and Henttonen1993a) concluded that helminth infections in bank vole from their study sites in Finland were mostly independently aggregated in hosts, but that rarely occurring species were restricted to hosts with concurrent common species. To some degree our data support this conclusion, both in that all the significant associations and correlations we detected were weak and in that cestodes, the taxon with the lowest prevalence in our data, showed higher prevalence and abundance in heligmosomid infected voles. Lewis et al. (Reference Lewis, Morley and Behnke2023) and Behnke et al. (Reference Behnke, Gilbert, Abu-Madi and Lewis2005) both also concluded that quantitative interactions are weak, context dependent, perhaps reflecting assemblages of helminth species, rather than interacting communities.
Our analyses demonstrate the powerful influence of controlling for both intrinsic and extrinsic factors, since when we failed to fit these explanatory factors, all of our analyses were highly significant. While it is necessary to carry out such adjustments to increase the inference of interactions between parasites, nonetheless we note that this may also be masking genuine signals of interaction, and reducing their observed effect size, when these are collinear with the adjusted variables. Taken together, our analyses support our prediction that bank voles carrying adult heligmosomid worms are likely to be infected also with other helminths, and that to some extent this relationship is quantitative. In no small measure, the large sample size of our dataset enabled the detection of the subtle (confounder-adjusted) effect of heligmosomid nematodes on other concurrently residing helminth taxa. In this respect, our current findings concur with earlier reports seeking evidence for interactions between helminths in wild caught rodents. Behnke et al. (Reference Behnke, Gilbert, Abu-Madi and Lewis2005) concluded that while associations and interactions between helminths in wood mice in S. England could be detected in cross-sectional surveys, they were highly context dependent and unlikely to play a dominant role among processes involved in structuring helminth component communities in this host species. However, in this study and in a subsequent work analysing data from wood mice from sites in England and Portugal (Behnke et al., Reference Behnke, Eira, Rogan, Gilbert, Torres, Miquel and Lewis2009), the likelihood of being infected with non-heligmosomid helminth species, as reflected in non-heligmosomid species richness, was more likely in heligmosomid infected mice and a highly predictable element of co-infections in wood mice. Our current results support this view and extend it to infections in bank voles. Indeed, given the consonance of these results in gastrointestinal communities containing heligmosomid species, and similar results observed in spirurid-dominated communities (Behnke et al., Reference Behnke, Jackson, Gilbert, Mohallal and Bajer2024), it might be predicted that many vertebrate gastrointestinal communities are structured by immunosuppressive gastrointestinal nematodes that act as ecosystem engineers (Wright et al., Reference Wright, Jones and Flecker2002), impacting species richness and heterogeneity in community assemblages.
Supplementary material
The supplementary material for this article can be found at https://doi.org/10.1017/S0031182025101376.
Acknowledgements
We thank Mgrs Grzegorz Górecki and Anna Zaborowska, the wardens of the field station at Urwitałt, for the use of the University of Warsaw’s field station facilities. We wish to acknowledge the support of the forestry departments responsible for the woodland sites utilized in our study (Nadleśnictwa Giżycko, Mikołajki and Orzysz). We are grateful to Em. Prof F. Gilbert for his advice on earlier versions of this manuscript. Finally, we thank all the undergraduate student helpers from both NU and WU who helped with fieldwork during the 4 expeditions to our study sites.
Author contributions
The fieldwork was conceived, planned and implemented by JMB and AB, and supported by MG, MA and JBB and undergraduate students. The data were analysed by JMB and JAJ, and the paper was written by JMB, JAJ, MG and AB.
Financial support
We thank the Universities of Nottingham and Warsaw for enabling this study through financial support for travel and consumables. Our work was supported additionally by travel grants from the Royal Society (JMB), the British Ecological Society (JMB), the Grabowski Fund (JMB), the Polish State Committee for Scientific Research (KBN, AB), and by the KBN and British Council’s Young Scientist Programme (AB).
Competing interests
None
Ethical standards
The data analysed herein were derived from fieldwork carried out in strict accordance with the recommendations in the Guide for the Care and use of Laboratory Animals of the National Ethics Committee for Animal Experimentation. Our project was approved by the First Warsaw Local Ethics Committee for Animal Experimentation (protocol no 73/2010), which also has overarching responsibility for fieldwork involving the trapping and culling of wild vertebrates for scientific purposes.
Open access
