Winter Tick Burdens for Moose Are Positively Associated With Warmer Summers and Higher Predation Rates

Climate change is expected to modify host-parasite interactions which is concerning because parasites are involved in most food-web links, and parasites have important influences on the structure, productivity and stability of communities and ecosystems. However, the impact of climate change on host–parasite interactions and any cascading effects on other ecosystem processes has received relatively little empirical attention. We assessed host-parasite dynamics for moose (Alces alces) and winter ticks (Dermacentor albipictus) in Isle Royale National Park over a 19-year period. Specifically, we monitored annual tick burdens for moose (estimated from hair loss) and assessed how it covaried with several aspects of seasonal climate, and non-climatic factors, such as moose density, predation on hosts by wolves (Canis lupus) and wolf abundance. Summer temperatures explained half the interannual variance in tick burden with tick burden being greater following hotter summers, presumably because warmer temperatures accelerate the development of tick eggs and increase egg survival. That finding is consistent with the general expectation that warmer temperatures may promote higher parasite burdens. However, summer temperatures are warming less rapidly than other seasons across most regions of North America. Therefore, tick burdens seem to be primarily associated with an aspect of climate that is currently exhibiting a lower rate of change. Tick burdens were also positively correlated with predation rate, which could be due to moose exhibiting risk-sensitive habitat selection (in years when predation risk is high) in such a manner as to increases the encounter rate with questing tick larvae in autumn. However, that positive correlation could also arise if high parasite burdens make moose more vulnerable to predators or because of some other density-dependent process (given that predation rate and moose density are highly correlated). Overall, these results provide valuable insights about interrelationships among climate, parasites, host/prey, and predators.


INTRODUCTION
Climate change influences organisms directly (i.e., by influencing physiology and behavior) and indirectly by modifying trophic interactions, such as consumer-resource relationships (Ottersen et al., 2001;Stenseth et al., 2002;Walther et al., 2002). Parasitism is the most common consumer strategy (Poulin and Morand, 2000;De Meeûs and Renaud, 2002), and parasites play an important role in shaping the composition, structure, productivity and stability of communities and ecosystems (Minchella and Scott, 1991;Poulin, 2002a,b, 2005;Wood and Johnson, 2015). There are growing concerns about how climate change will modify host-parasite interactions. For example, new host-parasite relationships seem likely to form as species shift their geographic ranges in response to climate change and some host species will be exposed to parasites to which they have no immunity or coevolutionary history (Brooks and Hoberg, 2007). Additionally, outbreaks of parasites have already been linked to extreme weather (Hudson et al., 2006;Wegner et al., 2008) and extreme weather events are becoming more frequent and intense as the climate warms (National Academies of Science [NAS], 2016). Moreover, warming temperatures could favor earlier emergence dates, faster development, increased survival, and extended periods of activity for many parasite species (Ogden et al., 2006;Poulin, 2006;Calero-Torralbo et al., 2013). For those reasons, climate change is generally expected to favor parasites over hosts with more severe effects of parasites predicted in the future (Patz et al., 2003;Pounds et al., 2006;Barber et al., 2016). Any such impacts of climate change on host-parasite relationships are likely to be substantial and widespread because parasites account for a large portion of the world's biodiversity (Dobson et al., 2008;Kuris et al., 2008) and are involved in the majority of food-web links (Lafferty et al., 2006;Dunne et al., 2013).
One challenge of predicting how climate change may modify host-parasite interactions is that the rate of climate change varies substantially among seasons. For example, temperatures are warming more rapidly in winter than in summer across most regions of North America (Vose et al., 2017). Additionally, winter and spring precipitation are expected to increase by up to 20% across large parts of midwestern United States, due to more intense heavy rain or snowfall events, but no significant changes in precipitation are expected during the summer and autumn (Hayhoe et al., 2010). These seasonal differences are significant because parasite-host relationships may be more sensitive to weather conditions during some seasons, especially for parasites that live freely during a portion of their life-cycle (Harvell et al., 2002;Cizauskas et al., 2017).
The impact of climate change on host-parasite interactions and their ecosystem consequences have received limited empirical attention Møller et al., 2013). Moreover, few studies have adequately considered the influence of climate in the presence of other dynamic factors, such as anthropogenic changes in land-use or exploitation of hosts, some of which are highly correlated with climate change (Rohr et al., 2011). This shortcoming is presumably due to the difficulty of simultaneously monitoring parasites, hosts, weather, and other potentially confounding non-climatic influences over sufficiently long periods of time and large spatial scales.
Here we assess the dynamics of an ectoparasite of North American cervids -the winter tick (Dermacentor albipictus)over a 19-year period (2001-2019). Specifically, we monitored interannual variation in tick burden for a population of moose (Alces alces) in Isle Royale National Park, and assessed the extent that tick burden covaried with seasonal temperatures and precipitation, moose (host) density, the rate of predation by wolves (Canis lupus) on the host and with wolf abundance. We primarily focused on weather conditions during April through September, which is the time of year when winter ticks are free-living (Drew and Samuel, 1987;Samuel, 2004). We focused on moose density and wolf predation because they are two of the dominant non-climatic factors in this portion of the food web (details below). Moreover, host density is believed to have an important influence on parasite abundance and transmission rates (Mugabo et al., 2015;Stringer and Linklater, 2015;Wang et al., 2016). There is also a growing appreciation that predators may influence parasite abundance (Samková et al., 2019). For example, predators may trigger changes in host behavior or physiology which can affect the prevalence and intensity of parasite infections (Zukowski et al., 2020). Predators also disproportionally prey on hosts infected with parasites (Hudson et al., 1992;Møller and Nielsen, 2007), which could reduce parasite burdens the following year. Furthermore, predation may affect parasite burdens via its effect on host population dynamics. The specific hypotheses that we test are outlined in Tables 1, 2.

Study System
The life cycle of winter ticks begins with females laying eggs near the soil surface in June (Drew and Samuel, 1987;Samuel, 2004). Each female can produce several thousand eggs which hatch between August and September. For winter ticks, warmer summer temperatures are associated faster development of eggs into larvae and increased egg survival (Drew and Samuel, 1987). For other tick species, warmer summer temperatures during the oviposition period tend to reduce the time it takes a female to form and lay eggs, reduce oviposition failure, and increase the total number of eggs produced (Ogden et al., 2004;Lysyk, 2014). Those relationships with temperature plausibly apply to winter ticks as well.
After hatching, clumps of larvae ascend vegetation and quest for ungulate hosts typically sometime between September and October. Questing activity ends when temperatures drop below freezing for an extended period of time. As such, warmer temperatures in autumn may extend the duration of the questing period, which may increase the probability of larvae finding a host. Larvae cannot survive the winter without a host.
Upon finding a host, larvae feed and then molt into nymphs between October and November. Nymphs remain dormant until 1 | Hypotheses pertaining to the influence of non-climatic factors on tick burdens for moose as indicated by the extent of hair loss (see Figure 1).

Variable
Description Expected relationship

Moose
Moose density (moose per km 2 , estimated in mid-January to early February -a few months before hair loss was estimated) Positive/Negative A positive relationship could arise if high moose densities increase the likelihood of questing larvae finding a host. Alternatively, the risk of being parasitized may decline when more hosts are available, due to an "encounter-dilution effect" (Mugabo et al., 2015) Moose lag Moose density, estimated in previous winter Positive Previous studies observed tick abundance to be positively associated with moose density following a 1-year lag (Samuel, 2004(Samuel, , 2007. Greater moose densities may equate to greater availability of hosts upon which ticks can develop leading to higher parasite burdens the following year Predation Predation rate on moose by wolves (estimated in between mid-January and early March -a few months before hair loss was estimated) Positive/Negative Increasing predation risk can trigger changes in ungulate behavior, such as movement and habitat selection (Montgomery et al., 2013;Kohl et al., 2018); habitat selection during the autumn may influence encounter rates between hosts and questing tick larvae (Healy et al., 2018;Blouin et al., 2021) Predation lag Predation rate on moose by wolves, estimated in previous year Negative Predators disproportionally prey on hosts infected with parasites (Hudson et al., 1992;Møller and Nielsen, 2007), and high predation rates on such individuals may reduce the number of female ticks surviving to produce eggs the following spring. Additionally, predation by wolves tends to reduce moose population growth rates (Vucetich et al., 2011), which may subsequently influence tick abundance the following year by reducing the number of hosts available Moose density is an indicator of host availability. Predation rate is the proportion of the moose population killed by wolves. Predation rate and wolf abundance are both indicators of temporal variation in predation risk for moose at the population level. However, moose population growth rates are more closely associated with predation rate than predator densities (Vucetich et al., 2011). Note that all the non-climatic predictor variables are highly correlated with one another (Supplementary Table 1).
late January and then feed on their hosts and molt into adults. Adult ticks of both sexes take their final blood meal and mate on their host between March and April. Adults then detach from their hosts in late April. If ticks "drop-off " hosts onto snow they are unlikely to survive (Drew and Samuel, 1987;Wilton and Garner, 1993). The timing of the winter tick's life-cycle (e.g., questing, feeding on hosts, egg laying, detaching from hosts) is thought to be regulated by photoperiod (Wright, 1971;Drew and Samuel, 1985). The main hosts for winter ticks are North American ungulates, especially Cervidae (Samuel, 2004). Winter ticks impact moose more severely than other host species, presumably because moose only arrived to North America about 10,000-24,000 years ago and have less co-evolutionary history with winter ticks than other North American host species (Samuel, 2004). Whereas other cervid hosts appear to be effective in limiting the number of larvae acquired by grooming in the fall, moose don't start intensively grooming until mid-winter when nymphs end their dormancy (see Figure 6.4 in Samuel, 2004). As a result, individual moose routinely acquire tens of thousands ticks during a single winter (Samuel, 2007).
The irritation caused by tick bites causes moose to groom intensely enough to result in significant damage to (and loss of) their winter coat, which increases energetic demands for thermoregulation (Glines and Samuel, 1989;Samuel, 2004). Blood consumption by ticks often leads to chronic anemia, protein deficits, and substantial energetic costs for moose (Glines and Samuel, 1989;Musante et al., 2007;Wünschmann et al., 2015). These impacts often result in lower fecundity and an elevated risk of death, especially for calves (Musante et al., 2007;Samuel, 2007;Jones et al., 2019;Pekins, 2020). Winter ticks are cited as being an important cause of recent declines in moose populations in the north-eastern United States (Ellingwood et al., 2020). As the climate warms, milder temperatures in fall and early winter may extend the questing period for winter tick larvae, and are therefore expected to result in higher tick burdens for moose (Jones et al., 2017(Jones et al., , 2019Ellingwood et al., 2020). However, there is an apparent lack of long-term studies assessing how moose-tick interactions are influenced by climate warming.
Isle Royale National Park (IRNP) is an archipelago in Lake Superior, North America (47 • 50 N, 89 • 00 W), comprised of a large island (544 km 2 ) and dozens of smaller islets (most of which are <2 km 2 ). Isle Royale is also known as Minong by local indigenous communities and is under the stewardship of the Grand Portage Anishinaabe and U.S. National Park Service. The climate is temperate, characterized by short summers (July-August) and long winters with snow cover typically starting at the end of October and lasting until April (Supplementary Material 2 for details). IRNP has a lower diversity of mammals than the mainland (Peterson, 1977). Moose are the only large herbivore and gray wolves (Canis lupus) are their only predator. Moose population dynamics are thought to be strongly influenced by wolf predation (Vucetich et al., 2011).
An advantage of studying moose-tick interactions in this system is the absence of alternative host species or human TABLE 2 | Hypotheses pertaining to the influence of weather variables on tick burdens for moose as indicated by the extent of hair loss (see Figure 1).

Variable name Variable description Expected relationship
Underlying mechanism

Snow
Average snow depth (mid-January to early March) Negative Winters with deep snow increase the likelihood of snow cover in late April when ticks detach from hosts, which may reduce adult tick survival (Drew and Samuel, 1987) Temp APR Mean daily maximum temperatures, April Positive Similar to the mechanism described for snow, warmer April temperatures may reduce the likelihood of snow cover when ticks detach from hosts (Wilton and Garner, 1993) Spring Day of the year that daily minimum temperatures first remain above 0 • C for 7 consecutive days Negative Similar to the mechanisms described above, a later start to spring may reduce the survival of adult female ticks once they have detached from their hosts Temp JUL Mean daily maximum temperatures, July Positive Warmer temperatures in summer may promote faster egg development and increased egg survival (Drew and Samuel, 1987) Precipitation Total precipitation, July Positive Wetter summers may reduce the risk of tick eggs and larvae becoming desiccated and dying (Knülle, 1966;Yoder et al., 2015) Temp SEPT Mean daily maximum temperatures, September Positive Warmer temperatures in September may increase the activity of questing larvae (Drew and Samuel, 1985) Autumn Day of the year that daily maximum temperatures first remain below 0 • C for 7 consecutive days Positive A later end to autumn (hence a later start to winter) may extend the questing period for larvae and thereby increase the number of larvae finding hosts (Drew and Samuel, 1985) Winter The number of days between the start of winter and the start of spring Negative Shorter winters are expected to benefit ticks by either extending the duration of the autumnal questing period or by increasing the likelihood that detached female ticks will survive to produce eggs (Jones et al., 2019) Note that all weather variables refer to year t-1 as we expect a 1-year lag in their effect on hair loss (e.g., snow depth in year t-1 will influence hair loss in year t). For context, the mean daily maximum temperatures are highly correlated with the daily minimum and overall mean daily temperatures.
influences on the habitat or demography of the host and its parasite or predator populations. More precisely, there are no other cervid species on Isle Royale to serve as alternative hosts for winter ticks. Additionally, neither the moose population, wolf population nor forest have been harvested for over a century and no major forest fires occurred during the study.

Estimating Hair Loss
Surveys of hair loss in spring have been found to provide a good indication of interannual variation in the number of winter ticks on moose (Samuel, 2007), and the number of ticks is indicative of the impacts of ticks on moose (Musante et al., 2007). Therefore, to assess interannual variation in tick burden, we estimated the proportion of individuals' winter coat with lost or damaged hair (hereafter, hair loss) in spring every year between 2001 and 2019. Although there is anecdotal evidence of winter ticks impacting moose in this population prior to 2001, we excluded those early observations from this assessment because they were not recorded in a systematic or consistent way. We estimated hair loss from images of the side profiles of moose taken between the start of May until early June. We started collecting images in early May because adult ticks have detached from their hosts by that time and because moose grooming declines dramatically after ticks detach (see Figure 6.4 in Samuel, 2004). We did not collect images after early June because hair loss is obscured by the growth of summer pelage, which begins during the third week of June. Dates for collecting images did not change substantively among years. Because of this consistency, we assumed that estimates of interannual variability in hair loss were not confounded by the timing of sampling.
To collect hair loss images, researchers visited locations known to be regularly used by moose (e.g., inland lakes, mineral licks). Images of individuals were also collected opportunistically whilst researchers conducted other fieldwork and by remote cameras. We collected hair loss images of any individual moose for which we could get a clear view of its side profile, irrespective of the individual's age, sex or severity of hair loss. Prior to 2008, some images were taken by drawing patterns of hair loss on datasheets with a blank profile, whereas after 2008 all images were taken using digital cameras.
For each image, we calculated the proportion of an individual's torso and neck where its winter coat was: undamaged, damaged, obscured from view (for example obscured by vegetation), and missing (i.e., bare skin, see Figure 1 and Supplementary  Figure 1). We only evaluated hair loss on the torso and neck because legs were often obscured by vegetation or by water. We excluded the head because hair loss is typically only observed on the back of the ears which are challenging to completely photograph. If more than 10% of the neck or torso were obscured from view then we excluded the image from the dataset. The response variable (hair loss) was estimated as the proportion of the profile where hair was damaged or lost, excluding small areas obscured from view.
The risk of including repeat samples of the same individual moose within a given year was low because most moose are distinguishable on the basis of body size (adult or yearling), sex, the size and shape of antlers and the dewlap, size and distribution of distinguishing marks and injuries (fibropapillomas, limps, skin wounds, torn ears), and pattern of hair loss. The likelihood of repeatedly sampling the same individuals in multiple years is not known because some of these distinguishing features (e.g., shape FIGURE 1 | Example profile image of a moose used to estimate hair loss on the animal's neck and torso due to winter ticks. Hair loss occurs when moose groom themselves to get rid of winter ticks and was estimated from images of individual moose. In the upper image, undamaged hair is brown, damaged guard hairs are light gray/white in color, and areas where the hair is completely missing (i.e., bare skin) are dark brown or black in color (see also Supplementary Figure 4). The lower image is identical to the one above, except that the three hair categories are indicated with block colors (red is undamaged hair, green is damaged hair, and blue is missing hair). In all statistical analyses, the response variable was the logit transformation of the proportion of hair lost or damaged. That proportion was calculated as the sum of the green and blue areas divided by the sum of the red, green and blue areas.

Statistical Analyses
For all analyses, we logit-transformed the response variable, the proportion of hair lost or damaged (hairloss) because it is bounded between 0 and 1. All analyses were performed in Program-R version 4.0.5 (R Core Team, 2016).
Tick burden may vary between sexes and age-classes of moose (Bergeron, 2011) due to age-class and sex-specific differences in movement and habitat use (Ruckstuhl and Neuhaus, 2002;Montgomery et al., 2013) which could affect encounter rates with questing larvae. To assess this possibility, we built linear mixed-effect models where hairloss was the response variable and age-class (yearling or adult) and sex were included as fixed effects and year included as a random effect, using lme4 package (Bates et al., 2015). We also evaluated a model which included an interaction between age-class and sex as well as the main effects of age-class and sex. This analysis excluded 56 individuals whose sex and age were not recorded.
Next, we assessed the extent to which temporal variation in mean annual hair loss was associated with eight climate variables and six non-climatic variables. The non-climatic variables that we considered were: moose density during year t (moose) and t−1 (moose lag ); predation rate during year t (predation) and t−1 (predation lag ); and wolf abundance during year t (wolf ) and t−1 (wolf lag ). A detailed description of the rational for each of these six non-climatic variables is outlined in Table 1 and in Supplementary Material 2. We used estimates of predation rate, wolf abundance and moose density collected as part of previous studies (Hoy et al., 2019). Predation rate is both a cause specific mortality rate for moose and a useful indicator of predation risk for moose at the population level (Hoy et al., 2019). Predation rate was estimated as Predation = KR × P/N, where KR (kill rate) is an estimate of the number of moose killed per predator, per time unit, N represents prey (moose) abundance and P represents predator (wolf) abundance using the methods described in Vucetich et al. (2011). Moose abundance, wolf abundance and kill rate were estimated annually from aerial surveys conducted between late January and February each year throughout the study period (Gasaway et al., 1986;Peterson and Page, 1988).
The climatic variables that we considered were: the average snow depth in winter (snow), total precipitation in July (precipitation) and the mean daily maximum temperatures in April (temp APR ), July (temp JUL ), and September (temp SEPT ). For context, April, July, and September correspond to critical points in the life-cycle of the winter tick. April is when mated female ticks detach from their hosts, July is when tick eggs are on the ground and developing into larvae, and September is when tick larvae are ascending vegetation and questing for hosts (Samuel, 2004). Additionally, we included the first day of the year that daily maximum temperatures remained below freezing (0 • C) for seven consecutive days (autumn) to represent the end autumnal questing period (i.e., the start of winter). Moreover, we included the first day of the year that daily minimum temperatures remained above freezing (0 • C) for seven consecutive days (spring) to represent the start of spring. Lastly, we estimated the number of days between the end of autumn and the start of spring to represent the length of winter (winter).
We obtained daily and mean monthly temperatures and precipitation from a weather station in Grand Portage, Minnesota, located approximately 40-60 km from Isle Royale (Western Regional Climate Center [WRCC], 2016). This is the closest weather station to IRNP that has continuous records of weather data for the duration of the study period. We measured snow depth on a daily basis between mid-January and early March each year at a single site at the west end of Isle Royale, near the location of the basecamp used to conduct aerial survey counts of moose abundance in winter (Hoy et al., 2019(Hoy et al., , 2020b.
We then averaged those daily measurements to estimate the mean annual snow depth. Note, that we built models where hair loss in calendar year t was a function of climatic variables in calendar year t−1. The rational for including each of these eight weather variables is outlined in Table 2. The extent that each of the climatic and non-climatic variables varied over time and were correlated with one another is detailed in Supplementary  Figures 4, 5 and Supplementary Table 1.
To evaluate which candidate variables were useful predictors of hair loss we used the dredge function of the MuMIn package in Program-R (Bartoń, 2018). The dredge function assessed models with all possible combinations of variables included in a global model and ranked models on the basis of Akaike's Information Criterion corrected for small sample size (AICc). The global model contained all climatic and non-climatic predictor variables (moose, moose lag , predation, predation lag , wolf, wolf lag , precipitation, snow, temp APR , temp JUL , temp SEPT , autumn, spring, winter). Because this analysis is focused on understanding factors leading to temporal variation in tick-burdens, we first ran a population-level analysis where the response variable was the mean annual proportion of hair loss (i.e., hair loss averaged for all individuals sampled in a given year and then logit transformed) using linear models. We also repeated this analysis where the response variable was hair loss of an individual using linear mixed-effect models with year included as a random effect (hereafter, individual-level analysis). We limited the maximum number of variables included in any candidate model to four because the duration of the timeseries is 19 years. Because we did not have a priori reason to think that any particular two-way interactions would be significant, we did not include any two-way interactions in the global model. Instead, we built ad hoc models evaluating all two-way interactions involving main effects that were included in the most parsimonious model identified by the dredge function.
We report the best model identified by the dredge function and all models with AICc < 2. For additional context, we also report the null (intercept only) model. Importantly, some of the predictor variables were highly correlated with one another (see Supplementary Table 1). However, strong a priori reasoning supports the assessment of hypotheses linked to each variable. Rather than pre-judge which of the correlated variables are likely to be the most important, we let model performance statistics (R 2 and AIC) guide our inferences and then checked for multicollinearity in the best model, as identified by the model selection procedure. After identifying the most parsimonious model, we visually inspected plots of model residuals to check assumptions of homoscedasticity and normality and whether any observations had high leverage. We also formally tested for homoscedasticity (Breusch-Pagan test) and deviations from normality (Shapiro-Wilk and Kolmogorov-Smirnov tests).

Demographic Classes
Mixed-effect models indicated that hairloss did not differ between sexes (Table 3). Although, there was some evidence to suggest that hairloss varied between age-classes, the difference was small, being on average 5% lower for yearlings than for adults (see Supplementary Figure 6). Moreover, age-class explained less than 1% of variance in hair loss among individuals. One reason for age-class appearing to be an important predictor (on the basis of AIC), yet also explaining only a small proportion of the variance in hairloss is that the dataset is composed primarily of adults (85% were adults) which is reflective of the population's age structure (Hoy et al., 2020a). There was no evidence that an interaction between age-class and sex was statistically significant ( Table 3). For these reasons, we ignored differences in hairloss between the age-classes for the subsequent analyses.

Predictors of Hair Loss
The best predictors of temporal variation in hairloss were temp JUL and predation, with hairloss tending to be greater following warmer summers and during winters with high predation rate (Figure 2 and Table 4). More precisely, the model identified as having the lowest AICc in the population-level analysis, included three predictors, predation, temp JUL and wolf lag , and explained 79% of the variance in hair loss. However, that model performed equivalently to a simpler model that included only temp JUL and predation ( AICc = 1.38) and which explained 75% of the variance in hair loss (Table 4). Moreover, although multicollinearity is not a concern for the bivariate model because temp JUL and predation were not correlated, it is a concern for the tri-variate model given that predation and wolf lag are highly correlated (see Supplementary Table 1). Additionally, examination of univariate models suggests that although temp JUL and predation are useful predictors of hairloss, wolf lag explains little variation in hairloss. For these reasons, we considered the bivariate model including only temp JUL and predation to be the most parsimonious model for predicting temporal variation in hair loss. Plots of model residuals and formal tests indicate that assumptions about homoscedasticity and normally distributed errors were not violated (Breusch-Pagan test: 0.65, p = 0.42; Shapiro-Wilk: 0.96, p = 0.64; Kolmogorov-Smirnov: 0.13, p = 0.88), and no observations had unduly high leverage.
Because tick abundance was observed to be positively associated with moose abundance in other study systems (Samuel, 2007;Bergeron and Pekins, 2014;Ball, 2017), but TABLE 3 | Mixed-effect models assessing how the proportion of winter coat hair lost or damaged (an indicator of tick burden, see Figure 1 and Supplementary  Figure 1) for individual moose varied between different age-classes (yearling and adult) and sex (bulls and cows) for a population in Isle Royale National Park over a 19-year period (2001-2019 Year was fitted as a random effect.  Table 4. More specifically, the lines in (A) represent predictions when predation rate was fixed at the median value (0.09) and in (B) when mean maximum temperatures were fixed at the median value (25.44 • C). Gray shaded area indicates the 95% upper and lower confidence intervals around predictions.
variables of moose density were not included in the most parsimonious model, we built two ad hoc models, one included moose and the other included moose lag . The better of the two models included moose density with no time lag and explained 28% of the variance in hairloss with hairloss tending to decline as moose density increased (Table 4 and Supplementary Figure 7).   Results of the individual-level analysis (using mixed-effect models) were very similar to the population-level analysis. More precisely, the most parsimonious model included temp JUL and predation (see Supplementary Table 2). Furthermore, model coefficients suggest that hairloss was greater following warmer summers and during winters with high predation rate (temp JUL : 0.24 ± SE 0.04; predation: 3.98 ± SE 0.71).
Lastly, average hair loss exceeded 55% in six of the 19 years which may represent epizootic events. Five of these six putative epizootic events occurred following summers where July temperatures were in the upper 75% percentile of observations. Additionally, five of these six putative epizootic events occurred in years when predation rate was in the upper 75th percentile.

DISCUSSION
Summer temperatures were associated with higher tick burdens the following spring, as indicated by hair loss, with July temperatures explaining nearly half the interannual variation in hair loss (Figure 2A and Table 4). The underlying mechanism is probably that warmer temperatures promote faster egg development and increased egg survival (Drew and Samuel, 1987;Ogden et al., 2004;Lysyk, 2014). That result (Figure 2A) is consistent with general expectations that warmer temperatures are likely to favor higher parasite burdens (Patz et al., 2003;Pounds et al., 2006;Barber et al., 2016).
Prior to performing this analysis, we assumed that other weather-related hypotheses ( Table 1) were equally plausible. For example, earlier studies hypothesized that shorter winters might allow for a longer autumnal questing period and thereby result in higher tick burdens for moose (Jones et al., 2019;Ellingwood et al., 2020). However, we found no evidence to suggest that hair loss was correlated with either the length of the previous winter (winter), mean annual snow depth during the previous winter (snow), or the date that temperatures first drop below freezing for an extended period of time (autumn, see Supplementary Table 1).
Because we observed interannual variation in hair loss over two decades, these results (Table 4 and Figure 2) represent relevant clues for anticipating the influence of climate change on moose-tick interactions. In particular, across most areas of North America, climate warming has been least pronounced during summer and most pronounced during the winters (Vincent et al., 2015;Vose et al., 2017). Indeed, average temperatures in July have not increased significantly over the study period or since the mid-20th century (Supplementary Figure 8). Therefore, our findings suggest that tick burdens are more closely associated with aspects of climate that are changing less rapidly across most regions of North America. One notable exception is that summer temperatures have increased in Maine and New Hampshire and winter tick epizootic events are being observed more frequently for moose populations in this region (Jones et al., 2019). Although some moose populations at the southern limit of the species geographic range have exhibited long-term declines [e.g., in northern Minnesota, New Hampshire, and north-western Ontario (Lankester, 2010;Lenarz et al., 2010;Jones et al., 2019)] other populations have increased or remained stable [i.e., in IRNP, southern Ontario and Utah (Murray et al., 2012;Ruprecht, 2016)]. Those observations suggest that the influence of climate change on moose populations is importantly mediated by regionspecific factor(s).
Tick burden was also observed to be positively associated with predation rate ( Table 4, Supplementary Table 1, and Figure 2B). That positive relationship could have arisen due to a number of different mechanisms. First, it plausible that a positive relationship between predation rate and tick burden may be caused by moose exhibiting risk-sensitive habitat selection during years with high predation risk if it increases the use of habitats (during autumn) where questing ticks are more common. More precisely, in years when predation risk is high, it may be preferable for moose to use habitats where the risk of being killed by wolves is low, even when the risk of being parasitized by ticks is high -given that winter ticks typically do not kill their hosts. While risk-sensitive habitat selection has been observed in Isle Royale moose (Montgomery et al., 2013(Montgomery et al., , 2014, it is unclear whether it occurs in a manner that results in more exposure to questing larvae. This mechanism would represent a case where predation exacerbates the impacts of parasitism. Any such exacerbating influences are liable to have important consequences because of the role that parasites play in determining the composition, structure, and stability of ecosystems (Minchella and Scott, 1991;Poulin, 2002a,b, 2005;Wood and Johnson, 2015). Therefore, one might expect the stability of host-parasite systems to differ when predators are also present.
Second, a positive association between tick burden and predation rate could also be indicative of some underlying density-dependent process. The reason to think this is that predation rate has a strong tendency to be inversely related to moose density (see Supplementary Table 1; Vucetich et al., 2011) and hair loss is negatively correlated with moose density ( Table 3 and Supplementary Figure 9). Additionally, inspection of temporal trends in hair loss, predation rate and moose density suggest that with the exception of the first few years of the study period, hair loss, and predation rate (and wolf abundance) steadily decline over the study period, whilst moose abundance steadily increases (Supplementary Figures 4, 9). However, if changes in moose density (as opposed to predation rate) were the primary factor influencing tick burdens, then one might expect the predictors moose and moose lag to be a significantly better predictors of hair loss than predation, but this was not the case ( Table 3 and Supplementary Material 2).
Third, it is plausible that the direction of the causal relationship is reversed and that higher tick burdens may result in higher predation rates by increasing the vulnerability of moose. Evidence for that mechanism includes high tick burdens being known to significantly reduce body condition of moose (Pekins, 2020). Furthermore, predators are known to disproportionately kill prey in substandard condition (Mech and Boitani, 2003;Wright et al., 2006;Wilmers et al., 2020;Hoy et al., 2021) and wolf per capita kill rates also tend to be higher when vulnerable prey are more common (Sand et al., 2008(Sand et al., , 2012. More generally, parasites have been observed to make hosts more vulnerable to predators in other parasite-host/predator-prey systems (Hudson et al., 1992;Møller and Nielsen, 2007). The extent to which this mechanism operates is the extent to which parasitism exacerbates the impact of predation (by predisposing moose to predation who might otherwise have survived without high tick burden). If so, then the influence of parasitism on host dynamics may be partially masked by the apparent influence of predation in places where predation is an important influence.
Observing a negative correlation between tick burden and host density (see Table 4 and Supplementary Figure 9) is consistent with general expectations of an "encounter-dilution effect" occurring when parasites wait for hosts in the environment (Mooring and Hart, 1992;Mugabo et al., 2015). Plausible as that mechanism may seem, tick burden was positively associated with moose abundance in other study systems (Samuel, 2007;Bergeron and Pekins, 2014;Ball, 2017). This discrepancy between our result and previous assessments may be due to those other studies occurring in systems where moose are subject to humancaused mortality (i.e., hunting and culls), but are not exposed to wolf predation. Hunters tend to be less selective for vulnerable prey than wolves (Wright et al., 2006) and hunter mortality tends to be positively density-dependent (whereas wolf predation is inversely density-dependent). Those differences in cause-specific mortality for moose -and any concomitant effects on moose habitat usage -may account for the discrepancy between our results and previous assessments.
Overall, our results (Table 4 and Figure 2) provide unique insights about interrelationships among climate change, parasites-host dynamics and predator-prey dynamics. Although summer temperatures have not increased dramatically over the last half-century, they are expected to increase significantly in the future. For example, under a high emissions scenario, average summer temperatures are predicted to increase 4-8 • C by the end of the century in the Great Lakes region (Hayhoe et al., 2010). Even under a low emission scenario, summer temperatures are expected to increase 2-4 • C by the end of the century in this region (Hayhoe et al., 2010). If such temperature increases occur, then our results (Figure 2A) suggest that average hair loss may regularly exceed 50-60% which is indicative of severe tick burdens for moose (Samuel, 2004). The ramifications of any such climate related changes in parasite dynamics are likely to be substantial and widespread given that parasites play an important role in shaping the composition, structure, productivity and stability of communities and ecosystems (Minchella and Scott, 1991;Poulin, 2002a,b, 2005;Wood and Johnson, 2015). Lastly, we suggest that future research might profitably focus on better understanding the mechanisms underlying the relationship between parasitism and predation. Such research is likely to become especially valuable if apex predators, such as wolves, continue to increase in abundance and recolonize their former ranges.

DATA AVAILABILITY STATEMENT
The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.

ETHICS STATEMENT
The animal study was reviewed and approved by the IACUC Committee at Michigan Technological University. Specifically, protocol number L0093R includes monitoring moose populations by taking photographs of moose during the spring/summer. for RP also came from the Robbins Chair in Environmental Sustainability. This work was also supported by private donations to the Michigan Tech Fund.

ACKNOWLEDGMENTS
We are grateful to the many individuals who contributed to the collection of field data, especially Erik Freeman. We are also grateful to the U.S. National Park Service for permitting us to conduct research in Isle Royale National Park.