Introduction
Parasite infections are widely known to follow seasonal patterns, with variation in prevalence and abundance of parasites apparent within years (Altizer et al., Reference Altizer, Dobson, Hosseini, Hudson, Pascual and Rohani2006). This seasonality is known to be driven by a number of factors, including climatic effects on parasite development and mortality, host physiology or behaviour, or seasonal variability in the availability of naïve hosts (Haukisalmi et al., Reference Haukisalmi, Henttonen and Tenora1988; Montgomery and Montgomery, Reference Montgomery and Montgomery1988; Haukisalmi and Henttonen, Reference Haukisalmi and Henttonen1990; Halvorsen et al., Reference Halvorsen, Stien, Irvine, Langvatn and Albon1999; Pelletier et al., Reference Pelletier, Page, Ostiguy and Festa-Bianchet2005; Rossin et al., Reference Rossin, Malizia, Timi and Poulin2010; Turner and Getz, Reference Turner and Getz2010; Cizauskas et al., Reference Cizauskas, Turner, Pitts and Getz2015). Much less well understood is how parasite communities change across years, and specifically the extent to which parasites show long-term temporal trends in prevalence (the proportion of the host population that they infect) or abundance (the number of individual parasites harboured per host). Long-term studies of wild animal populations have frequently found evidence for variation and trends in population size across time (Clutton-Brock et al., Reference Clutton-Brock, Grenfell, Coulson, Maccoll, Illius, Forchhammer, Wilson, Lindstrom, Crawley, Albon, Clutton-Brock and Pemberton2004; Armitage, Reference Armitage2012; Andreassen et al., Reference Andreassen, Sundell, Ecke, Halle, Haapakoski, Henttonen, Huitu, Jacob, Johnsen, Koskela, Luque-Larena, Lecomte, Leirs, Mariën, Neby, Rätti, Sievert, Singleton, Van Cann, Vanden Broecke and Ylönen2021; Crawley et al., Reference Crawley, Pakeman, Albon, Pilkington, Stevenson, Morrissey, Jones, Allan, Bento, Hipperson, Asefa and Pemberton2021), and trends in phenological traits such as timing of breeding (Charmantier et al., Reference Charmantier, Mccleery, Cole, Perrins, Kruuk and Sheldon2008; Bonnet et al., Reference Bonnet, Morrissey, Morris, Morris, Clutton-Brock, Pemberton and Kruuk2019), morphological traits (Ozgul et al., Reference Ozgul, Tuljapurkar, Benton, Pemberton, Clutton-Brock and Coulson2009) and sexually selected traits (Evans and Gustafsson, Reference Evans and Gustafsson2017). Given these long-term trends in host dynamics, one would expect similar temporal effects to apply to their parasites, but the link between long-term host and parasite dynamics is not well-understood. Similar long-term monitoring of the prevalence and abundance of parasites is important for addressing their potential impact on host life histories, immunological profiles and their influence as a selective force on host populations (Hudson et al., Reference Hudson, Dobson and Newborn1998; Babayan et al., Reference Babayan, Liu, Hamilton, Kilbride, Rynkiewicz, Clerc and Pedersen2018). Such long-term changes may be linked to a number of drivers, including climatic variation (Haukisalmi and Henttonen, Reference Haukisalmi and Henttonen1990; Huntley et al., Reference Huntley, Fürsich, Alberti, Hethke and Liu2014), since the development of transmission stages and vectors may be temperature- and moisture-dependent (Nordenfors et al., Reference Nordenfors, Höglund and Uggla1999; O'Connor et al., Reference O'Connor, Walkden-Brown and Kahn2006; Rose et al., Reference Rose, Hoar, Kutz and Morgan2014); changes in host demography, for example, if more vulnerable host age or sex classes change in relative abundance; or changes in selection on host immune function and/or related traits. Only long-term studies offer the power to identify temporal trends in parasite prevalence and abundance and determine their drivers, and indeed the power of studies to make inferences about the impact of diseases on host health, fitness and population dynamics is positively associated with their duration (Barroso et al., Reference Barroso, Acevedo and Vicente2021).
The prevalence and abundance of many common and rare parasite species have been shown to be remarkably stable across years in studies of wild rodent populations (Haukisalmi et al., Reference Haukisalmi, Henttonen and Tenora1988; Bajer et al., Reference Bajer, Behnke, PawełCzyk, KuliŚ, Sereda and Sinski2004; Behnke et al., Reference Behnke, Bajer, Harris, Newington, Pidgeon, Rowland, Sheriff, Kulis-Malkowska, Sinski, Gilbert and Barnard2008; Knowles et al., Reference Knowles, Fenton, Petchey, Jones, Barber and Pedersen2013; Behnke et al., Reference Behnke, Bajer, Behnke-Borowczyk, Clisham, Gilbert, Glover, Jeffery, Kirk, Mierzejewska, Mills, Mohallal, Padget, Wainer and Zalat2018). Where temporal variation does occur, it appears to consist of fluctuations between years rather than consistent trends across time (Behnke et al., Reference Behnke, Lewis, Zain and Gilbert1999; Haukisalmi and Henttonen, Reference Haukisalmi and Henttonen2000; Irvine et al., Reference Irvine, Stien, Halvorsen, Langvatn and Albon2000; Lachish et al., Reference Lachish, Knowles, Alves, Wood and Sheldon2011; Grzybek et al., Reference Grzybek, Bajer, Bednarska, Al-Sarraf, Behnke-Borowczyk, Harris, Price, Brown, Osborne, Siński and Behnke2015; Peet et al., Reference Peet, Kirk and Behnke2022). In some cases, however, both colonization of hosts by novel species and extinction of previously prevalent species have been reported (Kennedy et al., Reference Kennedy, Shears and Shears2001; Lyndon and Kennedy, Reference Lyndon and Kennedy2001; Grzybek et al., Reference Grzybek, Bajer, Bednarska, Al-Sarraf, Behnke-Borowczyk, Harris, Price, Brown, Osborne, Siński and Behnke2015). The majority of these studies take place across a maximum of 10 years, and use destructive sampling of a relatively small number of hosts in order to identify and enumerate the species present. Two recent studies have, however, found evidence for consistent temporal trends in parasite prevalence and abundance. In a 26-year study of the helminth fauna of a population of wood mice (Apodemus sylvaticus) in NE England, the parasite community was dominated by the nematode Heligmosomoides polygyrus, which showed a marked decline in prevalence and abundance across the study period, patterns that were not observed for other parasite species in the system (Behnke et al., Reference Behnke, Rogan, Craig, Jackson and Hide2021). Further, a similar decline in H. polygyrus fecal egg count (FEC) was observed over 6 years in A. sylvaticus populations in NW England (Sweeny et al., Reference Sweeny, Albery, Venkatesan, Fenton and Pedersen2021). These studies may have been able to detect these temporal trends because of their sampling regimes, which occurred over a long time scale (Behnke et al., Reference Behnke, Rogan, Craig, Jackson and Hide2021) or were non-destructive (Sweeny et al., Reference Sweeny, Albery, Venkatesan, Fenton and Pedersen2021), both leading to a larger sample size and hence potentially greater power to detect consistent temporal trends. Perhaps most importantly, both studies returned to the same sampling sites every year, while many destructive sampling studies sample intermittently, leaving gaps in the data in order to enable host populations to recover. The lack of evidence for temporal trends may therefore stem from a combination of relatively small sample sizes, short time spans and intermittent sampling. Studies sampling sufficient individuals consistently at the same sites across a longer period would therefore seem to offer the greatest ability to discern temporal variation.
In addition to changes in the prevalence and mean abundance of parasites, temporal changes in the distribution of parasites among individuals across time may also be important (Grenfell et al., Reference Grenfell, Wilson, Isham, Boyd and Dietz1995; Haukisalmi and Henttonen, Reference Haukisalmi and Henttonen1999). Generally, only hosts with the highest parasite burdens will experience fitness effects (although tolerance of infection may play a role in mitigating this), and the number and proportion of hosts to which this applies will depend on the distribution of parasites across the host population (Poulin and Vickery, Reference Poulin and Vickery1993). Parasites tend to follow a negative binomial distribution among hosts, the shape of which dictates that the majority of hosts have relatively few parasites and a long ‘tail’ of a few hosts with high parasite burdens (Haukisalmi, Reference Haukisalmi1986; Woolhouse et al., Reference Woolhouse, Dye, Etard, Smith, Charlwood, Garnett, Hagan, Hii, Ndhlovu, Quinnell, Watts, Chandiwana and Anderson1997; Wilson et al., Reference Wilson, Bjornstad, Dobson, Merler, Poglayen, Randolph, Read, Skorping, Hudson, Rizzoli, Grenfell, Heesterbeek and Dobson2002). The negative binomial distribution is characterized by the aggregation parameter k, and as k increases, the distribution becomes essentially Poisson or random. Thus, as k increases a greater proportion of the hosts fall into the tail of the distribution with high burden and parasites may be more likely to act as a regulatory factor for hosts (Anderson and May, Reference Anderson and May1978; May and Anderson, Reference May and Anderson1978). Variation in k across time may be generated through a number of processes that may serve to reduce aggregation, such as parasite mortality and parasite-induced host mortality, and increase it, such as heterogeneity in host resistance to infection as a result of genes or the environment (Anderson and Gordon, Reference Anderson and Gordon1982). Determining if, for example, parasite distributions are becoming more or less aggregated across time may help us to determine which of these processes are occurring in a given host–parasite system, and lead us to identify the ecological or evolutionary drivers.
In this study, we determined how parasite prevalence and abundance have changed over 31 years in the wild Soay sheep population living in the St. Kilda archipelago, NW Scotland. Since the study of this population began in 1985, longitudinal data have been collected on the prevalence and abundance of several gastrointestinal parasite species, including apicomplexan and helminth parasites, with an abundance of associated data on host phenotype and environmental conditions (Clutton-Brock and Pemberton, Reference Clutton-Brock and Pemberton2004b). The sheep population has experienced significant temporal change during this time, with average body size declining (Ozgul et al., Reference Ozgul, Tuljapurkar, Benton, Pemberton, Clutton-Brock and Coulson2009) and population size increasing along with grassland productivity (Crawley et al., Reference Crawley, Pakeman, Albon, Pilkington, Stevenson, Morrissey, Jones, Allan, Bento, Hipperson, Asefa and Pemberton2021). The known impacts of gastrointestinal parasites on condition and fitness in the Soay sheep (Gulland, Reference Gulland1992; Craig et al., Reference Craig, Tempest, Pilkington and Pemberton2008; Hayward et al., Reference Hayward, Wilson, Pilkington, Clutton-Brock, Pemberton and Kruuk2011) suggest that parasites may play a role in these dynamics, and the long-term nature of the study and abundance of metadata make the Soay sheep an excellent system in which to study long-term temporal trends in parasite infection and its drivers. Here, we take the first step by describing temporal changes in 4 parasite taxa across 31 years.
Materials and methods
Study population and data collection
We collected data from the wild population of Soay sheep living on Hirta in the St. Kilda archipelago, NW Scotland (57°49′N, 08°34′W). The population are descendants of early European domestic sheep that were introduced to the island of Soay several thousand years ago (Clutton-Brock and Pemberton, Reference Clutton-Brock, Pemberton, Clutton Brock and Pemberton2004); in 1934–1935, 108 sheep were moved from Soay onto Hirta and the population has since expanded across the island. The population living in the Village Bay area of Hirta has been intensively studied since 1985, with data collected on births, deaths, morphometrics, genetics, parasite infections and abundance, spatial use and environmental variation (Clutton-Brock and Pemberton, Reference Clutton-Brock, Pemberton, Clutton Brock and Pemberton2004). Summer population density (PD) is estimated by 10 censuses of the population in August of each year; any animal seen in any census is included in the population estimate (Fig. S1).
Since 1988, fecal samples have been collected and gastrointestinal parasite infection quantified as FEC, using a modified version of the McMaster egg-counting technique, accurate to 100 eggs per gram of feces (M.A.F.F., 1986). Around 85% of the parasite data included in this study were collected by one of the authors (JGP). The McMaster egg-counting technique has been shown to be a good index of actual parasite burden in Soay sheep, both on St. Kilda and elsewhere (Wilson et al., Reference Wilson, Grenfell, Pilkington, Boyd, Gulland, Clutton-Brock and Pemberton2004). The present study exploits data collected in August of each year when animals were captured for weighing, skeletal measurement and blood and fecal sampling. Fecal samples were collected per rectum and stored at 4°C until processing, all of which took place on Hirta. The parasite fauna in Soay sheep is very similar to that of domestic sheep on the mainland, with the notable exception of the nematode Haemonchus contortus. Among the principal taxa are strongyle nematodes, a group of 6 species of the gastrointestinal tract, namely Teladorsagia circumcincta, Trichostrongylus axei, Trichostrongylus vitrinus, Chabertia ovina, Bunostomum trigonocephalum and Strongyloides papillosus. Since the eggs of these species are indistinguishable to the naked eye, they are counted as a single ‘strongyle’ FEC. The sheep also experience infection with 11 species of apicomplexan microparasites of the genus Eimeria, the oocysts of which are again difficult to distinguish by the naked eye and which are therefore incorporated as a single ‘coccidia’ fecal oocyst count (FOC). Strongyles and coccidia are by far the most prevalent and abundant parasite taxa present on St. Kilda (Craig, Reference Craig2005). In addition, eggs of the nematodes Nematodirus spp., Trichuris ovis and Capillaria longipes are enumerated, and the presence or absence of eggs of the cestode Moniezia expansa is recorded. Data on strongyles were collected 1988–2018, while collection for the other taxa commenced in 1993; the exceptionally low prevalence (<1%) of T. ovis and C. longipes meant that we did not include them in any of our analyses (Table 1). Across the history of the project, several experiments have given anthelminthic treatments to a small number of animals in order to determine the effects on parasite burden and particularly survival and reproduction (Gulland, Reference Gulland1992; Craig et al., Reference Craig, Jones, Pilkington and Pemberton2009). Most of these were long-lasting intraruminal boluses that released albendazole across several months and which were shown to be associated with lower worm burden at death in the following spring (Craig et al., Reference Craig, Jones, Pilkington and Pemberton2009). Our analyses aimed to quantify changes in the annual prevalence of 4 of these parasite taxa (strongyles, coccidia, Nematodirus spp., M. expansa), and abundance of strongyles and coccidia across the 31 years of the study.
Number of samples (N), host individuals (IDs) and summary statistics for the parasite taxa analysed
The mean abundance of M. expansa is NA because only presence or absence is noted.
Statistical analyses
All statistical analyses were conducted using R ver. 4.1.2 (R Core Team, 2021).
For each parasite group (strongyles, coccidia, Nematodirus, M. expansa), we first modelled changes in prevalence across time. We fitted generalized additive mixed-effects models (GAMMs) using the R package ‘gamm4’ (Wood and Scheipl, Reference Wood and Scheipl2020). GAMMs fit non-parametric smoothing terms for continuous explanatory variables, meaning that the response variable is not restricted to following parametric forms (e.g. linear, quadratic). To determine the temporal dynamics of infection, we fitted binomial models such that each sample scored 0 if the parasite was absent and 1 if it was present. We then tested a series of models to determine how prevalence changed across the years of the study, and whether changes in prevalence differed between age or sex classes. All models included individual animal identity and year as random effects in order to account for repeated sampling of individuals across years. We also included a categorical fixed effect (‘treatment’) in each model to account for anthelminthic treatments that were given as part of experiments in some years of the long-term study. Animals given an anthelminthic treatment in the year prior to data collection were scored as 1, while animals that did not received a 0; only 133/6873 samples (~2%) were associated with a treatment in our dataset. We took this conservative approach because of the long-lasting effect of treatments on worm burden (Craig et al., Reference Craig, Jones, Pilkington and Pemberton2009). All models also included a 3-level categorical fixed effect for age group, distinguishing lambs (animals aged ~4 months) from yearlings (animals aged ~16 months) and adults (animals aged ~28 months and older), due to the generally higher prevalence and abundance of parasites in younger animals (Wilson et al., Reference Wilson, Grenfell, Pilkington, Boyd, Gulland, Clutton-Brock and Pemberton2004); 37% of our data were from lambs, 15% from yearlings and 48% from adults. All models included sex as a fixed categorical variable with 2 levels (67% female; 33% male), and the interaction between age class and sex.
To this ‘base’ model, we then added fixed effects describing temporal trends as described in detail in Table 2. Given the increase in PD across time (Fig. S1), and the role of host density in driving parasite transmission, accounting for PD in our analyses was essential. In models 1–8, we explored main effects of both year and PD, fitting models with linear and/or non-parametric smooth terms for year and/or PD. Linear models tested for directional trends, while smooth terms tested for non-linear variation across time. In all models, year and PD were both standardized to mean = 0 and s.d. = 1 in order to aid model convergence. In models 9–23, we then introduced interactions between smooth terms for year and/or PD with age class, sex or the combined age/sex category; we used smooth rather than linear terms since the smooth simply reduces to a linear function if a linear association is present. Each of models 9–23 always included smooth terms for both year and PD, ensuring that any effects of either that we observed were independent of the other variable. We compared the Akaike information criterion (AIC) values for each model, and considered models with lower AIC values to be better supported by the data. Models within 2 AIC units of the best model were considered to receive some support from the data, and this is reflected in our presentation of results.
A comparison of models describing temporal trends in the prevalence of strongyles and coccidia in St. Kilda Soay sheep
In the model terms list, ‘s(Year/PD)’ indicates a non-parametric smooth function was fitted to the variable inside parentheses. Note that both year and PD (population density) were standardized to mean = 0 and s.d. = 1 to aid model convergence. ‘ΔAIC’ is the difference in AIC value between the model in question and the model with the lowest AIC value, which is highlighted in bold. Models shaded in grey fall within AIC = 2 of the best-supported model and are considered to receive some support from the data.
For strongyles and coccidia, we then repeated a similar analysis in order to identify temporal patterns of abundance measured as FEC and FOC respectively, once again accounting for PD. For both parasite groups, we used negative binomial models with a log link function in ‘gamm4’, where the variance is modelled as σ 2 = μ + μ 2/k, where μ is the mean and k is the dispersion parameter. The dispersion parameter was estimated from the data and provided to the model as a starting value. We fitted the same 23 models as described above for prevalence and is shown in Table 3.
A comparison of models describing temporal trends in the abundance of strongyles and coccidia in St. Kilda Soay sheep
In the model terms list, ‘s(Year/PD)’ indicates a non-parametric smooth function was fitted to the variable inside parentheses. Note that both year and PD (population density) were standardized to mean = 0 and s.d. = 1 to aid model convergence. ‘ΔAIC’ is the difference in AIC value between the model in question and the model with the lowest AIC value, which is highlighted in bold.
Finally, we tested for changes in 3 parameters describing the distribution of strongyle FEC and coccidian FOC across time and with PD. For each year, we calculated the variance [$\sigma ^2 = \mathop \sum \limits_i ( {x_i-\;\mu } ) ^2/( n-1)$
], the dispersion index (I = σ 2/μ) and the negative binomial aggregation parameter (k = μ 2/σ 2 − μ), where xi, μ and σ are the values for an individual sample i, the population mean and the population standard deviation of FEC/FOC, respectively, for that year. We then fitted a linear model testing for linear and quadratic changes in each of these parameters with both time and PD. We simplified this initial model using Wald F-tests to determine whether quadratic, linear, or neither, associations were supported. We then calculated age class-specific (lamb, yearling, adult) parameters and tested whether changes in the parameters with year or density differed between age classes, and then did the same for sex-specific parameters. If quadratic year or density was supported in the initial model, we included the quadratic term in the age- and sex-specific models; if linear or neither term was supported, we used the linear term in the age- and sex-specific analysis.
Results
Coccidia and strongyles are the most prevalent and abundant parasites infecting the sheep population, with both found in >70% of animals (Table 1). Meanwhile, Nematodirus spp. and M. expansa are found in ~10% of individuals and C. longipes and T. ovis in <1%.
Temporal variation in prevalence
The prevalence of strongyles was influenced by year and PD independently. Of the main effect only models 0–8, the model with smooth terms for both year and PD had the lowest AIC (model 8 ΔAIC = −2.30 compared to model 4 with same terms but no smoothing, Table 2). The model with the lowest AIC overall was model 15, which fitted interactions between the smooth terms for both year and PD with age class (Table 2). The model predicted an approximately linear increase in the prevalence of strongyles in lambs and adults over time, with a curvilinear association in yearlings, predicting an increase in 1988–1995, a decline towards 2010 and then stasis thereafter (Fig. 1A). In addition, the model predicted higher strongyle prevalence at higher host densities, although the shape of this association varied among age classes: the model predicted an approximately linear increase with PD in lambs, while in yearlings and adults the model predicted steep increases in prevalence from low-to-moderate PD and relative stasis at higher PD (Fig. 1B). Models 16 and 17 had ΔAIC < 2 and so were considered to receive some support; they both supported the interaction between year and age, but also provided support for interactions between PD and sex (model 16) and between PD, sex and age (model 17). These results are presented in Fig. S2.
Variation in strongyle prevalence across both (A) the study period and (B) population density varied between age classes; (C) the prevalence of coccidia increased across the study period in an approximately linear fashion. Points show raw mean prevalence, while lines and shaded areas show predictions and standard errors from model 15 for strongyles and model 3 for coccidia in Table 2.
The prevalence of coccidia was only weakly influenced by year and PD. The model with the lowest AIC value was model 1, which fitted just a linear effect of year (Table 2). The null model, however, had a ΔAIC value of only +0.88, suggesting that the effect of year was only marginal. Also within ΔAIC of 2 were models 3 and 4 (Table 2), both of which included main effects of year. Models that included PD and interactions with age and sex fitted poorly. Overall, there appeared to be a weak increase in the prevalence of coccidia across the study period (Fig. 1C).
For both Nematodirus spp. and M. expansa, the best-supported model was the null model, suggesting no strong influence of either year or PD (Table S1). Although the prevalence of both varied between around 0.05 and 0.20 across the study period, neither linear trends nor non-parametric smooth terms were supported (Fig. S3). As described above, both T. ovis and C. longipes had mean prevalence of <1% and were not analysed (Fig. S3).
Temporal variation in abundance
Strongyle FEC was largely associated with PD rather than year (Table 3). For example, while models including year (e.g. models 1 and 4) generally had lower AIC than the null model, this was not true once PD was added to the model (e.g. compare models 5 and 8). The best-supported model, with AIC −2.55 compared to the next-best model, suggested that strongyle FEC varied with PD in a manner that varied between age and sex classes (model 14). Broadly, males and younger animals had higher FEC than females and adults, and the increase in FEC with increasing PD was greatest in lambs, while in yearlings and adults there was a relatively weak positive association (Fig. 2A). There were also subtle differences among the sexes, which may explain why this model fitted better than the solely age-specific PD model 12 (Table 3). For example, among lambs, the difference between the sexes was greatest at higher PD, while along yearlings, the difference between the sexes was greater at lower PD (Fig. 2A).
Generalized additive mixed-effects model-predicted variation in strongyle FEC and coccidian FOC from the best-supported models in Table 3. (A) Variation in strongyle FEC with population density in each of 6 age and sex classes; (B) variation in coccidian FOC across time in the same age and sex classes. Points show raw mean abundance, while lines and shaded areas show predictions and standard errors from model 14 for strongyles and model 11 for coccidia in Table 3.
In contrast to strongyle FEC, coccidian FOC was chiefly associated with year rather than PD (e.g. compare models 4, 5 and 8 in Table 3). The best-supported model (model 11) had AIC −4.20 relative to the next-best model and supported variation in coccidian FOC across years that varied between age and sex classes (Table 3). The most notable pattern was the increase in FOC among lambs, with the model predicting around a 3-fold increase in FOC across the course of the study (Fig. 2B). Meanwhile, FOC in yearlings and adults appears relatively stable across time. Looking more closely at the model predictions, however, sex differences become apparent, which explain why model 11 was best-supported. For example, the increase in FOC across time appears to have been slightly stronger in female lambs than in male lambs, while in yearlings and adults there appear to have been slight declines across time in males and slight increases in females, with a non-linear pattern apparent in adult females (Fig. S4).
Abundance distribution parameters
Variance in strongyle FEC did not change across time, but it did increase with PD (Table S2), specifically in lambs (Fig. 3A), while the dispersion parameter also increased with density (Table S3; Fig. 3B). Aggregation in FEC varied in a non-linear way with both year and density, although neither was involved in interactions with year or sex (Table S4). Aggregation was predicted to increase from low-to-moderate PD and then plateau (Fig. 3C), while this was the opposite for year (Fig. 3D), although this may have been driven by an influential data point in 1988. Variance in coccidian FOC varied with year but not density as a main effect (Table S5), but interactions suggested that variance did decline with density in lambs but not the other age classes (Fig. 3E), while it increased with year in lambs as well (Fig. 3F). The dispersion parameter for FOC also increased across time (Fig. 3G), but it did not vary with year and was not involved in any interactions (Table S6). Finally, although aggregation did not vary with either year or density as a main effect (Table S7), there was an interaction between year and age, with FEC becoming less aggregated in lambs (higher k) and more aggregated in yearlings (lower k) with time (Fig. 3H).
Variation in distribution parameters for (A–D) strongyle fecal egg count (FEC) and (E–H) coccidian fecal oocyst count (FOC) with year and population density. Points show raw parameter estimates from the data and lines show results of linear models ± 1s.e.; lambs are in orange, yearlings in red, and adults in blue.
Discussion
Long-term studies of wildlife disease are vital to enable us to understand how these systems respond to individual, population and environmental change, and to separate the drivers of host–pathogen/parasite dynamics (Barroso et al., Reference Barroso, Acevedo and Vicente2021). In this study, we describe the nature of temporal trends in parasite infections in the Soay sheep of St. Kilda across 31 years. Our results show that the prevalence and/or abundance of 2 prevalent and pathogenic types of parasite have changed across the study period in age- and sex-dependent ways.
Strongyles and coccidia are the most prevalent and abundant parasites infecting St. Kilda Soay sheep, and here we show that the prevalence of both, and the abundance of coccidia, increased across time in an age and/or sex-dependent manner. These increases were over-and-above any effect of host PD, which was positively associated with prevalence and abundance of most of our parasite taxa, and which itself increased across time by around 5 sheep per year between 1988 and 2018 (Fig. S1; linear regression of density on year estimate 4.95 ± 2.25s.e., F = 4.85, P = 0.036). Furthermore, we accounted for host age and sex in our models, so these year-on-year trends are independent of these factors. Previous studies of wildlife disease have determined that changes in host PD can play a role in generating temporal patterns of disease prevalence, such as the case of tuberculosis in European badgers (Meles meles), red deer (Cervus elaphus), fallow deer (Dama dama) and wild boar (Sus scrofa), which has increased across time in these populations at least partly in response to host density (Delahay et al., Reference Delahay, Walker, Smith, Wilkinson, Clifton-Hadley, Cheeseman, Tomlinson and Chambers2013; Vicente et al., Reference Vicente, Barasona, Acevedo, Ruiz-Fons, Boadella, Diez-Delgado, Beltran-Beck, González-Barrio, Queirós, Montoro, de la Fuente and Gortazar2013; Barroso et al., Reference Barroso, Barasona, Acevedo, Palencia, Carro, Negro, Torres, Gortázar, Soriguer and Vicente2020). Due to the particularly long-term nature of our study, we were able to disentangle effects of density and year in our analyses, and found that both influenced prevalence and/or abundance of strongyles and/or coccidia in an age-dependent manner. While prevalence of coccidia increased with time in a roughly linear fashion, prevalence of strongyles increased linearly in lambs and adults, but followed a non-linear trajectory in yearlings, with an increase in prevalence early in the study followed by stasis thereafter. Meanwhile, although prevalence of strongyles was high in lambs (94% overall), it still increased with host PD in a linear fashion; in yearlings and adults, the pattern was more of an increase between low and moderate density and stasis at higher density. Transmission is expected to be higher at higher density in this population (Hayward et al., Reference Hayward, Wilson, Pilkington, Pemberton and Kruuk2009) and lambs at 4 months are likely to be more susceptible to infection due to their less effective parasite-specific immune responses compared to adults and yearlings (Sparks et al., Reference Sparks, Watt, Sinclair, Pilkington, Pemberton, McNeilly, Nussey and Johnston2019). The linear association between density and prevalence in lambs may reflect the fact that naïve lambs act as ‘sentinels’ of infection pressure (Hayward et al., Reference Hayward, Garnier, Watt, Pilkington, Grenfell, Matthews, Pemberton, Nussey and Graham2014), but once transmission reaches a certain level in older animals, their acquired immunity may prevent further increases in infection. This result was also reflected in that for strongyle abundance: lambs experienced increased strongyle FEC at higher density, while the effect of density was much weaker in adults and yearlings (Fig. 2A). By far the most striking result was the marked increase in coccidian FOC across the study period in lambs that was largely absent from yearlings and adults (Fig. 2B). Coccidian parasites can cause considerable clinical signs and pathology in lambs from domestic sheep flocks. In these livestock systems, infection levels generally peak just after weaning (Chartier and Paraud, Reference Chartier and Paraud2012), which is just before we sampled the St. Kilda Soay Sheep in August. If lambs also act as sentinels of coccidian infection, our results suggest that infection pressure applied by coccidia has increased across the course of the study independently of host density. Explaining the temporal trends we observed, and particularly the striking increase in coccidia abundance in lambs, remains a task for further, detailed, analysis, but several potential mechanisms could be involved.
When explaining long-term temporal trends in ecology, climate and specifically global warming naturally spring to mind (Huntley et al., Reference Huntley, Fürsich, Alberti, Hethke and Liu2014). Since the start of the Soay sheep study in 1985, mean temperature across all parts of the year on St. Kilda has increased by ~0.03°C per year (Crawley et al., Reference Crawley, Pakeman, Albon, Pilkington, Stevenson, Morrissey, Jones, Allan, Bento, Hipperson, Asefa and Pemberton2021). Both strongyles and coccidia are directly transmitted, and both temperature and moisture are associated with the survival and development of strongyle larvae (O'Connor et al., Reference O'Connor, Walkden-Brown and Kahn2006; Rose et al., Reference Rose, Hoar, Kutz and Morgan2014) and coccidian oocysts (Waldenstedt et al., Reference Waldenstedt, Elwinger, Lundén, Thebo and Uggla2001; Makau et al., Reference Makau, Gitau, Muchemi, Thomas, Cook, Wardrop, Fèvre and de Glanville2017) in the external environment. These observations are reflected in the associations between, for example, warmer temperatures and greater parasite abundances in other wild (Holand et al., Reference Holand, Jensen, Kvalnes, Tufto, Pärn, Sæther and Ringsby2019; Mennerat et al., Reference Mennerat, Charmantier, Perret, Hurtrez-Boussès and Lambrechts2019) and domestic animal populations (Kenyon et al., Reference Kenyon, Sargison, Skuce and Jackson2009; van Dijk et al., Reference van Dijk, Sargison, Kenyon and Skuce2009). Indeed, a warming climate has been predicted to cause an increase in infection pressure by gastrointestinal nematodes in domestic sheep (van Dijk et al., Reference van Dijk, David, Baird and Morgan2008; Morgan and van Dijk, Reference Morgan and van Dijk2012; Rose et al., Reference Rose, Wang, van Dijk and Morgan2015; Rose et al., Reference Rose, Caminade, Bolajoko, Phelan, van Dijk, Baylis, Williams and Morgan2016). Projected changes in precipitation are also likely to play a role (O'Connor et al., Reference O'Connor, Walkden-Brown and Kahn2006). The availability of long-term climate data on St. Kilda (Crawley et al., Reference Crawley, Pakeman, Albon, Pilkington, Stevenson, Morrissey, Jones, Allan, Bento, Hipperson, Asefa and Pemberton2021), means that the role of climate in driving temporal trends in parasitology can be thoroughly investigated in future analyses. Simulation models showing how parasite epidemiology may change under various climate change scenarios have been developed for domestic sheep (Rose et al., Reference Rose, Wang, van Dijk and Morgan2015; Rose Vineer et al., Reference Rose Vineer, Verschave, Claerebout, Vercruysse, Shaw, Charlier and Morgan2020), and their adaptation for the Soay sheep population could enable projections of future epidemiology in this population.
Another factor potentially explaining the temporal trends in parasite prevalence and abundance could be changes in selection pressure in the population. Higher strongyle burdens (in particular) are associated with reduced winter survival in Soay sheep, an association that is restricted to lambs and absent in adults (Gulland, Reference Gulland1992; Hayward et al., Reference Hayward, Wilson, Pilkington, Clutton-Brock, Pemberton and Kruuk2011). If selection has weakened such that lambs with less effective immune responses and higher parasite burdens are now more likely to survive, increases in parasite abundance may become apparent. This hypothesis is, however, somewhat troubled by the observation that, despite the much lower strongyle-specific antibody levels in lambs compared to adults (Sparks et al., Reference Sparks, Watt, Sinclair, Pilkington, Pemberton, McNeilly, Nussey and Johnston2019), strongyle-specific antibodies are associated with survival in adult females, but not in lambs (Sparks et al., Reference Sparks, Watt, Sinclair, Pilkington, Pemberton, Johnston, McNeilly and Nussey2018). Selection on August FEC in lambs does vary among years, with FEC being associated with survival in low-mortality years but not high-mortality years (Hayward et al., Reference Hayward, Wilson, Pilkington, Clutton-Brock, Pemberton and Kruuk2011). Temporal changes in patterns of selection have been observed in other studies of wild populations: warming temperature can explain a shift from selection for larger to selection for smaller forehead patch size in male collared flycatchers (Ficedula albicollis) (Evans and Gustafsson, Reference Evans and Gustafsson2017) and weakening of selection on hatch date in pied flycatchers (Ficedula hypoleuca) (Visser et al., Reference Visser, Gienapp, Husby, Morrisey, de la Hera, Pulido and Both2015). Despite all of these findings, no formal analysis of temporal or environment-dependent variation in selection on immunoparasitological traits has been conducted on the Soay sheep population.
We saw limited evidence for changes in the distribution of strongyles across time, which was perhaps unsurprising given that we found no change in mean abundance of strongyles across time. We did, however, see variation in all 3 parameters with density (Tables S2–S4). Among lambs, variance increased with density as did the mean, but so does the dispersion index I, suggesting that the variance to mean ratio also increased (Fig. 3B). Meanwhile, at low PD, FEC was more aggregated (lower k). We observed variation in all 3 distributional parameters for coccidia (Fig. 3), and these were more closely linked to year than to density (Tables S5–S7). An intriguing result from the point of view of parasite epidemiology is the age-specific trend in the aggregation parameter k for coccidian FOC (Fig. 3F), which suggests that FOC is becoming less aggregated among lambs (increasing k) and more aggregated in adults and especially yearlings (decreasing k), although the trends in lambs and adults are relatively weak. This result broadly supports theory that predicts that aggregation should decrease as encounters with parasites (i.e. the force of infection) increase (Gourbière et al., Reference Gourbière, Morand and Waxman2015). Given that k increases (and aggregation decreases) with the mean, it is perhaps unsurprising that we saw this result in lambs, and visual inspection of the distribution of FOC in lambs suggests a year-on-year reduction in aggregation (Fig. S5). Where parasites are abundant, increasing k places a greater proportion of hosts under potentially health-limiting parasite burdens (Poulin and Vickery, Reference Poulin and Vickery1993), suggesting that lambs are increasingly facing selection pressures from coccidian infection. A major advantage of the Soay sheep system is that hypotheses about selection on parasite resistance can be tested, thanks to the detailed data on survival, reproduction and parasitology that have been collected across the decades (Clutton-Brock and Pemberton, Reference Clutton-Brock, Pemberton, Clutton Brock and Pemberton2004), plus a cache of stored blood samples that has been interrogated for parasite-specific antibody responses (Hayward et al., Reference Hayward, Garnier, Watt, Pilkington, Grenfell, Matthews, Pemberton, Nussey and Graham2014; Froy et al., Reference Froy, Sparks, Watt, Sinclair, Bach, Pilkington, Pemberton, McNeilly and Nussey2019; Sparks et al., Reference Sparks, Watt, Sinclair, Pilkington, Pemberton, McNeilly, Nussey and Johnston2019).
A potential caveat of this study is our use of FEC/FOC, rather than actual parasite numbers, and their relatively low resolution of 100 eggs/oocysts per gram. The advantages of our method are that it requires basic equipment and is easy to perform in the field; it is non-invasive and can be measured repeatedly on the same animals without the need for culling; it accounts for the effects of host immunity on both worm number and fecundity; and, in the case of Soay sheep on St. Kilda and elsewhere, it is linearly associated with actual worm burden (Grenfell et al., Reference Grenfell, Wilson, Isham, Boyd and Dietz1995; Boyd, Reference Boyd1999; Wilson et al., Reference Wilson, Grenfell, Pilkington, Boyd, Gulland, Clutton-Brock and Pemberton2004). As such, while this measure may not be a true reflection of worm burden, it has many desirable properties that enable us to estimate the salient drivers of variation in parasite infection in the population. We also only analysed data collected in August, when animals are captured for sampling and measurement, despite considerable amounts of data being available across the year. August, however, is the main source of data and the only month where data have been collected in every year; in other months, data have been collected sporadically for various specific research projects in a limited number of years. While marked seasonal variation in parasite infection exists in this population, with variation in this seasonality associated with age, sex and reproductive status (Sweeny et al., Reference Sweeny, Corripio-Miyar, Bal, Hayward, Pilkington, McNeilly, Nussey and Kenyon2022), focusing on August data enabled us to provide the best representation of changes over time.
The strength of our study lies in its long-term nature: data have been collected at the same time in every year since 1988, making it a rarity among studies of wildlife disease (Barroso et al., Reference Barroso, Acevedo and Vicente2021). The longitudinal nature of the study has enabled us to identify temporal trends and to separate them from effects of host density and demography. The next challenge is to make use of the associated meta-data in order to identify the forces that are driving the changes we have observed. On the other hand, the main weakness of our study lies in the lack of parasite species-specificity in our data, which arises because of the difficulty of distinguishing individual parasite species based on egg/oocyst morphology alone. The focus has been on fecal sampling throughout the St. Kilda Soay sheep study largely because of its non-invasive nature. Previous work on the Soay sheep has identified helminth species based on adult worm morphology by necropsy of animals that died over winter (Craig et al., Reference Craig, Pilkington and Pemberton2006) and Eimeria species through labour-intensive sporulation of oocysts and identification based on morphology (Craig et al., Reference Craig, Pilkington, Kruuk and Pemberton2007). These studies have provided intriguing results: for example, while T. axei and T. vitrinus dominate the strongyle species in lambs, T. circumcincta is much more common in adults (Craig et al., Reference Craig, Pilkington and Pemberton2006); meanwhile, Eimeria species diversity is high in lambs but declines with age (Craig et al., Reference Craig, Pilkington, Kruuk and Pemberton2007). It is currently impossible to determine which species are responsible for the temporal trends observed here, but elucidating this is important, because of the different locales of the parasites within the host, and their potentially different impacts on host health. Application of new meta-barcoding techniques to nematodes (Avramenko et al., Reference Avramenko, Redman, Lewis, Yazwinski, Wasmuth and Gilleard2015) that have already been applied in sheep (Redman et al., Reference Redman, Queiroz, Bartley, Levy, Avramenko and Gilleard2019), and those developed for coccidia (Heitlinger et al., Reference Heitlinger, Ferreira, Thierer, Hofer and East2017), which are in the process of being developed for sheep, could help in this regard.
Overall, our results provide evidence for temporal trends in the prevalence and abundance of fitness-limiting parasites across 3 decades in a wild ungulate population. These trends were generally stronger in the youngest and most susceptible animals in the population than in older animals, and were absent for rarer parasite species. Decades of intensive, individual-based study of this population, plus new data on parasite species identity, will enable the role of environmental change and selection pressures in driving these dynamics to be determined, potentially providing unrivalled insight into the drivers of temporal variation in parasite dynamics in the wild.
Supplementary material
The supplementary material for this article can be found at https://doi.org/10.1017/S0031182022001263.
Data availability
Data used in this publication can be accessed via Mendeley Data at https://data.mendeley.com/.
Acknowledgements
We thank Dan Nussey for insightful comments on earlier drafts of the manuscript. We acknowledge the National Trust for Scotland for permission to work on St. Kilda and QinetiQ, Eurest and Kilda Cruises for logistics and support. We also thank I. Stevenson and many volunteers who have collected field data and samples and all those who have contributed to keeping the project going.
Author's contributions
A. D. H. and J. M. B. conceived the study; J. M. P. and J. G. P. led coordination of fieldwork and data collection on St. Kilda; A. D. H. analysed the data with input and suggestions from D. Z. C., A. F., F. K., R. J. P. and A. B. P.; A. D. H. led writing of the manuscript with edits and suggestions from J. M. B., D. Z. C., Y. C. M., A. F., M. D. F., F. K., T. N. M., R. J. P., A. B. P., J. M. P., A. R. S. and K. W.
Financial support
Field data collection has been supported principally by NERC over many years, with some funding from the Wellcome Trust. A. D. H. is supported by a Moredun Foundation Research Fellowship.
Conflict of interest
None.
Ethical standards
All sampling was carried out in accordance with UK Home Office regulations under Project License PP4825594.
Open access
