Identifying sources of variation in parasite aggregation

Sample collection and the symbiont assemblage

Samples of 79–102 ptarmigan were collected by hunting in each autumn of 12 years (collections organized by the Icelandic Institute of Natural History) to ensure sufficient sample size for reliable estimates of infection measures (Wilson et al., 2002). Furthermore, ca. 20 adult male, 20 adult female, and 60 juveniles (30 of both sexes) were collected each year which minimized any effects of differential representation of age or sex classes of birds between years on infection measures (see Nielsen et al., 2020). In eight of 12 years, sample size was at least 100. Processing of samples and standardized protocols to enumerate parasites are detailed elsewhere (Skírnisson, Thorarinsdottir & Nielsen, 2012; Stenkewitz et al., 2016). However, some salient points bear repeating. Only parasites identifiable to species and that were not rare (<2% prevalence) were included in analyses. In total, our study included four taxa of symbionts, comprising 13 species: four insects (three lice and one fly), five mites, two species of coccidians, and two nematodes.

All parasite species were observed in all years, except for Trichostrongylus tenuis which was absent in 2007, 2008, and 2012, and Mironovia lagopus, for which there were no abundance measures greater than zero in 2006 (though this mite was still noted as present in some feathers in this year, it was not represented in samples based on filtered vacuuming).

The coccidians included the specialist Eimeria muta whose annual prevalence was shown to be correlated with time-lagged population cycles of ptarmigan (Stenkewitz et al., 2016), and the less common E. rjupa. Coccidians are often considered microparasites because they multiply within their hosts during the endogenous phase of their life cycle, though elsewhere they may be grouped with macroparasites because they cause intensity-dependent pathology (Norton, 1967; Trigg, 1967). The two other endoparasites were nematodes—Capillaria caudinflata and Trichostrongylus tenuis—which live in the small intestine and ceca, respectively. C. caudinflata is associated with overwintering mortality of juvenile birds (Stenkewitz et al., 2016), but the far less abundant T. tenuis was not associated with serious disease outcomes in these ptarmigan (Stenkewitz et al., 2016). In other studies, T. tenuis is much more prevalent and numerous in wild grouse and it is inversely related to fecundity (e.g., Dobson & Hudson, 1992).

The ectoparasites included skin mites Metamicrolichus islandicus (known to cause mange) and Myialges borealis, the two often appearing in co-infections (Stenkewitz et al., 2015). The quill mite (Mironovia lagopus), another ectoparasite, was infrequently sampled, whereas the two remaining mite species Strelkoviacarus holoapsis and the wing mite Tetraolichus lagopi were often sampled on birds (T. lagopi nearly always showed 100% prevalence in samples). These two species are suspected of being paraphages or mutualistic symbionts rather than true parasites (Stenkewitz, 2017; Doña et al., 2019). The three ectoparasitic lice included the chewing louse Amyrsidea lagopi which is thought to be harmful to birds by damaging plumage (Stenkewitz et al., 2016) and the highly specialized Goniodes lagopi and Lagopoecus affinis which feed on keratin. The fourth insect species, the hippoboscid fly Ornithomya chloropus, is an ectoparasite that sucks blood from its bird host and which acts as a vector of mites of ptarmigan (Guðmundsson, Skírnisson & Nielsen, 2021).

Statistical methods

Our chosen measure for quantifying aggregation, Poulin’s D, is calculated as: (1)${\displaystyle D=1-\frac{2\sum _{i=1}^N\left(\sum _{j=1}^i{x}_j\right)}{\bar{x}N\left(N+1\right)}}$

where N is the total number of hosts in the sample, $\bar{x}$ is the average number of parasites, and x is the number of parasites infecting host j, with hosts ordered from least to most heavily infected (Poulin, 1993). D is thus constrained to fall between approximately zero (minimum aggregation; an even distribution of parasites) and one (maximum aggregation; all parasites infecting a single host) regardless of parasite mean abundance. The measure therefore lends itself to analysis of aggregation across samples, studies, and/or species.

Modeling was conducted in a Bayesian framework using Stan software (Stan Development Team, 2019), which implements an adaptive method of Hamiltonian Monte Carlo (HMC) sampling of posterior distributions via the No-U-Turn sampler (NUTS; Hoffman & Gelman, 2014; Stan Development Team, 2019). The package ‘cmdstanr’ (Gabry & Cešnovar, 2021) was used as an interface between Stan and the statistical programming language R (Version 3.6.1; R Core Team, 2019).

We modeled D using beta regressions implemented as generalized linear mixed models (GLMMs; beta-distributed response, with a logit link function) because D is constrained between zero and one. Use of beta regression to model parasite aggregation appears unprecedented but provides an effective tool for assessing determinants/correlates of aggregation. Additionally, the ease of implementing custom hierarchical model structure afforded by Bayesian modeling when using software like Stan also allows effective information “sharing” across predictor factor levels through the incorporation of random effects, and the specification of even weakly-informative priors may permit the fitting of otherwise prohibitively complex, though biologically justified, models. A Bayesian approach could directly incorporate uncertainty in the response and/or predictors, a further benefit to analyzing predictors of aggregation across studies. We nonetheless recommend beta regression as an effective tool for modeling aggregation as measured by D regardless of whether the implementation is Bayesian, as it is here, or not.

Beta distributions are generally described by two parameters, though there are different choices of parameterizations. Here, we used the mean (µ) and variance (ν) parameterization rather than the more common α and β parameterization (to switch between the two: α = µν, and β = (1 − µ) ν). This choice allowed us to directly estimate mean values for D between zero and one, following the linear model. When ν is lower than two, the probability density of the beta distribution begins getting “pushed” towards the extremes of zero and one, rather than being spread (diminishing) around the mean (McElreath, 2020). Therefore, the implemented beta regression sets a lower limit of two on the variance parameter ν via a half-Cauchy prior to ensure that modeled aggregation levels are concentrated around the means (McElreath, 2020).

The predictors of aggregation in fitted models included at least one intercept term, possible slope terms relating sample mean abundance to aggregation, and a random effect of sampling year. Often, mean abundance is expected to relate negatively to aggregation, including when measured using Poulin’s D (Poulin, 1993; Shaw, Grenfell & Dobson, 1998; Johnson & Wilber, 2017). Mean abundance was log₂-transformed due to high overall positive skew, and then centered and standardized before modeling. In addition to a single-intercept (null) model, we considered candidate models with intercept (θ) terms indexed either by parasite species, taxonomic grouping (insect, mite, coccidian, nematode), or the ectoparasite-endoparasite dichotomy. The potential mean abundance relationship, added to relevant models as a slope term (β), was either included as a single parameter or indexed by the same potential groupings as listed above (i.e., species, broader taxonomic grouping, or ecto-/endoparasite). Therefore, we compared models to test whether parasite biology, either on its own or through its modulation of potential aggregation-mean abundance relationships, was an important predictor of aggregation. A single intercept, single slope candidate model follows the form: ${\displaystyle D_i\sim \text{Beta}\left({\mu }_i,\nu \right)}$ ${\displaystyle \text{logit}\left({\mu }_i\right)=\theta +\beta \times m_i+{\lambda }_{YEAR\left[i\right]}}$

where m_i is the mean abundance of sample i, and λ _YEAR[i] is the random effect of the relevant year. Selected weakly-informative priors for the θ, β, ν, and λ parameters were: ${\displaystyle \theta \sim \text{Normal}\left(0,1.5\right)}$ ${\displaystyle \beta \sim \text{Normal}\left(0,2.2\right)}$ ${\displaystyle \nu \sim \text{Half-Cauchy}\left(2,2.5\right)}$

${\lambda }_j\sim \text{Normal}\left(0,{\sigma }_{\lambda }\right)$ for j = 1..12 ${\displaystyle {\sigma }_{\lambda }\sim \text{Exponential}\left(1\right)}$

Here, the notation Normal(x, y) refers to a normal distribution with mean x and standard deviation y; Half-Cauchy(a, b) refers to a half-Cauchy distribution with location parameter (and lower limit) a and scale parameter b; and Exponential(c) refers to an exponential distribution with rate parameter c (mean = 1/c). Prior predictive simulations exploring the relevance of the selected intercept and slope priors are presented in Article S1. The exponential prior with rate parameter equal to one on the standard deviation term follows McElreath (2020).

When the intercept or slope parameters were indexed by factors with more than two groups (i.e., either parasite species or broader parasite taxonomic grouping), these were also included as random effects. For example, the linear model and θ and β priors above would, in the candidate species-specific intercept and group-specific slope model, instead become: ${\displaystyle \text{logit}\left({\mu }_i\right)={\theta }_{SPECIES\left[i\right]}+{\beta }_{GROUP\left[i\right]}\times m_i+{\lambda }_{YEAR\left[i\right]}}$

${\theta }_k\sim \text{Normal}\left(\theta ,{\sigma }_{\theta }\right)$    for k = 1..13

${\beta }_n\sim \text{Normal}\left(\beta ,{\sigma }_{\beta }\right)$    for n = 1..4 ${\displaystyle \theta \sim \text{Normal}\left(0,1.5\right)}$ ${\displaystyle \beta \sim \text{Normal}\left(0,2.2\right)}$ ${\displaystyle {\sigma }_{\theta }\sim \text{Exponential}\left(1\right)}$ ${\displaystyle {\sigma }_{\beta }\sim \text{Exponential}\left(1\right)}$

It is possible that intercept and slope parameters of a linear model may covary when indexed by the same grouping variable: higher intercepts may relate to higher—or lower—slope terms. Therefore, whenever θ and β random effects were indexed by the same variable, these were both modeled as following a single multivariate normal distribution, as this model structure allows underlying correlation to be incorporated, and to inform parameter estimation (McElreath, 2020). The correlation matrix for θ and β, R, can be factored out of the multivariate normal distribution covariance matrix, S, resulting in this formulation (for the sake of example, given the model wherein both intercepts and slopes are species-specific): ${\displaystyle D_i\sim \text{Beta}\left({\mu }_i,\nu \right)}$ ${\displaystyle \text{logit}\left({\mu }_i\right)={\theta }_{SPECIES\left[i\right]}+{\beta }_{SPECIES\left[i\right]}\times m_i+{\lambda }_{YEAR\left[i\right]}}$ ${\displaystyle \left[\begin{array}{c}\hfill {\theta }_k\hfill \\ \hfill {\beta }_k\hfill \end{array}\right]\sim \text{MVNormal}\left(\left[\begin{array}{c}\hfill \bar{\theta }\hfill \\ \hfill \bar{\beta }\hfill \end{array}\right],S\right)}$ ${\displaystyle S=\left(\begin{array}{cc}\hfill {\sigma }_{\theta }\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill {\sigma }_{\beta }\hfill \end{array}\right)R\left(\begin{array}{cc}\hfill {\sigma }_{\theta }\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill {\sigma }_{\beta }\hfill \end{array}\right)}$ ${\displaystyle \bar{\theta }\sim \text{Normal}\left(0,1.5\right)}$ ${\displaystyle \bar{\beta }\sim \text{Normal}\left(0,2.2\right)}$ ${\displaystyle \nu \sim \text{Half-Cauchy}\left(2,2.5\right)}$ ${\displaystyle {\lambda }_j\sim \text{Normal}\left(0,{\sigma }_{\lambda }\right)}$ ${\displaystyle {\sigma }_{\theta }\sim \text{Exponential}\left(1\right)}$ ${\displaystyle {\sigma }_{\beta }\sim \text{Exponential}\left(1\right)}$ ${\displaystyle {\sigma }_{\lambda }\sim \text{Exponential}\left(1\right)}$ ${\displaystyle R\sim \text{LKJcorr}\text{(2)}.}$

The weakly-informative prior on the correlation matrix R is the Lewandowski-Kurowicka-Joe distribution, with its η parameter set to two, as this prior is skeptical of extreme correlations (i.e., those close to −1 or 1), and has its probability density centred around zero correlation (Lewandowski, Kurowicka & Joe, 2009; Stan Development Team, 2019).

We fit each Bayesian model using three Markov chains, with 2,500 warm-up iterations and 1,500 sampling iterations (total samples = 3 × 1,500 = 4,500). Effective sample sizes for estimated parameters were assessed to ensure accurate estimation (ESS >100 per chain; Vehtari et al., 2021), and $\hat{R}$ convergence diagnostics for all parameters were checked to ensure they were less than 1.01 (an indication that chains had mixed well). The R package ‘shinystan’ was used to help assess model convergence (Gabry, 2018), and functions from the package ‘tidybayes’ provided tools for exploring posterior distributions (Kay, 2021). Posterior predictive checks were conducted to ensure that model estimates were consistent with observed data.

All fitted models were compared by their estimated out-of-sample predictive performance during leave-one-out cross validation, based on estimated log pointwise predictive density (ELPD; higher values indicate more predictive models). This method considers the entire posterior distribution, and relies on Pareto-smoothed importance sampling (PSIS) to calculate estimates (Vehtari, Gelman & Gabry, 2017). ELPD estimates, as well as the estimates of the difference in ELPD between a given model and the overall best-performing model, are provided with calculated standard errors, allowing the assessment of levels of confidence (e.g., if a model has higher ELPD compared to another, but the standard error in that estimated difference is larger than the difference itself, then there is very little confidence that the two models differ meaningfully in out-of-sample predictive performance). Additionally, PSIS calculations provide Pareto k values for each observation, which can indicate when ELPD estimates are unreliable, and/or that particular observations have disproportionate influence on the posterior distribution; k values should be less than 0.7 (Vehtari, Gelman & Gabry, 2017). The ‘loo’ package provided functions for the calculation of ELPDs and their comparisons (Vehtari et al., 2020). All reported credible intervals (CIs) are 95% highest density continuous intervals.