2.1 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 five red deer (Cervus elaphus) stags were kept (32.5 km-2). Deer were supplementary fed during the winter and treated orally against intestinal worms with a broad-spectrum anti-parasite agent (Ivermectin) once a year in the autumn. Due to the low level of treatment and its timing late in the year, it would have had minimal effects on tick burdens on deer (30,31). Within this moorland enclosure, eight plots (four fenced/four unfenced) of 0.23 ha were set up in 2004/2005 (fence mesh size: 20 x 20 cm). 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 (SI 1). As part of a different study, each plot was divided into five subareas each of approximately 15m x 15m to create five 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. To our knowledge, the only hosts permanently present at this site were small mammals and birds; other mammals such as hares, mustelids, foxes or badgers were either not present or rare.
2.2 Quantifying questing ricinus nymph density (DON)
Questing nymphs were surveyed in 2013, 2014, 2018 and 2019 using a standard dragging method (32) 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. To reduce edge effects, ticks were surveyed at least 2m away from the fence and ticks tend not to be able to walk horizontally more than this distance (33). Air temperature and relative humidity were recorded for each transect to be taken account of in statistical models as these factors affect tick activity (34). 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.
1.3. Quantifying rodent activity and tick burdens
To test the effect of vegetation on rodents and the effect of rodents on NIP (M2- dilution mechanism and M3 – ecological cascade mechanism) as well as to quantify tick burdens on rodents (M1- increase in transmission mechanism), 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 five habitats within each plot, at least 2m away from the edge of the 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 four of the plots (two fenced/two unfenced) and for two consecutive nights in the remaining four plots. 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 (M2) 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.
1.4. Measurement of B. burgdorferi s.l. prevalence (NIP)
As part of testing all three mechanisms 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 (35). 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 (36). 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 (37) was optimized using the IQTM Supermix (Bio-Rad Laboratories, Hercules, USA) in a Stratagene Mx3005P thermal cycler (Agilent, Santa Clara, US). Each reaction contained IQTM 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 confirmed 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).
2.1.1. 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 inflation factors (VIFs) and variables for which the VIF was above 4 were discarded from the model (38). Model selection was done using the dredge function from the MuMIn package (39) based on the corrected Akaike Information Criterion (AICc) (40). 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) with B. burgdorferi s.l.
Our three mechanisms are all concerned with a potential effect of high deer density on 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). An observation level random effect was included to account for overdispersion (41).
Mechanism 1: increase in transmission potential
We investigated whether high deer density led to increased DON, which we hypothesised to be a consequence of high deer density resulting in higher questing larval and nymph density and higher tick burdens on rodents. We used a hurdle generalized linear mixed effect model (GLMM) with a Poisson distribution. We tested for zero-inflation using the DHARMa package (42) and we chose a hurdle model as the number of zeros were underestimated with a Poisson GLMM (43). 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. Plot number and an observation level random effect were included, the latter to account for overdispersion (44,45).
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 M3) so this analysis could not be performed.
Mechanisms 2: dilution effect
The dilution effect mechanism 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 M3), which precluded us from fitting 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 mechanism, we therefore report descriptive statistics for tick burdens on rodents and relative host abundance.
Mechanism 3: ecological cascade
To test the first part of the cascade, that high deer density negatively affects ground vegetation height (M3), 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 fixed effects and a random effect of plot number.
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 first, the full model included vegetation height, habitat and year as fixed covariates and a random effect of plot number. 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 (17). For the second analysis, testing for a direct impact of deer on rodent activity, the full model included deer treatment, habitat and year as fixed covariates and a random effect of plot number.
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 included an observation level random effect to account for overdispersion (41).
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 (43). 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 surveyed. The full model included deer treatment, month, year, habitat, temperature and whether the ground was dry during tick survey as fixed variables and the plot number and an observation level random effect were included as random terms.