Experimental Evidence for Opposing Effects of High Deer Density on Tick-Borne Disease Prevalence and Hazard

Identifying the mechanisms driving disease risk is challenging for multi-host pathogens, such as Borrelia. burgdorferi s.l., the tick-borne bacteria causing Lyme disease. Deer are tick reproduction hosts but do not transmit B. burgdorferi s.l., resulting in potentially opposing effects on transmission. Here, we use a deer exclosure experiment to test three hypotheses for how high deer density shapes B. burgdorferi s.l. prevalence in ticks: (H1) high transmission on rodents due to higher tick densities; alternatively, (H2) low B. burgdorferi s.l. prevalence because more ticks feed on deer rather than transmission-competent rodents (dilution effect); (H3) ecological cascades, whereby lower vegetation decreases rodent abundance thus reducing transmission. Although we found support for all three mechanisms, prevalence was reduced almost 3-fold in high deer density plots compared to exclosures, suggesting that the dilution (H2) and cascade (H3) mechanisms outweighed the increased opportunities for transmission (H1). High deer density led to lower vegetation and fewer rodents, providing evidence for an ecological cascade. However, Lyme disease hazard (density of infected ticks) was increased 5-fold at high deer densities due to an 18-fold rise in tick density. This demonstrates that reproduction hosts like deer can drive up vector-borne disease hazard at high densities, despite simultaneously reducing pathogen prevalence.


Introduction
Ecological trophic cascades occur when a change in population or activity of a keystone species precipitates a cascade of alterations through trophic levels in an ecosystem (Paine 1980). For example, in coastal kelp forests in the Northeast Paci c, sea otters (Enhydra lutris) prey on herbivorous sea urchins (Strongylocentrotus spp.), which releases grazing pressure on kelp and has further consequences for benthic plant and animal communities (Carter et al. 2007; Watson and Estes 2011).
While there are many examples for trophic cascades in ecology, tests of their extension to infectious disease risk have received less attention. This concept could be particularly signi cant for vector-borne diseases, which rely upon vectors that carry and transmit pathogens from one host to another, as trophic effects could manifest from changes in both vector and reservoir host populations. For instance, one study showed that the introduction of Burmese pythons (Python bivittatus) resulted in a decline in some mammal species but not in those acting as transmission hosts for the Everglades virus. This affected the feeding behaviour of mosquitoes and increased blood meals taken from transmission hosts (Hoyer et al. 2017). Other researchers found a positive association between coyote (Canis latrans) density and Lyme disease incidence because coyotes have negative effects on red fox (Vulpes vulpes) populations, which caused an increase in small mammals that transmit Lyme disease pathogens (Levi et al. 2012). While these studies show how a change in predators can affect vector-borne pathogen risk through predation of transmission hosts, no study to our knowledge has tested for trophic cascades effects on pathogens arising from changes in large herbivore populations that alter habitat (e.g. vegetation structure) for transmission hosts. It is already established that large herbivores alter habitat (Buesching et al. 2011), with cascading effects on small mammal communities (Flowerdew and Ellwood 2001; van Wieren and Bakker 2008), which are transmission hosts for many pathogens, including tick-borne pathogens (e.g. Tick-borne encephalitis virus and the bacteria causing Lyme disease). This potential mechanism of vector-borne disease risk seems likely and requires testing since it has not been studied before.
Lyme disease, the most prevalent vector-borne disease in humans in the northern hemisphere, is an emerging disease in Europe and North America (Stanek et al. 2012). It is caused by the Borrelia burgdorferi sensu lato complex of bacteria (hereafter referred to as B. burgdorferi s.l.) and transmitted by Ixodid ticks. In Europe, the primary vector is Ixodes ricinus (Stanek et al. 2012), a generalist three-host tick that feed on a wide range of vertebrate species (Hofmeester et al. 2016). There are several genospecies of the bacteria within the B. burgdorferi s.l. complex, each with different transmission host associations: in Northern Europe, rodents transmit B. afzelii (Hanincová et al. 2003) while birds are associated with B. valaisiana and B. garinii (Heylen et al. 2014). B. burgdorferi sensu stricto, the fourth most abundant genospecies, can be transmitted by both mammals and birds. Transmission is therefore predicted to increase as more immature ticks feed on infected hosts such as rodents (Daniels and Fish 1995;Perkins et al. 2006 H1-Increase in transmission potential high deer density will produce a high density of questing nymphs (DON) and thus a higher tick burden on transmission hosts (Daniels and Fish 1995) and higher NIP in questing ticks than areas without deer.
H2-Dilution effect lower prevalence in areas of high deer density due to a higher proportion of larval ticks feeding on deer instead of transmission hosts such as rodents. prevalence. Rodent activity might also be affected by high deer density directly through disturbance.
These different mechanisms are not mutually exclusive and the effect of high deer density on NIP will depend on the relative strength of the rst (transmission potential) hypothesis compared to the other two (dilution and ecological cascade hypotheses) (Fig. 1). A further key aim was to test the impact of high deer density on Lyme disease hazard, which is de ned as the density of infected nymphs (DIN) in the environment and is the product of NIP and DON. It is di cult to predict the effect of deer on Lyme disease hazard because it will depend on the relative strengths of the effect of deer on DON and the mechanisms driving NIP (Fig. 1).
The use of experimental deer exclosures is particularly suitable for testing the ecological cascades hypothesis as it maximises habitat impacts while avoiding noise being introduced from heterogeneities in land use, habitat, topography and climate that typify landscape-scale surveys. In addition, the exclusion of deer has high applied relevance since deer fencing is a common land management tool and is increasingly being used in mitigating the impacts of ticks (Gilbert et al. 2012).

Experimental design
The experiment took place at Glensaugh research farm in Aberdeenshire, Scotland (56.914217N, -5.532070E), in a 15 ha enclosure of upland moorland where ve red deer (Cervus elaphus) stags were kept (32.5. km − 2 ). Within this moorland enclosure, four replicates of fenced-unfenced pairs of plots of 0.23 ha were set up in 2004/2005. The four fenced plots (hereafter referred to as deer exclusion plots) excluded deer while the four unfenced plots (hereafter referred to as high deer density plots) were accessible to the deer and subjected to high grazing pressure. As part of a different study, each plot was divided into ve subareas each of approximately 15mx15m to create ve habitats: high density birch, low density birch, a single birch in the centre of the plot, high density pine, and a control which was not planted and consisted of Calluna vulgaris-dominated heathland. Tree saplings were 9 years old by the time this study started in 2013.

Quantifying questing I. ricinus nymph density (DON)
Questing nymphs were surveyed in 2013, 2014, 2018 and 2019 using a standard dragging method (Falco and Fish 1992) which consists of dragging a 1m x 1m square of woollen blanket material over the ground vegetation for 10m linear transects. In 2013 and 2014, 10 transects were surveyed in each of the heather, high density birch and pine habitats for all plots (n = 720 transects/year) while in 2018 and 2019, six transects were surveyed in each of the same three habitats (n = 432 transects/year). Tick surveys were done in May, July and September between 0900 hours and 1900 hours. Air temperature and relative humidity were recorded for each transect to be taken account of in statistical models as these factors affect tick activity (Gilbert 2010). Ground vegetation height was recorded at the beginning (0m), middle (5m) and end (10m) of each transects using a sward stick. Questing nymphs were counted, collected and stored at -20°C for later pathogen analysis.

Quantifying rodent activity and tick burdens
To test the effect of vegetation on rodents and the effect of rodents on NIP (H2-dilution hypothesis and H3 -ecological cascade hypothesis) as well as to quantify tick burdens on rodents (H1-increase in transmission hypothesis), rodent activity was estimated in July 2017 and 2018 using a live-trapping method. In 2017, four non-selective Sherman live traps (16 x 5 x 6.5 cm. Sherman Inc., Tallahassee, Florida) were baited with oats and placed in each of the ve habitats within each plot with 10m spacing, resulting in a total of 20 traps per plot (i.e. 200 trap nights (TN) per deer treatment). Traps were set for three consecutive nights in two of the paired fenced/unfenced plots (replicates 1 and 2) and for two consecutive nights in the other two paired fenced/unfenced plots (replicates 3 and 4). In 2018, four traps were installed for four consecutive nights in the three habitats in which ticks were surveyed (heather, high density birch and pine habitats) in each plot, resulting in 12 traps per plot (i.e. 192 trap nights (TN) per high deer treatment). Traps were activated after 1600 hours and checked every morning before 1000 hours. For all captures, species, sex, weight and approximate age (juveniles or adult) of each individual were recorded. Ticks attached to rodents were counted and collected from around the head and ears and stored in 100% ethanol. All captured rodents were released at the capture site. For analysis, we took the proportion of traps with a capture as our rodent activity index. Both tick burden and rodent captures data were needed to test the dilution effect (H2) which requires a comparison between high deer density plots and deer exclosures in the relative proportion of the larval tick population that feeds on rodents versus deer.

Measurement of B. burgdorferi s.l. prevalence (NIP)
As part of testing all three hypotheses of the effects of high deer density on NIP, and to estimate Lyme disease hazard (DIN), DNA was individually extracted from questing nymphs using an ammonia extraction (Gern et al. 2010). B. burgdorferi s.l. was detected using two methods for this study: nested PCR and qPCR. Ticks collected in 2013 and 2014 were tested using a nested PCR targeting the 5S-23S intergenic spacer region using the protocol described by (Rijpkema et al. 1995). Following an issue arising from the nested PCR protocol, samples collected in 2018 and 2019 were tested using a qPCR method. A qPCR protocol on fragments of OspA genes (Heylen et al. 2013) was optimized using the IQ™ Supermix (Bio-Rad Laboratories, Hercules, USA) in a Stratagene Mx3005P thermal cycler (Agilent, Santa Clara, US). Each reaction contained IQ™ Supermix, two primers at 200nM (B-OspA_modF: AATATTTATTGGGAATAGGTCTAA and B-OspA_borAS:-CTTTGTCTTTTTCTTTRCTTACAAG), the probe (B-OspA_mod-probe:-FAM-AAGCAAAATGTTAGCAGCCTTGA-BHQ-1™) at 100nM and 3µL of DNA. One positive and one negative control were added for every plate.
We con rmed the correspondence between the two PCR protocols by using 61 known positive samples (by nested PCR, mix of genospecies) and 344 known negative samples. In all cases, results from the qPCR matched those of the nested PCR. Positive samples from the qPCR protocol were subjected to the nested PCR protocol to identify B. burgdorferi s.l. genospecies. All PCR products from positive samples were Sanger sequenced to identify the genospecies of B. burgdorferi s.l. For analysis, we used the proportion of questing nymphs infected (NIP).

Statistical analyses
All statistical analyses were performed in R version 3.5.1 (R-core-team, 2013). For Generalized Linear Mixed Effects Models (GLMMs) and Generalized Linear Models (GLMs), we tested for potential collinearity between explanatory variables by calculating variance in ation factors (VIFs) and variables for which the VIF was above 4 were discarded from the model (Zuur et al. 2009). Model selection was done using the dredge function from the MuMIn package (Barton 2019) based on the corrected Akaike Information Criterion (AICc) (Brewer et al. 2016). When appropriate, we conducted post hoc Tukey tests to assess pairwise comparisons between levels of categorical variables.
The effects of high deer density on nymphal infection prevalence (NIP) withB. burgdorferis.l.
Our three hypotheses are all concerned with potential mechanisms through which high deer density might affect NIP. To test for this effect of high deer density, we used a binomial GLMM with a logit link and NIP for B. burgdorferi s.l. (3 prevalence estimates per habitat and per year) as the response variable. The full model included deer treatment (deer exclusion or high deer density) as our main predictor, as well as month (May, July or September), habitat (high density pine, high density birch or heather control) and year (2013, 2014, 2018 or 2019). Random effects of a combined plot and habitat parameter (e.g. high density birch in fenced replicate 1) and an observation level random effect were included, the latter to account for overdispersion (Harrison 2015).

Hypothesis 1: increase in transmission potential
To test the hypothesis that high deer density results in high questing nymph density, which could result in higher tick burdens on rodents, we investigated the effects of high deer density on DON. We used a hurdle generalized linear mixed effect model (GLMM) with a Poisson distribution. We tested for zero-in ation using the DHARMa package (Hartig 2020) and we chose a hurdle model as the number of zeros were underestimated with a Poisson GLMM (Lewis et al. 2011). The response variable was the number of questing nymphs collected per 10m transect and the full model included deer treatment, month, year, habitat, ground vegetation height, relative humidity, temperature and whether the ground was dry during collection. We also included the interactions between deer treatment and month and between deer treatment and habitat, as the effects of deer on ticks might vary between months and habitat type. A random effect for each plot habitat combination (e.g. high density birch in fenced replicate 1) and an observation level random effect were included to account for overdispersion (Elston et al. 2001;Harrison 2014).
Although we initially planned to compare tick burdens on rodents between deer density treatments, only two rodents were caught in high deer density plots (see H3) so this analysis could not be performed.

Hypothesis 2: dilution effect
The dilution effect hypothesis implies that the relative proportion of the larval tick population feeding on deer vs rodents (i.e. tick burden x relative abundance of each host type) is higher in high deer density plots compared to exclosures. Testing this formally was not possible because only two rodents were captured in high deer density plots (as predicted by H3), which precluded us from tting a model of the effects of high deer density on tick burden on rodents. We did not count tick burdens on deer. To evaluate this hypothesis, we therefore report descriptive statistics for tick burdens on rodents and relative host abundance.

Hypothesis 3: ecological cascade
To test the rst part of the cascade, that high deer density negatively affects ground vegetation height (H3), we used a GLMM with a Gaussian distribution and ground vegetation height as our response variable. The full model included deer treatment, month and year as xed effects and a random effect of each plot habitat combination.
To test the second part of the cascade, that lower vegetation negatively affects rodent activity, we conducted two analyses, both using a binomial GLMM with a logit-link and our rodent index as the response variable. For the rst, the full model included vegetation height, habitat and year as xed covariates and a random effect of each plot habitat combination. Vegetation height was averaged per habitat type and per year and we used vegetation height from the heather habitat for the single birch and sparse birch habitats, because there was no difference between these three (Gilbert et al. 2012). For the second analysis, testing for a direct impact of deer on rodent activity, the full model included deer treatment, habitat and year as xed covariates and a random effect of a random effect of each plot habitat combination.
To test the last part of the cascade, that rodent abundance positively affects NIP, we used a binomial GLMM. The response variable was NIP for B. afzelii and the full model included the number of rodents captured the previous year (data for 2017 and 2018), habitat, year and month. We added the combined plot and habitat parameter variable as a random effect and an observation level random effect was included to account for overdispersion (Harrison 2015).

The effects of high deer density on Lyme disease hazard
To test how high deer density affects Lyme disease hazard (DIN), we used a hurdle GLMM with a Poisson distribution, as the number of zeros was underestimated with a Poisson GLMM (Lewis et al. 2011). The response variable was the density of infected nymphs and, as this variable was averaged at the plot level (3 estimates per habitat and per year), we used an offset for the area to transform the values to integers (Zuur et al. 2009). The full model included deer treatment, month, year, habitat, temperature and whether the ground was dry during tick survey as xed variables and the combined plot and habitat level variable and an observation level random effect were included as random terms.

Results
Over four years of data collection, 12,310 questing nymphs were sampled (Table 1). A random subset of 1,000 nymphs was examined under the microscope for species identi cation using speci c keys (Hillyard 1996) and all were found to be Ixodes ricinus. Thus, we assumed that all questing ticks collected in this study were I. ricinus. We tested all 585 questing nymphs collected from deer exclusion plots and a subset of 1,042 from high deer density plots for B. burgdorferi s.l. Prevalence (NIP) was 0.9% [95%CI: 0.4-1.6] in high deer density plots and 2.7% [95%CI: 1.6-4.4] in deer exclusion plots ( Table 1). The dominant B. burgdorferi s.l. genospecies was the rodent-associated B. afzelii (92%, n = 23/25) followed by the birdassociated B. valaisiana (4%, n = 1/25) and B. garinii (4%, n = 1/23). Over two years of trapping, 24 individual bank voles were caught in deer exclusion plots (6.69/100TN, SD = 2.41), while only two were captured in high deer density plots (0.51/100TN, SD = 0.97) and no other rodent species were detected (Table 1). Table 1 Mean values across all plots and habitat-sub areas, separated by year and deer treatment (high deer density and deer exclusion plots), for the density of questing nymphs (DON), number of bank voles caught per 100 trap nights, prevalence of B. burgdorferi s.l. in questing nymphs (NIP) and density of nymphs infected with B. burgdorferi s.l. (DIN). Ticks were not collected in 2017. Rodent trapping was conducted in 2017 and 2018 only. SD -standard deviation; 95% CI − 95% con dence interval. compared to deer exclusion plots (2.2%, 95%CI: 0.6-7.6) (Fig. 2). NIP also varied across years, habitats and months and the full model results are provided in the Supplementary Information 1.

Hypothesis 1: increase in transmission potential
Effects of high deer density on DON: The nal model included deer treatment, month, habitat type, year, relative humidity, temperature and whether the ground was wet. On average, predicted DON was 18 times higher in high deer density plots (DON = 10.8, 95%CI: 0.1-36.2) compared to deer exclusion plots (DON = 0.6, 95%CI: 0-3.5) (Fig. 3). DON varied across months, years and habitats and was in uenced by temperature, relative humidity and whether the ground was wet (SI 1).
Effects of high deer density on tick burdens on rodents: Bank voles in deer exclusion plots (n = 24) harboured, on average, 14.2 larvae (SD = 11.66) while bank voles in high deer density plots (n = 2) had 17.5 larvae (SD = 10.61). Low captures in high deer density plots prevented us from statistically comparing tick burdens between treatments further. No nymphs were found attached to a bank vole.

Hypothesis 2: Dilution effect
Comparison between high deer density plots and exclosures on the relative proportion of the larval population feeding on rodents vs deer: In deer exclusion plots, bank voles fed 340.8 larvae (24 individuals feeding 14.2 larvae each on average). In high deer density plots, bank voles fed 35 larvae (2 individuals feeding 17.5 larvae each on average). These data suggest that 9.7 times more larvae fed on bank voles in deer exclusion plots compared to high deer density plots. In deer exclusion plots, no larvae fed on deer (there were no deer in exclosures). High deer density plots contained a deer density of 32.5. km − 2 but we were not able to measure tick burdens on these deer.  (Fig. 4A).
Ground vegetation height was also in uenced by month and habitats (SI 1).
Effects of vegetation height on rodents: The nal model included deer treatment and year and the predicted number of bank voles captured was positively correlated with ground vegetation height and increased by 1.17 captures per 100 TN for every 1 cm increase in vegetation height (Fig. 4B). Bank vole capture rate also varied across both years (SI 1).
Effects of rodent capture rate on NIP for B. afzelii: While the best model only included year and month, rodents were included as predictors in the second-best model (ΔAICc = 1.9) and predicted NIP increased by 0.01% for every vole caught per 100 TN (Fig. 4D) (SI 1).

Discussion
A replicated paired fencing experiment allowed us to test three hypotheses (increased transmission potential, dilution and ecological cascades) for how high deer density may affect B. burgdorferi s.l. infection prevalence in questing ticks (NIP). Furthermore, we could test how these effects on NIP interact with the role of deer in supporting tick populations and shaping Lyme disease hazard. We found that nymphal infection prevalence was lower in high deer density plots along with evidence pointing to a combination of dilution and trophic cascading effects. Despite the lower NIP, Lyme disease hazard (DIN) was ve times higher in high deer density plots compared to deer exclusion plots, due to a strong, positive association between deer and nymph density.
Our transmission potential hypothesis (H1) predicted that a high deer density would result in an increased tick density, causing higher tick burdens (and therefore increased transmission potential) on individual rodents. We did nd a strong positive effect of high deer density on DON, highlighting the importance of In contrast however, we found high deer density plots had half the predicted NIP of exclusion plots, suggesting that increased tick densities and transmission were counter-acted by other factors driving NIP down.
Several studies have reported positive correlations between DON and tick burdens on rodents (Daniels and Fish 1995;Hofmeester et al. 2017a) and other hosts (Gilbert et al. 2017). This may also have been the case for rodents in our experiment but, as only two bank voles were captured in high deer density plots (as predicted by the ecological cascade hypothesis), we were not able to test for differences between deer treatments in individual rodent tick burdens. This dramatic reduction in rodent captures in high deer density plots is a key nding and invalidates H1's assumption of similar rodent densities with deer density. Therefore, even if individual rodents had higher tick burdens in high deer density areas, rodent density was probably too low for maintaining effective pathogen transmission in this environment.
Our observation that B. burgdorferi s.l. prevalence was almost three times higher in deer exclusion plots, despite estimates being associated with considerable uncertainty, is consistent with predictions from both the dilution effect (H2) and ecological cascade (H3) hypotheses. The dilution effect predicts lower NIP at high deer densities due to a lower proportion of the larval tick population feeding on rodents, which are transmission hosts for B. burgdorferi s.l., than on deer that do not transmit the pathogens (Gray et al. 1992;Ostfeld and Keesing 2000;Vourc'h et al. 2016). This is challenging to test because, ideally, it requires density estimates and tick burdens of the key hosts. We could not obtain tick burden measurements on deer and we caught only two voles in high deer density plots, precluding statistical test.
However, based on differences in vole abundance between deer treatments, we estimated that voles fed almost 10 times more larvae in deer exclusion plots compared to high deer density plots. This compares to no larvae feeding on deer in exclusion plots (as deer were absent), in contrast to the situation in high deer density plots where it is highly likely that most larvae fed on deer. This is due to the (i) very high deer density (32.5. km − 2 ) and (ii) high densities of questing nymphs (a result of high larval survival) in high deer density plots. We therefore suggest that a dilution effect was one of the mechanisms operating on NIP in response to high deer densities, as it has been previously suggested for B. burgdorferi s.l. (Gray et al. 1992 . It is therefore likely that the populations of alternative hosts were low in high deer density plots and that the majority of larvae must have been feeding on red deer. In addition, the fact that almost all nymphs that tested positive (92%) were infected with the rodent associated pathogen, B. afzelii, suggests that other hosts (e.g. birds) did not contribute much to B. burgdorferi s.l. transmission in our system.
To the best of our knowledge, this is the rst test of the ecological trophic cascade hypothesis of B. burgdorferi s.l. prevalence, predicting that grazing pressure from high deer densities will reduce the vegetation, resulting in fewer rodents and therefore lower NIP. We found support for most of the expected trophic links: high deer density resulted in shorter vegetation, highlighting the effects of deer grazing on  (Fig. 4D). The weak association between vole abundance and NIP in this experiment could be due to several factors, including low prevalence of B. burgdorferi s.l. overall, providing insu cient signal for detecting unequivocal statistical effects. Another contributing factor might have been the small spatial scale of our experiment plots (0.23 ha), facilitating likely movements of rodents between fenced and unfenced plots, which were 30-100 m apart.
The small spatial scale of our plots, while a potential issue for con rming a link between rodent abundance and NIP, proved su cient to con rm strong effects of high deer density with DON, NIP, Lyme and Bakker 2008), such a combination may not commonly occur in nature. This demonstrates the importance of taking a systems approach, including considering potential ecological cascades when predicting disease risk. While our experiment used extreme contrasts in deer densities (zero vs 32.5. km − 2 ) it is notable that we were able to detect effects on infection dynamics even on very small spatial scales (0.23 ha plots).
Our results highlight the need for taking a systems approach to such a complex disease system where different host types might affect each other's densities through habitat modi cation, or other means such as predation (Levi et al. 2012). Gaps in knowledge that can now be addressed include testing for cascading effects from deer on B. burgdorferi s.l. prevalence at landscape scales, and including a full range of intermediate deer densities, to investigate non-linearities and thresholds of effects of deer and to test whether the mechanisms supported here operate in a non-experimental setting.
In conclusion, we found that high deer density could lower Lyme disease pathogen prevalence by a combination of dilution and trophic cascading effects. Despite this, Lyme disease hazard was ve times higher in high deer density plots due to a strong positive effect of deer on tick density. This study is, to our knowledge, the rst to test for cascading effects of deer on B. burgdorferi s.l. prevalence via grazing and suppressing rodent abundance. This study highlights the need for a systems approach to understand disease dynamics and risks that could arise from complex ecological interactions between host types and habitat.