Long-term temporal trends in gastrointestinal parasite infection in wild Soay sheep

Abstract

Monitoring the prevalence and abundance of parasites over time is important for addressing their potential impact on host life histories, immunological profiles and their influence as a selective force. Only long-term ecological studies have the potential to shed light on both the temporal trends in infection prevalence and abundance and the drivers of such trends, because of their ability to dissect drivers that may be confounded over shorter time scales. Despite this, only a relatively small number of such studies exist. Here, we analysed changes in the prevalence and abundance of gastrointestinal parasites in the wild Soay sheep population of St. Kilda across 31 years. The host population density (PD) has increased across the study, and PD is known to increase parasite transmission, but we found that PD and year explained temporal variation in parasite prevalence and abundance independently. Prevalence of both strongyle nematodes and coccidian microparasites increased during the study, and this effect varied between lambs, yearlings and adults. Meanwhile, abundance of strongyles was more strongly linked to host PD than to temporal (yearly) dynamics, while abundance of coccidia showed a strong temporal trend without any influence of PD. Strikingly, coccidian abundance increased 3-fold across the course of the study in lambs, while increases in yearlings and adults were negligible. Our decades-long, intensive, individual-based study will enable the role of environmental change and selection pressures in driving these dynamics to be determined, potentially providing unparalleled insight into the drivers of temporal variation in parasite dynamics in the wild.

Introduction

Parasite infections are widely known to follow seasonal patterns, with variation in prevalence and abundance of parasites apparent within years (Altizer et al., 2006). 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., 1988; Montgomery and Montgomery, 1988; Haukisalmi and Henttonen, 1990; Halvorsen et al., 1999; Pelletier et al., 2005; Rossin et al., 2010; Turner and Getz, 2010; Cizauskas et al., 2015). 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., 2004; Armitage, 2012; Andreassen et al., 2021; Crawley et al., 2021), and trends in phenological traits such as timing of breeding (Charmantier et al., 2008; Bonnet et al., 2019), morphological traits (Ozgul et al., 2009) and sexually selected traits (Evans and Gustafsson, 2017). 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., 1998; Babayan et al., 2018). Such long-term changes may be linked to a number of drivers, including climatic variation (Haukisalmi and Henttonen, 1990; Huntley et al., 2014), since the development of transmission stages and vectors may be temperature- and moisture-dependent (Nordenfors et al., 1999; O'Connor et al., 2006; Rose et al., 2014); 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., 2021).

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., 1988; Bajer et al., 2004; Behnke et al., 2008; Knowles et al., 2013; Behnke et al., 2018). Where temporal variation does occur, it appears to consist of fluctuations between years rather than consistent trends across time (Behnke et al., 1999; Haukisalmi and Henttonen, 2000; Irvine et al., 2000; Lachish et al., 2011; Grzybek et al., 2015; Peet et al., 2022). In some cases, however, both colonization of hosts by novel species and extinction of previously prevalent species have been reported (Kennedy et al., 2001; Lyndon and Kennedy, 2001; Grzybek et al., 2015). 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., 2021). 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., 2021). 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., 2021) or were non-destructive (Sweeny et al., 2021), 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., 1995; Haukisalmi and Henttonen, 1999). 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, 1993). 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, 1986; Woolhouse et al., 1997; Wilson et al., 2002). 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, 1978; May and Anderson, 1978). 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, 1982). 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, 2004b). The sheep population has experienced significant temporal change during this time, with average body size declining (Ozgul et al., 2009) and population size increasing along with grassland productivity (Crawley et al., 2021). The known impacts of gastrointestinal parasites on condition and fitness in the Soay sheep (Gulland, 1992; Craig et al., 2008; Hayward et al., 2011) 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, 2004); 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, 2004). 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., 2004). 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, 2005). 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, 1992; Craig et al., 2009). 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., 2009). 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.

Table 1. 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, 2020). 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., 2009). 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., 2004); 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.

Table 2. 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 σ ² =  μ +  μ ²/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.

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 = σ ²/μ) and the negative binomial aggregation parameter (k =  μ ²/σ ² −  μ), where x_i, μ 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.

Fig. 1. 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).

Fig. 2. 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).

Fig. 3. 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., 2021). 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., 2013; Vicente et al., 2013; Barroso et al., 2020). 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., 2009) 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., 2019). 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., 2014), 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, 2012), 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., 2014). 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., 2021). 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., 2006; Rose et al., 2014) and coccidian oocysts (Waldenstedt et al., 2001; Makau et al., 2017) 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., 2019; Mennerat et al., 2019) and domestic animal populations (Kenyon et al., 2009; van Dijk et al., 2009). Indeed, a warming climate has been predicted to cause an increase in infection pressure by gastrointestinal nematodes in domestic sheep (van Dijk et al., 2008; Morgan and van Dijk, 2012; Rose et al., 2015; Rose et al., 2016). Projected changes in precipitation are also likely to play a role (O'Connor et al., 2006). The availability of long-term climate data on St. Kilda (Crawley et al., 2021), 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., 2015; Rose Vineer et al., 2020), 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, 1992; Hayward et al., 2011). 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., 2019), strongyle-specific antibodies are associated with survival in adult females, but not in lambs (Sparks et al., 2018). 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., 2011). 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, 2017) and weakening of selection on hatch date in pied flycatchers (Ficedula hypoleuca) (Visser et al., 2015). 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., 2015). 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, 1993), 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, 2004), plus a cache of stored blood samples that has been interrogated for parasite-specific antibody responses (Hayward et al., 2014; Froy et al., 2019; Sparks et al., 2019).

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., 1995; Boyd, 1999; Wilson et al., 2004). 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., 2022), 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., 2021). 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., 2006) and Eimeria species through labour-intensive sporulation of oocysts and identification based on morphology (Craig et al., 2007). 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., 2006); meanwhile, Eimeria species diversity is high in lambs but declines with age (Craig et al., 2007). 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., 2015) that have already been applied in sheep (Redman et al., 2019), and those developed for coccidia (Heitlinger et al., 2017), 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.