High burdens of Ixodes scapularis larval ticks on white-tailed deer may limit Lyme disease risk in a low biodiversity setting.

An inverse relationship between biodiversity and human health has been termed the 'dilution effect' paradigm. In the case of tick-borne infections such as Lyme disease, the key assumption is that Borrelia burgdorferi sensu lato abundance is increased by the loss of less competent (dilution) hosts as biodiversity declines. White-tailed deer play a dual role in the pathogen cycle, as key reproductive hosts for adult ticks and incompetent hosts for the pathogen. While the role of deer as hosts of adult ticks is well established, the extent to which deer also feed immature ticks and reduce the proportion infected is unknown because of logistic constraints in measuring this empirically. We estimated the proportion of larvae that fed on deer in an extremely species-poor community on Block Island, RI, where tick nymphal infection prevalence was found to be lower than expected. In 2014, we measured the density, larval tick burdens, and realized reservoir competence of small mammal and bird hosts on Block Island, RI. In 2015, we measured the infection prevalence of host-seeking Ixodes scapularis nymphs resulting from larvae fed on available hosts in 2014. We back-estimated the proportion of larvae expected to have fed on deer in 2014 (the only unknown parameter) to result in the nymphal infection prevalence observed in 2015. Back-estimation predicted that 29% of larval ticks must have fed on deer to yield the observed 30% nymphal infection prevalence. In comparison, the proportion of larvae feeding on mice was 44% and 27% on birds. Our study identified an influential role of deer in reducing nymphal tick infection prevalence and a potential role as dilution hosts if the reduction in nymphal infection prevalence outweighs the role of deer as tick population amplifiers. Because both deer and competent hosts may increase in anthropogenic, fragmented habitats, the links between fragmentation, biodiversity, and Lyme disease risk may be complex and difficult to predict. Furthermore, a nonlinear relationship between deer abundance and Lyme disease risk would reduce the efficacy of deer population reduction efforts to control Lyme disease.

Introduction

The majority of pathogens infect multiple hosts that vary in their ability to acquire and transmit the various pathogens (reservoir host competence). For these multi-host pathogens, the addition of one or more host species to a community can make the pathogen less abundant and less likely to persist than in the presence of reservoir host competent species alone (Begon, 2008). When reduced pathogen abundance results in reduced human disease risk, this effect has been termed the ‘dilution effect’ (Norman et al., 1999; Schmidt and Ostfeld, 2001) or the ‘biodiversity-buffers-disease’ hypothesis (Randolph and Dobson, 2012; McCallum, 2015). The concept of buffering against disease as an ecosystem service of biodiversity has garnered considerable attention, both favorable (Allan et al., 2003; LoGiudice et al., 2003; Keesing et al., 2010; Johnson et al., 2013; Turney et al., 2014; Werden et al., 2014; Civitello et al., 2015) and unfavorable (Randolph and Dobson, 2012; Lafferty and Wood, 2013; Salkeld and Jones, 2013). There is, however, general agreement that the dilution effect may not just be an effect of species diversity per se, but rather be driven by the association between diversity, the specific identity, and relative abundance of competent or incompetent hosts in a community (Keesing et al., 2006).

In the case of tick-borne pathogens, the relationship between pathogen persistence/abundance and host density is complex. Lyme disease is the most commonly reported vector-borne disease in the United States (Centers for Disease Control and Prevention, 2015) and has been the focus of many studies linking host community diversity and human risk (Ostfeld and Keesing, 2000; LoGiudice et al., 2003, 2008). The etiologic agent of Lyme disease in the United States, Borrelia burgdorferi sensu stricto, is transmitted by the blacklegged tick, Ixodes scapularis Say (Acari:Ixodidae) and the Western blacklegged tick, I. pacificus Cooley and Kohls (Acari:Ixodidae) on the West Coast, to a wide range of vertebrate hosts (Donahue et al., 1987; Castro and Wright, 2007; Brinkerhoff et al., 2011). Ixodes scapularis ticks feed once per life stage; the immature stages (larva and nymph) feed on mammalian and avian reservoir hosts of varying levels of competence, as well as white-tailed deer, Odocoileus virginianus (hereafter ‘deer’) which are mostly reservoir incompetent (Telford et al., 1988; Luttrell et al., 1994). Most ticks acquire the pathogen as larvae and transmit it to new hosts as nymphs. Adult I. scapularis, however, depend mostly on deer to mate and for females to obtain a final blood meal (Piesman, 1979; Kilpatrick et al., 2014).

Deer thus play a dual role in the B. burgdorferi sensu lato (hereafter B. burgdorferi) life cycle, they are the reproductive hosts for adult ticks – thus contributing to transmission through increasing tick populations, but they are incompetent hosts for B. burgdorferi – thus reducing infection in the immature ticks. This dual role has been theoretically and empirically explored in European studies of other deer species role in tick-borne encephalitis virus (TBEV) (Bolzoni et al., 2012; Cagnacci et al., 2012) and louping ill virus (LIV) (Norman et al., 1999), which are transmitted by I. ricinus, the vector of B. burgdorferi in Europe (Tälleklint and Jaenson, 1994; Kurtenbach et al., 1995). Using an experimental exclosure, Perkins et al. (2006) inferred the role of deer (roe deer, Capreolus capreolus, and red deer, Cervus elaphus) as dilution hosts because they observed higher infection prevalence within deer exclosures, presumably due to increased larval feeding on mice where deer were removed. A similar phenomenon was observed with exclusion of fallow deer (Dama dama) in Ireland (Gray et al., 1992). Consistent with a dual role of deer, pathogen persistence and the density of infected nymphs were found, both theoretically and empirically, to have humped–shaped relationships with deer density (Norman et al., 1999; Perkins et al., 2006; Rosà and Pugliese, 2007; Bolzoni et al., 2012; Cagnacci et al., 2012). This nonlinear relationship is thought to arise because increases in deer density initially contribute to pathogen amplification by increasing tick abundance up to a critical point. After this threshold, the role of deer in diverting larval tick bites from competent small vertebrate hosts becomes dominant and thus pathogen infection/persistence decreases (Hudson et al., 1995; Dobson and Foufopoulos, 2001). Rosa and Pugliese (2007) referred to this reduction in infection at high deer densities as the ‘dilution effect’. Limited theoretical work has proposed this humped-shaped relationship may also apply to the eco-epidemiology of B. burgdorferi s.s. in the United States (Johnson et al., 2015), but this has not been examined empirically.

We investigated the potential role of deer in feeding larval ticks and reducing nymphal infection prevalence (i.e. the right-hand side of the humped-shaped relationship) on Block Island, Rhode Island, which features a species-poor vertebrate host community strongly dominated by white-footed mice, Peromyscus leucopus, and a low number of bird species (States et al., 2014). In such a host community, the nymphal infection prevalence (NIP), a measure of infection risk, was predicted by modeling studies to be significantly higher than in more biodiverse host communities (LoGiudice et al., 2003; Brisson and Dykhuizen, 2006). Instead, studies from Block Island found a NIP of 28% in 2010 and 24% in 2011 (States et al., 2014), which is similar or lower than most other studies on mainland communities with higher species diversity (Allan et al., 2003; Diuk-Wasser et al., 2012). To assess whether the high deer densities on Block Island contribute to the lower than expected NIP, we conducted a two-year study where we measured host density, larval burdens, and reservoir competence in year 1 (2014) for all relevant small mammal and bird hosts, as well as the density of deer. Larval ticks that fed in year 1 (2014) emerged as nymphs in year 2 (2015) of the study and were sampled to estimate year 2 NIP. We then back-estimated the proportion of larvae expected to have fed on whitetailed deer in year 1 to produce the NIP observed in year 2.

Methods

Field sampling

Study Area

The study was conducted on Block Island, a small island (25.9 km2) located 23 km off the southern coast of Rhode Island, where B. burgdorferi is endemic (Krause et al., 2003). Vegetation is characterized by tall coastal scrub species, deciduous natives, and exotic shrubs, such as Bayberry (Myrica cerifera), Red Maple (Acer rubrum), Black cherry (Prunus serotina), and Multiflora rose (Rosa multiflora) (The Rhode Island Natural History Survey, 2002). Small mammal and bird sampling was conducted in 2014 and host-seeking I. scapularis nymphs – which fed as larvae on hosts sampled in 2014, were sampled in 2015. Small-mammal sampling and collection of host-seeking nymphs occurred at three different locations across the island: the north plot (NP: 41°21′ N, 71°57′ W), south plot (SP: 41°15′ N, 71°58′ W), and midland plot (ML: 41°16′ N 71°58′ W). Birds were mist-netted in two representative habitat types on Block Island: Bayrose (BR: 41°12′ N, 71°33′ W) which is characterized by dense vegetation and relative isolation from human disturbance, and Ocean View Pavilion (OVP: 41°10′ N, 71°33′ W) which is an open, peridomestic habitat characterized by grasses and invasive plant species.

Small-Mammal Sampling

Block Island is characterized by low mammalian diversity, comprised exclusively of the house mouse (Mus musculus), muskrat (Ondatra zibethicus), meadow vole (Microtus pennsylvanicus), Norway rat (Rattus norvegicus), white-footed mouse (P. leucopus), and white-tailed deer (O. virginianus) (Comings, 2006). Only 0.6% of the small animals trapped in our surveys comprised non-P. leucopus species; these species were excluded from our analyses. In 2014, mark-recapture of small mammals occurred at three study plots (NP, SP, ML) and was conducted for three consecutive nights, with each plot visited every other week, totaling seven trapping sessions from late May to late August. Sherman live traps (7.62 × 8.89 × 22.86 cm, H.B. Sherman Traps, Inc. Tallahassee, FL) were positioned 10 m apart in a grid formation at NP (n = 58), SP (n = 110), and ML (n = 60), each trap was baited with peanut butter, oats, and sunflower seeds. Traps were set at dusk and checked the following morning at dawn. All captured white-footed mice were ear-tagged, morphometric data (sex, age, weight, breeding status) were collected, and each individual was checked for ticks across the entire body. Attached ticks were removed with forceps, placed in live collection tubes, and later frozen at −80 °C. A maximum of 10 engorged I. scapularis larva were collected per ear and body when available, the remaining body burden was determined through visual examination. All animal experiments comply with the National Institutes of Health guide for the care and use of laboratory animals (NIH publications No. 8023). All trapping and handling procedures were approved by the Yale University Institutional Animal Care and Use Committee (permit #07596) and Columbia University Institutional Animal Care and Use Committee (permit #AC-AAAL3656).

Bird Sampling

Bird sampling focused on resident species, rather than migratory species, by concentrating sampling between mid-June to late August, when limited migration occurs. At each sampling site (BR and OVP) 4–10 mist-nets (12 × 2.5 m) were located adjacent to one another. Mist netting and all procedures were performed under permit #09636 from the federal Bird Banding Laboratory. Bird sampling was conducted weekly during peak morning activity (0600–1000 h). For each individual captured: species, body mass, sex, age, weight, tail length, and wing length were recorded. Birds were banded with a federally issued band and systematically checked for ticks on the head and body before release, with the head and neck prioritized (Marsot et al., 2012). All attached ticks were removed with forceps and placed in 100% ethanol. Bird species abundance was estimated using distance sampling via point transect surveys (Buckland et al., 2009; Thomas et al., 2010). The surveys were conducted at 10 randomly selected sites located greater than 200 m apart across the island, following recommended practices (Huff et al., 2000). To perform the count, a stationary observer recorded the distance of all avian individuals seen or heard within a 100 m radius for a fixed 3 min interval. To adjust the counts for species detectability, the distribution of detection distances was used to estimate the number of birds present but not detected (see Host Density Estimation section) (Fewster et al., 2008; Thomas et al., 2010). Each site was visited three times in 2014, from early June to early August (Sauer et al., 2013), with approximately four weeks between each visit. The surveys were conducted in hours of peak activity from 0630 h to 1000 h.

White-tailed Deer Sampling

An aerial survey was conducted on 19 February 2015, when conditions provided for 100% snow cover to aid in the detection of white-tailed deer. The aerial survey was accomplished by flying a helicopter at low altitude (approximately 60 m) and at low speed (approximately 40 knots or 46 mph) to allow for good visual identification of deer. Transects of 16 km by 0.2 km wide were plotted on topographic and aerial photographic maps to permit visual identification of features during flights. The pilot was guided by compass and GPS to begin and end survey transects.

Host-Seeking I. scapularis Nymphal Density Sampling

In 2015, host-seeking I. scapularis nymphs were collected at the same three mammal trapping sites (NP, SP, and ML). These collected host-seeking nymphs in 2015 were mostly derived from larvae that fed on the local host community in 2014. Nymphs were collected by dragging a 1 m2 white corduroy sheet between each 10 m section of the entire grid (Daniels et al., 2000; Tälleklint-Eisen and Lane, 2000). Nymphs were collected from the cloth, counted, and preserved in 100% ethanol. For both years, each trapping location was sampled once every other week from late May to late August.

Borrelia burgdorferi Infection Estimation

All collected host-derived larvae were identified morphologically to species using taxonomic keys (Clifford et al., 1961; Durden and Keirans, 1996). Larvae identified as I. scapularis and determined to be fully or near fully engorged were used for estimation of transmission probabilities. Genomic DNA was extracted from both host-derived l. scapularis larvae and host-seeking I. scapularis nymphs using the QIAamp 96 DNA QIAcube HT robot and kits (Qiagen, Valencia, CA) following manufacturer’s protocols. Prior to DNA extraction, each tick was individually frozen in liquid nitrogen and homogenized using a sterile pestle in lysis buffer and proteinase K. The tick DNA extracts were screened for the presence of B. burgdorferi s.l. DNA by real-time quantitative PCR using an Applied Biosystems 7500 Real-Time PCR System (Applied Biosystems, Foster City, CA). A primer (16S-F: 5′−GGCGGC ACACTTAACACGTTAG−3′, 16S-R: 5′−GGCGGCACACTTAACACGT TAG−3′) and probe (6FAM-TTCGGTACTAACTTTTAGTTAA-MGBNFQ) combination that targeted the 16S rRNA region of the bacteria was used following a slightly modified protocol from (Barbour et al., 2009) (we added bovine serum albumin to the PCR master mix). This primer-probe combination distinguished the non tick-borne relapsing fever B. burgdorferi sensu lato spirochetes from B. miyamotoi. We expect most or all samples will be B. burgdorferi s.s. because other B. burgdorferi s.l. species are rarely reported from the Northeast U.S., but see (Margos et al., 2014).

Data analyses

Host Density Estimation

We assessed all host-associated parameters in 2014 and estimated the infection prevalence of I. scapularis nymphs in 2015. Mouse densities per hectare at each site for 2014 were estimated using the R package “Spatially Explicit Capture-Recapture” (SECR). This package estimates animal population density with trapping history data by fitting a spatial detection model by numerically maximizing the log likelihood (minimize the negative log likelihood) using the Nelder-Mead method (McKinnon, 1999; Powell, 1999). A Poisson distribution for mouse density was assumed because of the random nature of the spatial distribution of home range centers of mice. A trapping buffer was introduced to balance the edge effects and the computational expense. The detection probability is assumed to be zero when the distance between detector and mouse home range center is larger than the buffer size. Thus, a finite buffer size could introduce bias in density estimates. To ensure this bias was acceptable, we calculated the relative bias (defined as the ratio of bias of density estimate using the given buffer to the estimate with infinite buffer) for the default buffer (100 m) and confirmed the bias is less than 0.01% for all models considered.

Mouse density estimates were based on data from seven trapping sessions conducted every two weeks, each comprising three consecutive nights of trapping. To accurately estimate mouse capture probability; several covariates such as learned responses (recaptures of specific individuals throughout a season due to bait awareness) and trapping session were incorporated in the half-normal detection function parameterized by g0 (the probability of capture when the trap and center of the home range coincide) and (the spatial scale of the detection function). For each site, models were compared using AICc (corrected Akaike’s Information Criterion with small sample adjustment) values and the model with the lowest AICc was selected as the best fit model.

To fully capture the spatial variability in mouse densities across the island, we defined a uniform distribution for mouse density with the upper limit of the range being the maximum estimated density (among the three sites) and the lower limit being the minimum estimated mouse density (among the three sites). A random sample from this uniform distribution was used as the input parameter for a Poisson distribution to generate the possible values of mouse density in the numerical simulation.

Bird densities per hectare were estimated for all resident passerine species. These were estimated from the survey data using DISTANCE Software 6.2 Release 3.10 (Thomas et al., 2010). For all analyses, Conventional Distance Sampling (CDS) analysis engine was used and a log-normal distribution was assumed for density estimates to account for data including zeros or positive values. The single-observer surveys were filtered to include only the first 3-minutes of observation and were stratified by species (Alldredge et al., 2007); recorded observations were categorized into three categories (0 m–25 m, 25 m–50 m, 50 m–75 m). The CDS engine assumes all objects at zero distance are detected and applies the same detection function where the detection probability is a function of the distance from transects to all objects. The detection function and encounter rate were fitted for each estimation of density and abundance per species using various combinations of detection probability distributions (Uniform, Half-normal, Hazard rate, Negative binomial, Negative exponential) and model adjustments (cosine, simple polynomial, hermite polynomial) (Buckland et al., 2009; Thomas et al., 2010). The model with the lowest AICc was chosen as the most parsimonious model. Note that no model averaging was considered because no competing models, AIC < 2, were obtained (Anderson and Burnham, 2002). The resulting point estimate and standard error of avian density for each species were used to calculate the inputs of the log-normal distribution used in the numerical simulation. We omitted from the analyses bird species that were not observed frequently enough during survey counts to generate reliable density estimates (Buckland et al., 2009; Thomas et al., 2010) and/or were infested with ≤1 engorged larva, which did not allow estimation of realized reservoir competence.

The island-wide abundance of white-tailed deer in the summer of 2014 (pre-harvest) was estimated by adding the number of deer harvested during the 2014–2015 hunting season (assumed measured without error) to estimates of the deer population size derived from the aerial survey conducted post-harvest in 2015 (Tefft, 2016). The latter was obtained by multiplying the deer density estimated in the aerial census by the total habitable land area, excluding open freshwater and saltwater ponds (2386 ha) (Clough and Fulk, 1969). Based on Kilpatrick et al. (2001), a visual correction factor between 1.8 and 2.2 was applied to the observed census data to generate the maximum and minimum values in a uniform distribution.

Larval Burden Estimation on Reservoir Hosts

For white-footed mice and each avian host species, larval burden was estimated by collecting and counting larval ticks from captured individuals. This estimate was restricted to the peak larval period (July 15-August 28), when 96% of all the larvae were counted.

Back-estimation of Larval Burdens on White-tailed Deer

For deer, larval burden could not be directly estimated because of the lack of overlap between the questing larval peak in July and August and white-tailed deer hunting season, which starts in September. Larval burden was thus estimated from back-calculations of the proportion of larvae that needed to have fed on deer in 2014 to explain the NIP observed in 2015. With a value for the proportion of larvae feeding on deer, the average larval burden for deer was then calculated by dividing the estimated total larval abundance on deer to the estimated deer count. The general framework for the back-calculation is illustrated in Fig. 1; detailed procedures and equations are provided in the Supplementary Information.

Fig. 1.

Algorithm used to back-estimate the proportion of larvae fed by each host group and the average larval burden on white-tailed deer. We considered a simplified host community with white-footed mice, birds (seven species combined), and white-tailed deer. We simulated N = 100,000 engorged larvae sampled from mice and birds with probability r and r, which is proportional to the product of their densities d and average larval burdens b (Equation 1 in Supplementary materials). Larvae (L) acquire infection according to the measured realized reservoir competence (red tick symbols represent infected larvae) and molt to nymphs (N) according to a fixed probability derived from estimated molting success in LoGiudice et al. (2003). We estimated how many (uninfected) nymphs would have to be derived from deer (N) to result in the observed nymphal infection prevalence (NIP; Equation 2 in Supplementary information). The number of engorged larvae feeding on white-tailed deer L and the proportion of larvae feeding on different hosts p were calculated after applying the inverse processes of molting and infection (Equations 3 and 4 in Supplementary information). Finally, the number of deer present in an area producing 100,000 engorged larvae was applied to calculate the average larval burden on white-tailed deer (Equation 5 in Supplementary information). These steps were for 10,000 iterations to derive the distributions of density, average larval burden, and proportion of larvae feeding on different hosts, shown in Fig. 2. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article).

Realized Reservoir Competence Estimation

Realized reservoir competence was calculated for each species as the probability of a host transmitting the pathogen to a feeding vector (LoGiudice et al., 2003). This is an integrated measure of the infection prevalence of the host population and the probability of transmitting the infection to a feeding larva, if the host is infected. Realized reservoir competence for each host individual was calculated as the percentage of B. burgdorferi positive larvae compared to all larvae collected from hosts during the peak larval period (July 15-August 28). Realized reservoir competence for deer was assumed to be zero based on Luttrell et al. (1994) and Telford et al. (1988) (see local sensitivity analysis to evaluate this assumption). To account for the uncertainty in the estimates of larval burdens, densities, and realized reservoir competencies of the mouse and avian hosts, all of which contributed to the calculation of larval burden on deer, we randomly sampled from their distributions to generate a distribution for larval burden on white-tailed deer (see Supplementary Information).

Molting Success Parameter

Molting success parameters were derived from LoGiudice et al. (2003), which is the most detailed assessment available. In this study, LoGiudice et al. (2003) collected engorged larvae as they dropped off hosts during a 72 h period and counted the number of larvae that molted into nymphs. They then defined molting successes as the ratio of the number of nymphs to the number of engorged larvae.

Host-Seeking I. scapularis Nymphal Infection Prevalence (NIP) with B. burgdorferi Estimation

We obtained an island-wide estimate of B. burgdorferi NIP as the proportion of B. burgdorferi positive nymphs out of 570 randomly selected nymphs screened in 2015 (190 nymphs from each of the three trapping sites across the island). We used bootstrapping (10,000 samples with replacement) to estimate the standard deviation and minimum and maximum values for NIP.

Alternative NIP Scenarios

To illustrate how NIP varies under all possible larval feeding apportionments over mouse, avian, and white-tailed deer hosts, we expressed NIP as a function of the proportion of ticks feeding on mice, on avian hosts (combined), and on white-tailed deer (see Supplementary Information). These three proportions sum to 1 and the results are presented as a ternary plot (Fig. 2).

Fig. 2.

Distribution for observed and modeled parameters. Value distributions for observed parameters (Columns A, B and C), fixed values for molting success according to LoGiudice et al. (2003) (Column D) and modeled proportions of larvae on Block Island feeding on white-footed mice, seven avian hosts, and white-tailed deer in 2014 (Column E). The average larval burden on deer was not observed in this study, the distribution shown was inferred from the results of the numerical simulation.

Sensitivity Analysis

We calculated elasticities as a simple measure of (local) sensitivity of model output (average larval burden on deer) to changes in each input variable. We estimated elasticity as the percentage change in the proportion of larvae feeding on white-tailed deer in response to a 1% change in an input variable. The definition of elasticity assumes a linear relationship between the input and the output changes. Thus, the results were used to assess the relative importance of any potential biases in the input variables that may have been introduced by methodological limitations (Fig. 3).

Fig. 3.

Predicted nymphal infection prevalence of Borrelia burgdorferi under different larval feeding scenarios on Block Island, RI (varying the quantities p, p, and p). The corners of the ternary plot represent hypothetical scenarios wherein all larvae feed on a single host type (i.e. either exclusively deer, mice, or birds) while the edges of the ternary plot represent scenarios where larvae feed on only two of the three host types; points inside the triangle represent some combination of all three. The white dot (values indicated in black box) indicates the proportions of larvae estimated to have fed on birds, mice and deer to result in the observed nymphal infection prevalence (NIP) of 30%.

Results

Host Sampling

Small mammal sampling in 2014 resulted in 328 captures of 91 individuals from late May to late August. Of the 328 captures, 99.4% were white-footed mice (P. leucopus); the other 0.6% included a meadow vole (M. pennsylvanicus) and a Norway rat (R. norvegicus). The average larval burden of the mice caught during a trapping session (first occurrence if mouse was recaptured) between July 15 and August 28 was 39.7 larvae per individual.

Bird sampling resulted in 176 captures of 173 unique individuals including 18 different species. Of the 18 avian species caught, seven were included in this study: American Robin (Turdus migratorius), Carolina Wren (Thryothorus ludovicianus), Common Yellowthroat (Geothlypis trichas), Eastern Towhee (Pipilo erythrophthalmus), Gray Catbird (Dumetella carolinensis), Song Sparrow (Melospiza melodia), and Yellow Warbler (Setophaga petechia). The mean larval burden for the captured avian species between July 15 and August 28 ranged from 0 to 31.7 (Supplementary Table 1).

Host Density Estimates

White-footed mouse density estimates varied across the three study sites: 11.99 mice/ha (NP), 14.73 mice/ha (SP), and 7.99 mice/ha (ML), with an average density across the three sites of 11.57 mice/ha (Supplementary Table 2). Bird density estimates varied from11.60 birds/ha for Gray Catbirds (GRCA) to 0.59 birds/ha for Yellow Warblers (YEWA) (Supplementary Table 1). Parameter estimates of the best fit models for mouse and bird density are summarized in Supplementary Table 3. Aerial surveys estimated a total abundance of 460 white-tailed deer on the island (Tefft, 2016). This abundance was multiplied by a 1.8–2.2 correction factor, resulting in an estimated 828–1012 white-tailed deer across Block Island. To estimate pre-harvest deer abundance, we added the 387 deer harvested during the 2014–2015 hunting season (assumed measured without error) for a total abundance of 1215–1399 white-tailed deer. We divided deer total abundance by the 2386 ha of habitable land area to yield a density of 0.51–0.59 white-tailed deer ha−1.

Nymphal Infection Prevalence Estimate

The island-wide estimate for NIP was 30% ± 2% (SD), calculated by pooling all nymphs collected from three study sites and using boot-strapping to estimate the SD.

Distributions for the Proportion of Larvae Feeding on Mice, Birds, and White-tailed deer

The largest proportion (mean ± SD) of larvae determined by back calculation fed on mice (p, 0.44 ± 0.10), followed by white-tailed deer (p, 0.29 ± 0.07), and then birds which in order of contribution included: Carolina Wren (p, 0.17 ± 0.08), Gray Catbird (p, 0.05 ± 0.03), Common Yellowthroat (p, 0.02 ± 0.02), American Robin (p 0.01 ± 0.01), Eastern Towhee (p, 0.01 ± 0.01), House Sparrow (p, 0.01 ± 0.01), and Yellow Warbler (p, 0.00 ± 0.00) (Table 1). The variance of the distribution (note differences in x-axis scale) was highest for white-footed mice, Carolina Wrens, and white-tailed deer and was lowest for American Robins, Eastern Towhees, Song Sparrows, and Yellow Warblers. Intermediate distributions were observed for Common Yellowthroats and Gray Catbirds. The average estimated larval burden on deer - calculated from a fixed-point estimate of deer density for each run of the numerical simulation was determined to be 555.24 ± 258.82 (mean ± SD) larvae (Fig. 1). Of the 10,005 samples, there were 5 cases (< 0.1%) where the total infection prevalence in non-deer hosts was less than the bootstrapped NIP. Those cases were excluded from the analyses because the estimated values of N were negative and violated the biological reasonable constraint that N ⩾ 0.

Table 1

Sampling results of nymphal infection prevalence (NIP), host densities, mean larval burdens on hosts, and realized reservoir competences of hosts from fitted distributions and bootstrapping are summarized. Point values of molting success percentage for white-footed mice (Peromyscus leucopus), all avian species (Turdus migratorius, Thryothorus ludovicianus, Geothlypis trichas, Pipilo erythrophthalmus, Dumetella carolinensis, Melospiza melodia, Melospiza melodia), and white-tailed deer (Odocoileus virginianus) were used in the back-estimating simulation of the proportion of larvae feeding on white-tailed deer. Back-estimated results of mean larval burden on deer and proportion of larvae feeding on different hosts are calculated from 10,000 samples.

	Mean	Minimum	Maximum	SD
NIP	0.30	0.23	0.38	0.02
Density
Peromyscus leucopus	11.37	1.00	28.00	3.86
Turdus migratorius	0.71	0.11	3.08	0.32
Thryothorus ludovicianus	5.60	0.94	26.68	2.58
Geothlypis trichas	1.56	0.03	30.63	1.87
Pipilo erythrophthalmus	4.14	0.96	14.23	1.59
Dumetella carolinensis	11.64	2.45	42.23	4.68
Melospiza melodia	1.70	0.38	5.44	0.61
Setophaga petechia	0.60	0.03	7.13	0.51
Odocoileus virginianus	0.55	0.51	0.59	0.02
Burden
Peromyscus leucopus	39.67	23.68	60.07	4.94
Turdus migratorius	8.49	0.08	44.00	6.69
Thryothorus ludovicianus	31.65	8.43	92.71	12.68
Geothlypis trichas	9.91	3.21	21.21	2.62
Pipilo erythrophthalmus	2.34	0.17	6.33	0.94
Dumetella carolinensis	3.90	1.42	7.31	0.79
Melospiza melodia	4.87	0.80	13.27	1.64
Setophaga petechia	0.33	0	1.00	0.16
Odocoileus virginianus	555.24	0.35	3904.76	258.82
Realized reservoir competence
Peromyscus leucopus	0.49	0.3	0.65	0.05
Turdus migratorius	0.64	0	1.00	0.20
Thryothorus ludovicianus	0.59	0.10	0.87	0.13
Geothlypis trichas	0.35	0.08	0.64	0.09
Pipilo erythrophthalmus	0.38	0	1.00	0.21
Dumetella carolinensis	0.06	0	0.24	0.04
Melospiza melodia	0.58	0.05	1.00	0.16
Setophaga petechia	0.08	0	0.17	0.06
Odocoileus virginianus	0	0	0	0
Molting success
Peromyscus leucopus	0.42	0.42	0.42	0
Turdus migratorius	0.34	0.34	0.34	0
Thryothorus ludovicianus	0.34	0.34	0.34	0
Geothlypis trichas	0.34	0.34	0.34	0
Pipilo erythrophthalmus	0.34	0.34	0.34	0
Dumetella carolinensis	0.34	0.34	0.34	0
Melospiza melodia	0.34	0.34	0.34	0
Setophaga petechia	0.34	0.34	0.34	0
Odocoileus virginianus	0.56	0.56	0.56	0
Proportion (fed larvae)
Peromyscus leucopus	0.44	0.05	0.79	0.10
Turdus migratorius	0.01	0	0.09	0.01
Thryothorus ludovicianus	0.17	0.02	0.61	0.08
Geothlypis trichas	0.02	0	0.27	0.02
Pipilo erythrophthalmus	0.01	0	0.08	0.01
Dumetella carolinensis	0.05	0.01	0.34	0.03
Melospiza melodia	0.01	0	0.07	0.01
Setophaga petechia	0	0	0	0
Odocoileus virginianus	0.29	0	0.51	0.07

Alternative NIP Scenarios

The variable NIP scenarios were visualized within the ternary plot; the observed scenario on Block Island in this study is denoted by the black circle (Fig. 2). The remainder of the plot represents hypothetical scenarios with varying proportion of larva feeding on the three hosts. The corners of the triangle represent extremes where the community consists of only one host species, while the edges of the triangle represent hypothetical communities consisting of just two out of the three species (Fig. 2). Therefore, if there were only mice, or only birds, or a combination of mice and birds, the predicted NIP would be similar; birds and mice play similar roles with respect to the determination of NIP. The upper corner of the triangle represents a community consisting only of white-tailed deer (with NIP equal to 0), which emphasizes the strong influence of deer in reducing nymphal infection within this simplified three-host community.

Sensitivity Analysis

Local sensitivity analysis demonstrates that the proportion of larvae feeding on deer is sensitive to changes in realized reservoir competence of white-footed mice and the observed island-wide NIP. The proportion of larvae feeding on deer was moderately sensitive to changes in molting success and the realized reservoir competence of Carolina These are therefore highly influential parameters that should be carefully measured Wrens; the remaining parameters had little influence on the proportion of larvae feeding on deer when their uncertainties were considered in the numerical simulations (Table 2). Using elasticity analysis, we assessed the relative importance of potential biases introduced by some methodological limitations of the study. There were three potential biases in our study: (1) avian density may have been underestimated by including a point count survey in August, when male songbirds reduce singing frequency (Ralph et al., 1995); (2) larval burdens on all hosts may have been underestimated by counting engorged larvae while still attached to the host (Schmidt et al., 1999; Brunner and Ostfeld, 2008); and (3) limited studies indicate measuring infection in engorged larvae rather than molted nymphs may have resulted in an underestimation of realized reservoir competence because the molted nymphs have higher infection prevalence than the engorged larvae (Jacquet et al., 2017). Sensitivity analyses indicated that the first two factors had a relatively small impact on the estimation of the proportion of larvae feeding on deer, 0.0346% and 0.0001% respectively. An underestimation of realized reservoir competence would however lead to a larger change in the calculated proportion of larvae feeding on deer. For example, if the actual realized reservoir competence were 10% higher for each species, the proportion of larvae feeding on white-tailed deer would increase to 35% (i.e. a 6% increase).

Table 2

Local sensitivity indices of the proportion of larvae feeding on white-tailed deer on all model input parameters calculated by elasticity analysis. Note that average larval burden on deer is not an input parameter in our model.

	Peromyscus leucopus	Turdus migratorius	Thryothorus ludovicianus	Geothlypis trichas	Pipilo erythrophthalmus	Dumetella carolinensis	Melospiza melodia	Setophaga petechia	Odocoileus virginianus	Ixodes scapularis nymphs
Density	0.0345	0.004	0.0686	−0.0111	−0.0056	−0.0929	0.0028	−0.0004	0	NA
Average larval burden	0.0345	0.004	0.0686	−0.0111	−0.0056	−0.0929	0.0028	−0.0004	NA	NA
Realized reservoir competence	1.2513	0.0180	0.4820	0.0249	0.0168	0.0125	0.0221	0	0.0759	NA
Molting success	0.4692	0.0097	0.2369	0.0036	0.0035	−0.0500	0.0106	−0.0002	−0.6834	NA
NIP										−1.9036

Discussion

We aimed to investigate whether reduced infection in I. scapularis nymphs due to larval feeding on incompetent white-tailed deer and other hosts, could account for the lower than expected I. scapularis infection prevalence of B. burgdorferi in a species-poor host community. The reduced diversity of this island community allowed us to measure most key host parameters (i.e. density, tick burden, realized reservoir competence) - only assuming molting success from the literature, in order to estimate the parameter of interest, the larval burden on whitetailed deer. We back-estimated that about a third (29%) of larval ticks must have fed on deer to yield the observed 30% NIP, which is similar to our previous NIP estimates of 28% in 2010 and 24% in 2011 (States et al., 2014), and similar or lower than most other studies on mainland communities with higher species diversity (Allan et al., 2003; DiukWasser et al., 2012). These findings illustrate the potential dual role of white-tailed deer in the eco-epidemiology of B. burgdorferi. At high densities, deer may switch from an amplifying to a dilution host if their role in reducing tick infection overrides their role as amplifiers of the tick population, as the latter can become limited by other abiotic or biotic factors affecting tick survival (Lindsay et al., 1995, Lindsay et al., 1998). This study provides a proof of concept of the potential role of deer as dilution host, which needs to be further investigated in other ecological settings.

While the role of white-tailed deer as the key host for adult ticks has been widely recognized, their role as hosts for larval ticks has been poorly studied empirically. Only two studies obtained such data in the United States using special deer harvest permits around peak larval season (Telford et al., 1988; LoGiudice et al., 2003). These studies estimated larval burdens of 239 ± 99 larvae per white-tailed deer (LoGiudice et al., 2003) and 341.9 ± 115.5 larvae per white-tailed deer (Telford et al., 1988). Similarly we identified only two studies in Europe reporting larval burdens collected throughout the body from roe deer shot around the peak larval season in Sweden, reporting an average of 276 (range, 84–658) larval I. ricinus (n = 12 deer) (Jaenson and Tälleklint-Eisen, 1992) and an average of 265 larvae (n = 37 deer) (Tälleklint and Jaenson, 1997). The higher average larval tick estimate obtained in the current study (555.24 ± 258.82) is within the range of previous studies. Our higher value is consistent with potential underestimation in published studies due to their sampling partially off-season or not capturing the entire body burden.

An important assumption of the dilution effect hypothesis is that the most competent hosts are also the most resilient to extirpation under anthropogenic habitat modification (Ostfeld and Keesing, 2000; McCallum, 2015). This does not apply in the case of white-tailed deer, which are both reservoir incompetent and highly abundant (or aggregated) in human-modified habitats because of reduced vulnerability to hunter harvest and increased amounts of accessibility to forage (McAninch et al., 1993; Woolf and Roseberry, 1998; Williams et al., 2008, 2013). The dual role of deer in human-modified landscapes may explain the lack of consistent findings in studies of the association between habitat fragmentation and entomological risk with both positive (Allan et al., 2003; Brownstein et al., 2005) and non-significant (LoGiudice et al., 2008; Zolnik et al., 2015) findings. Because the density of white-footed mice and white-tailed deer both increase in more fragmented landscapes (McAninch et al., 1993; Nupp and Swihart, 1998), we caution against the use of forest fragmentation as a ‘proxy’ for biodiversity in studies of the association between biodiversity and human infection risk.

The small mammal community on Block Island, RI is markedly less diverse than those in comparable mainland settings (States et al., 2014), which facilitated estimation of the most relevant host parameters. White-footed mice, the most competent B. burgdorferi host in this setting, dominates the small mammal community on Block Island with a mean density of 11.57 mice per hectare in 2014 and contributing to 44% of fed larvae to the tick community; a similar contribution to that reported in more biodiverse mainland communities (LoGiudice et al., 2003; Tsao et al., 2004; Brisson and Dykhuizen, 2006). The similar NIP in this study (30% ± SD 2%) and in a previous study in this community (24% ± SE 27%) (States et al., 2014) compared to mainland communities (Allan et al., 2003; Diuk-Wasser et al., 2012) suggests that additional dilution hosts contribute to the lower than expected infection prevalence. Using comparable approaches, previous studies implicated flying squirrels (LoGiudice et al., 2008) or shrews (Brisson and Dykhuizen, 2006) as ‘missing’ dilution hosts. The absence of these mammalian hosts on Block Island allowed for the assessment of the (neglected) role of white-tailed deer in reducing nymphal infection prevalence. Deer density on Block Island was estimated to be 0.51–0.59 deer ha−1, which is similar to other intermediate density communities on the mainland (Adams et al., 2009), emphasizing the potential importance of deer as a larval host on the mainland as well as on island settings.

The extremely low mammalian community diversity also permitted the assessment of the role birds play as B. burgdorferi hosts, which is less emphasized than the role of mammals in other studies. It has been suggested that birds play a significant role in dispersal of infected larval ticks but are not considered as relevant amplification or dilution hosts (Brinkerhoff et al., 2011). However, limited studies have identified specific avian species as reservoir hosts for B. burgdorferi (LoGiudice et al., 2003, 2008; Hamer et al., 2011; Newman et al., 2015). We estimated birds contributed a combined total of 27% of fed larvae to the environment. In particular Carolina Wrens (CARW) contributed the most B. burgdorferi infected ticks (5 infected ticks per individual from field data; 79% of the avian-infected larvae in numerical simulations) and showed the highest mean larval burden (31.7 ticks per individual) compared to the other six avian species.

The methods used to collect data in this study may have caused inaccuracies in the estimations of avian density, larval burdens on mammals and birds, and the realized reservoir competence of hosts. Local sensitivity analysis showed that the proportion of larvae feeding on deer is highly sensitive to the change in realized reservoir competence of white-footed mice and NIP. These are therefore highly influential parameters that should be carefully measured. Moderately sensitive – or less influential parameters, are changes in molting success and the realized reservoir competence of Carolina Wrens. Limited literature indicates a potential underestimation of the calculated realized reservoir competence of white-footed mice, while there is no reason to assume any systematic biases in NIP. Therefore, our numerical results indicating that 29% of larval ticks should have fed on white-tailed deer would be the lower bound when these corrections are considered. An additional limitation of this study was that some avian species were sampled in point transect surveys but were not captured, therefore they were not included in the analyses. Many of these omitted species are not commonly associated with tick-borne pathogens because of the limited time they spend foraging on the ground. However, the Blue Jay (Cyanocitta cristata), Northern Cardinal (Cardinalis cardinalis), and Ring-necked Pheasant have shown moderate burdens and reservoir competence in previous studies (Anderson and Magnarelli, 1984; Kurtenbach et al., 1998; Hoodless et al., 2002; Hamer et al., 2012; Ginsberg et al., 2005). The impact of the absence of these species from our analyses cannot be directly estimated. Finally, our study did not account for heterogeneities across the island in the various estimates because of the limited availability of island-wide estimates of deer density. The multiple scales relevant to the various I. scapularis hosts is one of the challenges in testing hypotheses about community-wide effects on B. burgdorferi entomological risk (Killilea et al., 2008).

While our study has identified a potentially important role of white-tailed deer in reducing infection of immature ticks, an integrated, empirically-informed, assessment of the dual role of deer as a vector amplifying and pathogen dilution host is still pending. A modeling study concluded that either amplification or dilution may occur with the outcome depending on the precise mechanisms of competition, host contact rates with ticks, and acquired host resistance to ticks (Ogden and Tsao, 2009). Levi et al. (2016) analyzed long-term data from the Cary Institute to quantify dilution and amplification by various host species. However, parameters for white-tailed deer density and their immature tick burdens were based on a single fixed estimate for tick burden, therefore the effect of spatial and temporal variation in deer abundance and aggregation on pathogen transmission was not assessed. Empirical measurements of deer density or aggregation and larval burdens are essential to fully characterize the dual role of white-tailed deer in Lyme disease epidemiology (Kilpatrick et al., 2017).

Conclusions

Our study identified an influential role of white-tailed deer in reducing nymphal tick infection prevalence (NIP) and a potential role as dilution hosts if the reduction in NIP outweighs the role of deer as tick population amplifiers. We emphasized the importance of assessing larval feeding on deer, which has only been measured in two previous studies in the United States to a very limited extent (Telford et al., 1988; LoGiudice et al., 2003). Reduction of unfed nymphal infection by white-tailed deer at high deer densities could reverse the negative relationship between biodiversity and Lyme disease risk. Importantly, it would also reduce the efficacy of white-tailed deer population reduction efforts and thus should be considered in models predicting effective deer target densities for disease risk reduction.